H Berjamin
email: berjamin@lma.cnrs-mrs.fr

Bruno Lombard
email: lombard@lma.cnrs-mrs.fr

G Chiavassa

N Favrie

A finite-volume approach to 1D nonlinear elastic waves: Application to slow dynamics

A numerical method for longitudinal wave propagation in nonlinear elastic solids is presented. Here, we consider polynomial stress-strain relationships, which are widely used in nondestructive evaluation. The large-strain and infinitesimal-strain constitutive laws deduced from Murnaghan's law are detailed, and polynomial expressions are obtained. The Lagrangian equations of motion yield a hyperbolic system of conservation laws. The latter is solved numerically using a finite-volume method with flux limiters based on Roe linearization. The method is tested on the Riemann problem, which yields nonsmooth solutions. The method is then applied to a continuum model with one scalar internal variable, accounting for the softening of the material (slow dynamics).

Introduction

Rocks and concrete are known to behave nonlinearly when vibrating longitudinally, even at very low amplitudes [START_REF] Guyer | Nonlinear mesoscopic elasticity: evidence for a new class of materials[END_REF][START_REF] Tencate | Slow dynamics of earth materials: an experimental overview[END_REF]. Firstly, dynamic acoustoelastic testing (DAET) [START_REF] Renaud | Revealing highly complex elastic nonlinear (anelastic) behavior of Earth materials applying a new probe: Dynamic acoustoelastic testing[END_REF][START_REF] Rivière | Pump and probe waves in dynamic acoustoelasticity: comprehensive description and comparison with nonlinear elastic theories[END_REF] reveals that the speed of sound measured locally decreases with time, and recovers its initial value after the excitation is stopped. This softening occurs over a time scale larger than the period of the dynamic loading, which highlights the phenomenon of slow dynamics. Secondly, the evolution of this speed with respect to the strain presents an hysteretical behavior. Lastly, all these phenomena are enhanced when the forcing amplitude is increased. Such observations are not compatible with linear elastodynamics, where the speed of sound is a constant.

Polynomial nonlinear stress-strain relationships are widely used in nondestructive testing [START_REF] Van Den Abeele | Elastic pulsed wave propagation in media with second-or higher-order nonlinearity. Part I. Theoretical framework[END_REF][START_REF] Payan | Applying nonlinear resonant ultrasound spectroscopy to improving thermal damage assessment in concrete[END_REF]. However, nonlinear elasticity is not sufficient to represent the softening and the hysteresis revealed by DAET. Several models can be found in the literature to reproduce these phenomena (see e.g. [START_REF] Meurer | Wave propagation in nonlinear and hysteretic media-a numerical study[END_REF][START_REF] Li | Pseudo-spectral simulation of 1D nonlinear propagation in heterogeneous elastic media[END_REF]).

The soft-ratchet model by Vakhnenko et al. [START_REF] Vakhnenko | Strain-induced kinetics of intergrain defects as the mechanism of slow dynamics in the nonlinear resonant response of humid sandstone bars[END_REF][START_REF] Favrie | Fast and slow dynamics in a nonlinear elastic bar excited by longitudinal vibrations[END_REF] consists in introducing a scalar variable g to describe the softening of the material. A similar model with refinements was proposed by Lyakhovsky and coauthors in a series of papers [START_REF] Lyakhovsky | Distributed damage, faulting, and friction[END_REF][START_REF] Lyakhovsky | Non-linear elastic behaviour of damaged rocks[END_REF]. Recently, a 3D model of continuum has been developed by the authors [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internal-variable model[END_REF] in the framework of continuum thermodynamics with internal variables of state [START_REF] Maugin | Thermodynamics with internal variables[END_REF][START_REF] Maugin | The saga of internal variables of state in continuum thermo-mechanics (1893-2013)[END_REF]. This model, which generalizes the soft-ratchet model to 3D geometries and fixes thermodynamical issues, will be used along the present paper.

Regardless the fact that slow dynamics is taken into account or not, the equations of the model appear as a nonlinear hyperbolic system of conservation laws with relaxation, the solutions of which may be discontinuous (shock waves). The numerical computation of nonsmooth solutions requires a particular care. Here, dedicated finite-volume methods [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF][START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF][START_REF] Velasco-Segura | A finite volume approach for the simulation of nonlinear dissipative acoustic wave propagation[END_REF] are used to compute such solutions in a non-oscillatory manner. A 1D scheme based on the Roe linearization and flux limiters is adapted to the present system of equations. Due to the separation of timescales between the wave propagation and the softening of the material, the coupling between nonlinear elastodynamics and the slow dynamics is straightforward.

The article is organized as follows. For pedagogical purposes, the case of nonlinear elastic solids without slow dynamics is first considered (sections 2 and 3). Section 2 derives the equations of nonlinear elastodynamics, and the case of Murnaghan hyperelastic material [START_REF] Murnaghan | Finite deformations of an elastic solid[END_REF] is tackled. The use of this law is discussed, in particular the way it relates to polynomial stressstrain relationships [START_REF] Ostrovsky | Dynamic nonlinear elasticity in geomaterials[END_REF][START_REF] Payan | Determination of third order elastic constants in a complex solid applying coda wave interferometry[END_REF]. The numerical method is presented in section 3, including the construction of a Roe matrix. In section 3.4, the numerical method is validated with the analytical solution of the Riemann problem of nonlinear elastodynamics [START_REF] Berjamin | Analytical solution to 1D nonlinear elastodynamics with general constitutive laws[END_REF]. Section 4 introduces the modifications of the equations to account for the slow dynamics, as well as the modifications of the numerical method. Then, the propagation of a sinusoidal wave in the material is addressed (section 4.4). The results are in qualitative agreement with experimental observations.

Nonlinear elastodynamics

2.1. Governing equations Lagrangian hyperelasticity.. Let us consider an homogeneous continuum in which no heat transfer occurs. Furthermore, self-gravitation is neglected. A particle initially located at some position x 0 of the reference configuration moves to a position x t of the current configuration. The deformation gradient is a second-order tensor defined by (see e.g. [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF][START_REF] Holzapfel | Nonlinear Solid Mechanics: A Continuum Approach for Engineering[END_REF][START_REF] Drumheller | Introduction to wave propagation in nonlinear fluids and solids[END_REF][START_REF] Norris | Finite-amplitude waves in solids[END_REF])

F = grad x t = I + grad u , (1)
where u = x tx 0 denotes the displacement field and grad is the gradient with respect to the material coordinates x 0 (Lagrangian gradient). In the reference configuration, the deformation gradient ( 1) is equal to the metric tensor I. Here, the Euclidean space is described by an orthonormal basis (e 1 , e 2 , e 3 ) and a Cartesian coordinate system (O, x, y, z). In this case, the matrix of the coordinates of I is the identity matrix. Here, the Lagrangian representation of motion is used. The material derivative of the deformation gradient satisfies

Ḟ = grad v , (2)
where v(x 0 , t) is the velocity field. The conservation of mass implies ρ 0 /ρ = det(F ), where ρ denotes the mass density in the deformed configuration, and ρ 0 denotes the mass density in the reference configuration. The motion is also driven by the conservation of momentum

ρ 0 v = div det(F ) σ • F -⊤ + f v , (3)
where div denotes the divergence with respect to the material coordinates.

The tensor σ = σ ⊤ is the Cauchy stress tensor, and f v is an external volume force applied to the material. In hyperelasticity, the only variables of state are the specific entropy η and a strain tensor. Moreover, the dissipation in the material is zero, i.e. the thermodynamic process is isentropic. Here, the deformation of the material is represented by the Green-Lagrange strain tensor, but other choices are possible. The Green-Lagrange strain tensor is defined as

E = 1 2 (C -I)
where C = F ⊤ • F is the right Cauchy-Green strain tensor, or equivalently as a function of the displacement gradient tensor:

E = 1 2 grad u + grad ⊤ u + grad ⊤ u • grad u . (4)
The internal energy per unit of reference volume is W (E), where W is the strain energy density function. Under these assumptions, the expression of the Cauchy stress tensor is

σ = 1 det(F ) F • ∂W ∂E • F ⊤ . (5)
Thus, the first Piola-Kirchhoff tensor in (3) satisfies

det(F ) σ • F -⊤ = F • ∂W ∂E . (6)
In the isotropic case, the strain energy is a function of the invariants of E defined by

E I = tr(E) , E II = 1 2 tr(E) 2 -tr(E 2 ) , E III = det(E) . (7)
The derivative of the strain energy with respect to the strain tensor writes

∂W ∂E = α 0 I + α 1 E + α 2 E 2 , (8)
where

α 0 = ∂W ∂E I + E I ∂W ∂E II + E II ∂W ∂E III , α 1 = - ∂W ∂E II -E I ∂W ∂E III , α 2 = ∂W ∂E III , (9)
are functions of the invariants (E I , E II , E III ).

Longitudinal plane waves.. We make the assumption that the displacement field u has no component along e 2 and e 3 . Moreover, its component u along e 1 is independent on y and z. Therefore, the displacement gradient writes grad u = ε (e 1 ⊗ e 1 ), where ε = ∂ x u > -1 is the axial component of the displacement gradient. The invariants [START_REF] Meurer | Wave propagation in nonlinear and hysteretic media-a numerical study[END_REF] of the Green-Lagrange strain tensor are E I = ε + 1 2 ε 2 and E II = 0 = E III . Thus, the strain energy W is now a function of ε only. In the longitudinal case, the 11 -coordinate of the first Piola-Kirchhoff tensor ( 6) is equal to the coordinate σ of the Cauchy stress tensor. Equations ( 8)-( 9) and the expression of the invariants provide the stress-strain relationship

σ = (1 + ε) α 0 + E I α 1 + E I 2 α 2 , = ∂E I ∂ε ∂W ∂E I , = W ′ (ε) , (10)
where the ′ denotes the total derivative with respect to ε.

Constitutive laws.. The Murnaghan model of hyperelasticity [START_REF] Murnaghan | Finite deformations of an elastic solid[END_REF] is widely used in the communities of geophysics and nondestructive testing [START_REF] Mccall | Theoretical study of nonlinear elastic wave propagation[END_REF][START_REF] Johnson | Nonlinear elasticity and stressinduced anisotropy in rock[END_REF][START_REF] Ostrovsky | Dynamic nonlinear elasticity in geomaterials[END_REF][START_REF] Payan | Determination of third order elastic constants in a complex solid applying coda wave interferometry[END_REF]. Its strain energy density function is

W = λ + 2µ 2 E I 2 -2µ E II + l + 2m 3 E I 3 -2m E I E II + n E III , (11)
where (λ, µ) are the Lam parameters and (l, m, n) are the Murnaghan coefficients. Sometimes, Landau's law with parameters (A, B, C) is used instead, and the relationship with Murnaghan's law [START_REF] Lyakhovsky | Distributed damage, faulting, and friction[END_REF] is specified in [START_REF] Norris | Finite-amplitude waves in solids[END_REF]. With the above expression of the strain energy [START_REF] Lyakhovsky | Distributed damage, faulting, and friction[END_REF], the axial component of the Cauchy stress [START_REF] Favrie | Fast and slow dynamics in a nonlinear elastic bar excited by longitudinal vibrations[END_REF] is

σ = M 0 ε 1 + 3 2 + ϑ ε + 1 2 + 2ϑ ε 2 + 5ϑ 4 ε 3 + ϑ 4 ε 4 , (12)
where M 0 = λ + 2µ > 0 is the elastic modulus, and ϑ = (l + 2m)/(λ + 2µ). If the Murnaghan coefficients equal zero (ϑ = 0), only geometric nonlinearities remain, and the Saint Venant-Kirchhoff model of hyperelasticity is recovered. The later reduces to the classical Hooke's law in the case of infinitesimal strain ε ≃ 0. When geometric nonlinearities are neglected (e.g., when the Murnaghan coefficients are very large), the Green-Lagrange strain tensor is linearized with respect to the components of grad u, so that the strain tensor ( 4) is replaced by the infinitesimal strain tensor: E ≃ ε. Doing this, the first invariant of the strain tensor is E I ≃ ε in the longitudinal case. Moreover, the first Piola-Kirchhoff stress is linearized with respect to the components of grad u as well, i.e. F • ∂W/∂E ≃ ∂W/∂ε in [START_REF] Payan | Applying nonlinear resonant ultrasound spectroscopy to improving thermal damage assessment in concrete[END_REF]. Doing this, the longitudinal constitutive law [START_REF] Favrie | Fast and slow dynamics in a nonlinear elastic bar excited by longitudinal vibrations[END_REF] reduces to σ = α 0 + α 1 ε + α 2 ε 2 , and ( 12) becomes

σ = M 0 ε (1 + ϑε) . (13)
This constitutive law corresponds to a quadratic polynomial

σ = M 0 ε (1 -βε), with β = -ϑ.
A cubic polynomial constitutive law [START_REF] Meurer | Wave propagation in nonlinear and hysteretic media-a numerical study[END_REF][START_REF] Li | Pseudo-spectral simulation of 1D nonlinear propagation in heterogeneous elastic media[END_REF][START_REF] Favrie | Fast and slow dynamics in a nonlinear elastic bar excited by longitudinal vibrations[END_REF]]

σ = M 0 ε 1 -βε -δε 2 (14)
is also widely used to describe nonlinear elasticity in solids. If the geometric nonlinearities are negligible [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internal-variable model[END_REF], then the parameters β = -ϑ and δ = 0 correspond to Murnaghan's law. If the geometric nonlinearities are taken into account [START_REF] Lyakhovsky | Non-linear elastic behaviour of damaged rocks[END_REF], then the choice -β = 3 2 + ϑ and -δ = 1 2 + 2ϑ makes the two models coincide up to the third order. Contrary to [START_REF] Lyakhovsky | Non-linear elastic behaviour of damaged rocks[END_REF], orders 4 and 5 are not taken into account in [START_REF] Maugin | Thermodynamics with internal variables[END_REF], which shows that both models are not equivalent.

The major difference between the stress-strain relationships ( 12) and ( 14) is the number of independent parameters. In [START_REF] Maugin | Thermodynamics with internal variables[END_REF], the cubic term can be set independently of the quadratic term. In [START_REF] Lyakhovsky | Non-linear elastic behaviour of damaged rocks[END_REF], the cubic term is not independent on the quadratic term: if ( 12) and ( 14) are assimilated, then δ = 5 2 + 2β. However, experimental evidence shows that δ is larger than β by several orders of magnitude, so that a constitutive law of the Murnaghan type is not sufficient to represent accurately elastic nonlinearity in geomaterials [START_REF] Renaud | Revealing highly complex elastic nonlinear (anelastic) behavior of Earth materials applying a new probe: Dynamic acoustoelastic testing[END_REF]. From now on, the cubic polynomial law ( 14) is used.

Hyperbolic system of conservation laws

The equations of motion ( 2)-( 3) write as a non-homogeneous system of conservation laws with respect to the variables q = (ε, v) ⊤ , where v is the particle velocity component along e 1 :

∂ t q + ∂ x f(q) = s . ( 15
)
The flux function is f(q) = -(v, σ(ε)/ρ 0 ) ⊤ , and the forcing is s = (0, f v /ρ 0 ) ⊤ , where f v is the volume force component along e 1 . The Jacobian matrix of the flux is

f ′ (q) = - 0 1 σ ′ (ε)/ρ 0 0 , (16)
where σ ′ = W ′′ is the derivative of σ with respect to ε. The eigenvalues of f ′ (q) are {-c(q), c(q)}, where c(q) = σ ′ (ε)/ρ 0 is the speed of sound. The system of conservation laws is strictly hyperbolic if σ ′ (ε) > 0, i.e. over a domain where the strain energy W is a strictly convex function of ε. Some properties of the system (15) are listed below without proof. Interested readers are referred to standard textbooks for more details about hyperbolic systems [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF][START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF].

In the case of the cubic law ( 14), strict hyperbolicity is ensured if [START_REF] Berjamin | Analytical solution to 1D nonlinear elastodynamics with general constitutive laws[END_REF] ε ∈ 1

β -β 2 + 3δ , 1
β + β 2 + 3δ . ( 17
)
If β = δ = 0, then the polynomial law [START_REF] Maugin | Thermodynamics with internal variables[END_REF] amounts to Hooke's law, and the characteristic fields are linearly degenerate. If β = 0 and δ = 0, then the polynomial law is quadratic, and the stress is either a strictly convex or a strictly concave function of the strain. The characteristic fields are genuinely nonlinear. If δ = 0, then the polynomial law ( 14) is neither convex nor concave. Indeed, an inflexion point is located at ε 0 = -β/3δ. The characteristic fields are neither genuinely nonlinear nor linearly degenerate (i.e., they are nongenuinely nonlinear). In the case of Murnaghan's law [START_REF] Lyakhovsky | Non-linear elastic behaviour of damaged rocks[END_REF], a similar analysis is carried out in the appendix.

Finite volumes with flux limiters

Conservative scheme

In the examples presented later on, the physical domain is unbounded. We consider a finite numerical domain [x 0 , x N ]. It is discretized using a regular grid in space with step ∆x = (x N -x 0 )/N. Also, a variable time step ∆t = t n+1 -t n is introduced. Therefore, q(x i , t n ) denotes the solution to (15) at the abscissa x i = x 0 + i ∆x and at the time t n . The volume force f v is assumed to be a point load f v = δ s (x)ϕ(t), where δ s (x) = δ(x -x s ) is the Dirac delta located at the abscissa x = x s , and ϕ(t) is the source signal.

The non-homogeneous system of conservation laws ( 15) is integrated explicitly:

q n+1 i = q n i - ∆t ∆x (f n i+1/2 -f n i-1/2 ) + ∆t ∆x s n i , (18)
where q n i ≃ q(x i , t n ) approximates the solution at the grid nodes, and s n i = s(x i , t n ). The numerical flux f n i+1/2 of the flux-limiter method is specified later on. This scheme is stable under the classical CFL condition

κ = ∆t ∆x c n max 1 , ( 19
)
where κ is the Courant number, and c n max is the maximum sound speed that is encountered at time t n . If σ is convex or concave, then the maximum sound speed at time t n is

c n max = max 0 i N c(q n i ) . (20)
If the constitutive law is neither convex nor concave, sound speeds larger than (20) may be reached (see e.g. section 16.1 in [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF]). The local maximum sound speed is obtained by maximizing c(q) for q between q n i and q n i+1 . Then, one obtains

c n max = max 0 i<N 1 ρ 0 max ε∈D i σ ′ (ε) , (21)
where D i is the interval with bounds ε n i and ε n i+1 . Finally, the method has a variable time step satisfying ∆t = κ ∆x/c n max . Since the flux function f is nonlinear with respect to q, an initial-value problem (or Cauchy problem) of the homogeneous system defined by the data q(x, 0) at the time t = 0 can have several weak solutions (i.e. solutions of an integral form of ( 15)). It is not straightforward for a numerical method to converge towards the correct weak solution. In particular, the numerical fluxes f n i+1/2 in (44) must be computed carefully. Given that the eigenvalues ±c(q) of ( 16) have constant sign, no transsonic rarefaction can occur. Moreover, no slow-moving shock can occur either, due to the monotonicity of σ over the domain of hyperbolicity [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF]. Therefore, Roe linearization can be used to construct an accurate numerical scheme.

Roe linearization

The nonlinear flux f(q) is approximated locally by the linear flux A i+1/2 q. A Roe matrix A i+1/2 approximates the Jacobian f ′ (q n i+1/2 ) at the midpoint of [x i , x i+1 ] and the time t n , in such a way that 1. A i+1/2 is diagonalizable with real eigenvalues;

2. A i+1/2 q n i+1 -q n i = f(q n i+1 ) -f(q n i ) ; 3. lim q n i+1 →q n i A i+1/2 = f ′ (q n i ) .
According to section 15.3.2 of [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF], such a matrix may be expressed by the formula

A i+1/2 = 1 0 f ′ q n i + ξ (q n i+1 -q n i ) dξ , (22)
which ensures that the properties 2. and 3. are satisfied.

In the case of the system [START_REF] Maugin | The saga of internal variables of state in continuum thermo-mechanics (1893-2013)[END_REF], which corresponds to exercise 15.1.(a) p. 349 of [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF], the following matrix is obtained:

A i+1/2 = 0 -1 -a i+1/2 0 , (23)
where

a i+1/2 = W ′ (ε n i+1 ) -W ′ (ε n i ) ρ 0 ε n i+1 -ε n i . ( 24
)
To avoid divisions by zero when

ε n i = ε n i+1 , one computes a i+1/2 = W ′′ (ε n i ) ρ 0 , (25)
in this particular case. The eigenvalues of the matrix ( 23) are {-s i+1/2 , s i+1/2 }, where s i+1/2 = √ a i+1/2 . Since the strain energy function W is convex over the hyperbolicity domain, the coefficient a i+1/2 is positive. Therefore, the eigenvalue s i+1/2 is real, and the property 1. is satisfied. The matrix (47) is a Roe matrix.

For later use, we introduce the decomposition of q n i+1q n i in the basis of right eigenvectors of the Roe matrix

p 1 i+1/2 = 1, s i+1/2 ⊤ , p 2 i+1/2 = 1, -s i+1/2 ⊤ . ( 26
)
The jump of the solution can be expanded as

q n i+1 -q n i = 2 k=1 α k i+1/2 p k i+1/2 , = 2 k=1 W k i+1/2 , (27)
with the coefficients

α 1 i+1/2 = 1 2 ε n i+1 -ε n i + v n i+1 -v n i s i+1/2 , α 2 i+1/2 = 1 2 ε n i+1 -ε n i - v n i+1 -v n i s i+1/2 . ( 28
)

Flux limiter

We describe now the flux-limiter scheme. The numerical flux in (18) takes the form (section 15.4 in [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF])

f n i+1/2 = f L i+1/2 + f H i+1/2 , (29)
where f L i+1/2 is the flux of Godunov's method or one of its approximated versions, and f H i+1/2 is a higher-order correction deduced from the Lax-Wendroff method. Here, the Godunov flux is approximated linearly by the classical Roe flux, and the Roe matrix (47) is used. Thus,

f L i+1/2 = 1 2 f(q n i ) + f(q n i+1 ) - 1 2 s i+1/2 W 1 i+1/2 + W 2 i+1/2 , (30)
where W k i+1/2 is defined in [START_REF] Mccall | Theoretical study of nonlinear elastic wave propagation[END_REF]. The higher-order correction of the fluxlimiter method writes

f H i+1/2 = 1 2 s i+1/2 1 - ∆t ∆x s i+1/2 × φ(θ 1 i+1/2 )W 1 i+1/2 + φ(θ 2 i+1/2 )W 2 i+1/2 , (31)
where

θ 1 i+1/2 = W 1 i+3/2 • W 1 i+1/2 W 1 i+1/2 • W 1 i+1/2 , θ 2 i+1/2 = W 2 i-1/2 • W 2 i+1/2 W 2 i+1/2 • W 2 i+1/2 , ( 32
)
and φ is a limiter function. Here, the minmod limiter

φ(θ) = max{0, min{1, θ}} (33)
is used. The weights φ(θ k i+1/2 ) are designed to avoid spurious oscillations in the numerical solution. If the weights φ(θ k i+1/2 ) equal one in [START_REF] Berjamin | RiemannElasto1D Toolbox[END_REF], then the Lax-Wendroff method based on the Roe matrix (47) is recovered. In the case of nonconvex flux functions, the minmod limiter (33) has shown better convergence properties than other limiter functions. Similar observations are reported in [START_REF] Kurganov | Adaptive semidiscrete centralupwind schemes for nonconvex hyperbolic conservation laws[END_REF][START_REF] Voss | Exact Riemann solution for the Euler equations with nonconvex and nonsmooth equation of state[END_REF].

To carry out one iteration in time at some grid node i, the numerical values of q at the grid nodes i -2, . . . , i + 2 are required. Therefore, two "ghost cells" must be added on the left and on the right of the numerical domain. Here, a zero-order extrapolation of the numerical solution can be used to achieve outflow boundary conditions. Thus, we simply set

q n -2 = q n -1 = q n 0 , q n N +2 = q n N +1 = q n N , (34)
at each time step [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF].

Numerical test case

The orders of magnitude of the elastic parameters for Berea sandstone ρ 0 , M 0 , β, δ in table 1 have been taken in table 1 and figure 5-(a) of [START_REF] Renaud | Revealing highly complex elastic nonlinear (anelastic) behavior of Earth materials applying a new probe: Dynamic acoustoelastic testing[END_REF]. From [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF], it follows that the model ( 15) is hyperbolic if |ε| < 5.7 × 10 -5 . The numerical domain is [x 0 , x N ] = [-0.5, 0.5] m. To avoid instability, the Courant number ( 19) is set to κ = 0.95.

This test is carried out to validate the ability of the numerical scheme to represent accurately nonsmooth solutions. Here, no volume force is applied in the material. We consider piecewise constant initial data with a single discontinuity at the abscissa x = 0, i.e. a Riemann problem. Here, the ρ 0 (kg/m 3 ) M 0 (GPa) β δ γ (J/m 3 ) τ (ms) 2.2 × 10 3 14 50 10 8 1.0 10 initial data is ε = 2 × 10 -5 , v = 0 if x < 0, and ε = -2.6 × 10 -5 , v = 0.097 m/s if x > 0. The analytical solution to this initial-value problem is displayed in figure 1-(a). Details about its computation can be found in [START_REF] Berjamin | Analytical solution to 1D nonlinear elastodynamics with general constitutive laws[END_REF]. Moreover, an interactive application and a Matlab toolbox [START_REF] Berjamin | RiemannElasto1D Toolbox[END_REF] can be found at http://gchiavassa.perso.centrale-marseille.fr/RiemannElasto/.

The numerical solution (44) is computed up to t = 0.19 ms on a grid with N points. Figures 1-(b) and(c) compare the analytical solution with the numerical solution. On figure 1-(b), the left-going discontinuity is represented (shock wave). On figure 1-(c), the right-going compound wave with a continuous part followed by a discontinuous part is represented (rarefaction shock wave). Both waves are well-located and accurately computed by the numerical scheme. In particular, no spurious oscillations arise, and no convergence failure is noticed as N increases from 1000 to 2000.

Coupling with slow dynamics

As specified in the introduction, nonlinear elastodynamics presented in section 2 is not sufficient to describe the softening of the material. In [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internal-variable model[END_REF], a simple approach was proposed to account for such phenomena. The model was presented in 3D, and analytical computations were carried out in the longitudinal case. Here, we show how the previous numerical method is adapted to this model of slow dynamics.

Modified equations of motion

Internal-variable model.. An internal-variable g in [0, 1[ is introduced, to represent the softening of the material. The internal energy per unit of reference volume ρ 0 e = W (ε) becomes [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internal-variable model[END_REF]

ρ 0 e = (1 -g) W (ε) + Φ(g) , ( 35
)
where e is the specific internal energy. In (35), Φ(g) represents a storage energy with expression where γ > 0 is an energy per unit volume, but other choices are suitable. As discussed in [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internal-variable model[END_REF], the storage energy ( 36) is chosen such that equilibrium points of the model are unique, g = 0 is an equilibrium point, and g is bounded by 1. With such a choice, one observes that the internal energy per unit of reference volume ( 35) is equal to the strain energy W when g is equal to zero.

Φ(g) = - 1 2 γ ln(1 -g 2 ) , (36)
The second principle of thermodynamics yields the axial component of the Cauchy stress [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internal-variable model[END_REF] σ

= (1 -g) W ′ (ε) (37)
instead of W ′ (ε) [START_REF] Favrie | Fast and slow dynamics in a nonlinear elastic bar excited by longitudinal vibrations[END_REF]. The simplest thermodynamically admissible choice of evolution equation is [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internal-variable model[END_REF]

τ 1 ġ = W (ε) -Φ ′ (g) , (38)
where τ 1 = τ × 1 J/m 3 and τ > 0 is a time constant. With the laws (37) and ( 38), the dissipation per unit of reference volume is [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internal-variable model[END_REF]

ρ 0 T η = -ρ 0 ∂e ∂g ġ , = τ 1 ( ġ) 2 0 , (39)
where T > 0 is the absolute temperature, and η is the specific entropy.

Conservation laws with relaxation.. The variable g is added to the vector of unknowns q, so that q = (ε, v, g) ⊤ . The system of conservation laws (15) becomes

∂ t q + ∂ x f(q) = r(q) + s , (40)
where

f(q) = -(v, (1 -g) W ′ (ε)/ρ 0 , 0) ⊤ , r(q) = (0, 0, (W (ε) -Φ ′ (g)) /τ 1 ) ⊤ , s = (0, f v /ρ 0 , 0) ⊤ . (41)
The Jacobian matrix f ′ (q) of f has the eigenvalues {-c(q), c(q), 0}, where

c(q) = (1 -g) W ′′ (ε) ρ 0 (42)
is the speed of sound. The speed of sound ( 42) is real and nonzero -in other words, the system (40) is strictly hyperbolic -provided that g < 1 and W ′′ (ε) > 0. This second condition was already required in the elastic case [START_REF] Maugin | The saga of internal variables of state in continuum thermo-mechanics (1893-2013)[END_REF]. Now let us examine the spectrum of the relaxation function in (40). The Jacobian matrix r ′ (q) of r has the eigenvalues {0, 0, -Φ ′′ (g)/τ 1 }. The expression (36) ensures that r ′ (q) is negative semi-definite. Its spectral radius is

Υ = γ τ 1 1 + g 2 (1 -g 2 ) 2 , (43)
which involves the relaxation time τ 1 /γ characteristic of the slow dynamics [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internal-variable model[END_REF].

4.2. Modified numerical method Conservative scheme.. Applying the same explicit discretization as ( 18) to (40) yields

q n+1 i = q n i - ∆t ∆x (f n i+1/2 -f n i-1/2 ) + ∆t r(q n i ) + ∆t ∆x s n i , (44)
where the numerical flux f n i+1/2 is specified later on. According to the timemarching formula (44), explicit single-step time integration is used.

Numerical stability imposes a bound of the form

∆t min ∆x c n max , 2 Υ n max , (45)
where c n max is the maximum sound velocity (42) that is encountered at time t n , and Υ n max is the maximum spectral radius of the relaxation function (43). Since the relaxation time of the slow dynamics is much larger than the period of exciting signals, the condition (45) reduces to the classical CFL condition [START_REF] Murnaghan | Finite deformations of an elastic solid[END_REF]. If the constitutive law ε → W ′ (ε) is convex or concave, then the maximum sound speed at time t n is given by [START_REF] Ostrovsky | Dynamic nonlinear elasticity in geomaterials[END_REF], with the sound speed (42). If the constitutive law is neither convex nor concave, then (21) becomes

c n max = max 0 i<N 1 -min{g n i , g n i+1 } ρ 0 max ε∈D i W ′′ (ε) , (46)
where D i is the interval with bounds ε n i and ε n i+1 .

Roe linearization.. In the case of the system (40), the following matrix is obtained from [START_REF] Berjamin | Analytical solution to 1D nonlinear elastodynamics with general constitutive laws[END_REF]:

A i+1/2 =    0 -1 0 -a i+1/2 0 b i+1/2 0 0 0    , (47)
where

a i+1/2 = σ n i+1 -σ n i ρ 0 ε n i+1 -ε n i + g n i+1 -g n i ε n i+1 -ε n i b i+1/2 , b i+1/2 = W (ε n i+1 ) -W (ε n i ) ρ 0 ε n i+1 -ε n i , (48)
and the stresses (37) are

σ n i = (1 -g n i ) W ′ (ε n i ).
To avoid divisions by zero when

ε n i = ε n i+1 , one computes a i+1/2 = 1 - g n i + g n i+1 2 W ′′ (ε n i ) ρ 0 , b i+1/2 = W ′ (ε n i ) ρ 0 , (49)
in this particular case. The eigenvalues of the matrix (47) with the coefficients (48)-( 49) are {-s i+1/2 , s i+1/2 , 0}, where the expression s i+1/2 = √ a i+1/2 is unchanged. We rewrite the coefficient a i+1/2 from (48) as

a i+1/2 = W ′ (ε n i+1 ) -W ′ (ε n i ) ρ 0 ε n i+1 -ε n i -g n i W (ε n i+1 ) -W (ε n i ) -W ′ (ε n i ) ε n i+1 -ε n i ρ 0 ε n i+1 -ε n i 2 -g n i+1 W (ε n i ) -W (ε n i+1 ) -W ′ (ε n i+1 ) ε n i -ε n i+1 ρ 0 ε n i+1 -ε n i 2 . ( 50
)
One can note that the first term in (50) corresponds to the elastic case [START_REF] Holzapfel | Nonlinear Solid Mechanics: A Continuum Approach for Engineering[END_REF]. Moreover, the next terms vanish when g n i = 0 = g n i+1 , i.e. when no softening occurs. Since the strain energy function W is convex and g < 1 over the hyperbolicity domain, the coefficient a i+1/2 in (50) is positive. Finally, the eigenvalue s i+1/2 is real, and the property 1. is satisfied. The matrix (47) is a Roe matrix for (40). Now, the jump q n i+1q n i is decomposed in the basis of right eigenvectors of the Roe matrix (47)

p 1 i+1/2 = 1, s i+1/2 , 0 ⊤ , p 2 i+1/2 = 1, -s i+1/2 , 0 ⊤ , p 3 i+1/2 = 1, 0, a i+1/2 b i+1/2 ⊤ . ( 51
)
The jump of the solution can be expanded as a sum ( 27) of three

W k i+1/2 = α k i+1/2 p k i+1/2
, where the coefficients are

α 1 i+1/2 = 1 2 v n i+1 -v n i s i+1/2 + σ n i+1 -σ n i ρ 0 s i+1/2 2 , α 2 i+1/2 = 1 2 - v n i+1 -v n i s i+1/2 + σ n i+1 -σ n i ρ 0 s i+1/2 2 , α 3 i+1/2 = ε n i+1 -ε n i - σ n i+1 -σ n i ρ 0 s i+1/2 2 . (52)
Since the eigenvalue corresponding to k = 3 in the decomposition of the jump is zero, W 3 i+1/2 does not appear in the numerical flux [START_REF] Kurganov | Adaptive semidiscrete centralupwind schemes for nonconvex hyperbolic conservation laws[END_REF] of the fluxlimiter method. Therefore, the formulas ( 29) to [START_REF] Shu | High order weighted essentially nonoscillatory schemes for convection dominated problems[END_REF] can be applied without modification with the flux function (41), the eigenvalue s i+1/2 deduced from (48)-(49), and the decomposition of the jump (51)-(52).

Pulse propagation

Figure 3 of [START_REF] Renaud | Revealing highly complex elastic nonlinear (anelastic) behavior of Earth materials applying a new probe: Dynamic acoustoelastic testing[END_REF] provides an order of magnitude of the characteristic time of the slow dynamics, and an order of magnitude of the average softening ∆M/M for a given strain amplitude V . According to [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internal-variable model[END_REF], the characteristic time of the slow dynamics is τ 1 /γ, and the average softening satisfies ∆M/M ≃ -1 4 (M 0 /γ + 6δ) V 2 . If this scaling rule is combined with the values in [START_REF] Renaud | Revealing highly complex elastic nonlinear (anelastic) behavior of Earth materials applying a new probe: Dynamic acoustoelastic testing[END_REF], then a negative value of γ is obtained, which is not reliable. Thus, the parameters γ and τ in table 1 are not chosen to reach quantitative agreement with [START_REF] Renaud | Revealing highly complex elastic nonlinear (anelastic) behavior of Earth materials applying a new probe: Dynamic acoustoelastic testing[END_REF], but rather to obtain comparable orders of magnitude of the softening's characteristics.

The numerical domain and the Courant number are the same as in section 3.4, but here, the softening is taken into account. The initial data is zero, and a point load f v = δ s (x)ϕ(t) located at the abscissa x s = 0 generates a smooth pulse with angular frequency ω c = 2π f c , with f c = 6 kHz. The source signal in (40) satisfies

ϕ(t) = 2M 0 V 4 m=1 a m sin(2 m-1 ω c t) , (53)
where V = 6.64×10 -6 , a 1 = 1, a 2 = -21/32, a 3 = 63/768, and a 4 = -1/512. It is turned on from t = 0 to t = 1/f c , which corresponds to one fundamental period.

The reference solution is an oversampled numerical solution (44) computed at t = 0.19 ms on a grid with N ref = 2 15 points (figure 2-(a)). The pulse injected at x s = 0 has propagated both towards increasing x and towards decreasing x. By symmetry, only the right-going part is displayed here. No shock wave is observed: the waveform is slightly distorted but still smooth. Figure 2-(b) illustrates the convergence of the numerical method. For a given coarse spatial discretization where N ≪ N ref , the numerical solution is computed up to t ≈ 0.19 ms, and is compared to the reference numerical solution at the same final time. The evolution of the L 2 global error between both strain waveforms is represented in figure 2-(b) with respect to ∆x = 1/N. One can observe that the order of convergence is between one and two.

Dynamic acousto-elasticity

The setup is the same as in the previous section, but here, the point load generates a sinusoidal strain with amplitude V ≈ 10 -6 and angular frequency ω c = 2π f c , with f c = 5 kHz:

ϕ(t) = 2M 0 V sin(ω c t) . (54)
The source (54) is turned on from t = 0 to t = 40 ms, which corresponds to 200 periods of signal. A receiver records the numerical solution at the abscissa x r = 0.

M = ρ 0 c 2 -M 0 M 0 (55)
deduced from the sound speed (42), which are recorded at the position x r of the receiver. A slow decrease of the elastic modulus combined with fast oscillations is observed until the source is stopped. Simultaneously, the strain signals are smooth sinusoids, so that N = 80 is sufficient. After the source is sopped, the elastic modulus recovers slowly its initial value, while the strain is equal to zero. Figure 3-(b) focuses on the steady-state solution. Here, the last 80 numerical values before t = 40 ms are used, which corresponds to two periods of signal at the frequency f c . When ∆M/M from (55) is represented with respect to the strain recorded at the position of the receiver x r , a hysteresis curve is obtained. The orders of magnitude of these phenomena -duration of the transients, magnitude of the softening, size and shape of the hysteresis curves -are very similar to those reported in [START_REF] Renaud | Revealing highly complex elastic nonlinear (anelastic) behavior of Earth materials applying a new probe: Dynamic acoustoelastic testing[END_REF] for Berea sandstone. In particular, one can note that the average softening ∆M/M is proportional to V 2 , as predicted in [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internal-variable model[END_REF]. However, only qualitative agreement with figure 5-(a) of [START_REF] Renaud | Revealing highly complex elastic nonlinear (anelastic) behavior of Earth materials applying a new probe: Dynamic acoustoelastic testing[END_REF] is obtained.

Conclusion

A finite-volume method with flux-limiters for nonlinear longitudinal elastodynamics is implemented, with various polynomial constitutive laws. Based on Roe linearization, this numerical method is well-suited to the present system of partial differential equations, and has been validated with the Riemann problem of nonlinear elastodynamics. Nonsmooth solutions such as shock waves are well-captured, even in the case of nonconvex stress-strain relationships. The method is adapted to a case with an additional scalar evolution equation, which is deduced from the internal-variable model [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internal-variable model[END_REF] of slow dynamics. The numerical model reproduces qualitatively experimental observations related to dynamic acoustoelasticity.

We mention here a few improvements to be introduced. Boundary conditions such as free edges and oscillating walls need to be implemented so as to reproduce real configurations in a more realistic way (section 7 of [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF]). Due to resonance, viscoelastic attenuation needs then to be incorporated in the model, which has been carried out similarly to [START_REF] Favrie | Fast and slow dynamics in a nonlinear elastic bar excited by longitudinal vibrations[END_REF] in recent works [START_REF] Berjamin | Modeling longitudinal wave propagation in nonlinear viscoelastic solids with softening[END_REF]. Currently, higher-order shock-capturing methods [START_REF] Shu | High order weighted essentially nonoscillatory schemes for convection dominated problems[END_REF] are developed, as well as similar methods in multiple space dimensions.

Figure 1 :

 1 Figure 1: (a) Analytical solution to a Riemann problem. Comparison with the numerical solution (b) on the left-going strain wave, (c) and on the right-going strain wave.

  3 m. The numerical solution (44) is computed up to t = 80 ms on a grid with N = 80 points, which corresponds to 40 points per wavelength at the angular frequency ω c . The computations are performed in C++. Each simulation

Figure 2 :

 2 Figure 2: Assessment of the numerical method's convergence on the propagation of a smooth pulse. Reference solution at t = 0.19 ms (a): strain, particle velocity, and softening variable. (b) Error measurement.

20

 20

V = 1 . 5 × 10 - 6 Figure 3 :

 151063 Figure 3: (a) Strain signal and softening of the material, as recorded by the receiver at the abscissa x r = 0.3 m. (b) Hysteresis curves in steady-state. The arrow indicates how time increases along the curve.

Figure . 4 :

 4 Figure .4: Sketch of the hyperbolicity domains (white) of Murnaghan's law (12) with respect to the parameter ϑ = (l + 2m)/(λ + 2µ).

Table 1 :

 1 Physical parameters of Berea sandstone.

Acknowledgments

This work was supported by the interdisciplinary mission of CNRS (IN-FINITI). The project leading to this publication has received funding from Excellence Initiative of Aix-Marseille University -A*MIDEX, a French "Investissements d'Avenir" programme. It has been carried out in the framework of the Labex MEC. Opération réalisée avec le concours du Programme d'investissements d'avenir du Gouvernement français dont la gestion a été confiée l'Andra.

Appendix

In the case of Murnaghan's law [START_REF] Lyakhovsky | Non-linear elastic behaviour of damaged rocks[END_REF], strict hyperbolicity of ( 15) is ensured if

Therefore, one must have

where

The constitutive law ( 12) is convex if σ ′′ (ε) 0, i.e.

Hence, Murnaghan's law is locally concave at small strains (ε ≃ 0) provided that ϑ -3/2, otherwise it is locally convex at small strains. The inflexion point

is represented in figure .4, as well as the hyperbolicity domains (.2) of Murnaghan's law.


Kenneth Sörensen
email: kenneth.sorensen@ua.ac.be

Marc Sevaux
email: marc.sevaux@univ-ubs.fr

A practical approach for robust and flexible vehicle routing using metaheuristics and Monte Carlo sampling

Keywords: stochastic vehicle routing, robustness, flexibility, Monte Carlo sampling, metaheuristics, memetic algorithm

A practical approach for robust and flexible vehicle routing using metaheuristics and Monte Carlo sampling

Robust and flexible vehicle routing

The objective of vehicle routing problems is to determine the order in which to visit a spatially distributed set of customers using a set of vehicles so that the total travel cost (distance, time) is minimized. Standard vehicle routing formulations however, assume that all data concerning customer demand, travel costs, etc. are known with perfect certainty at design time. For many reasons, e.g. the uncertainty related to traffic conditions, these assumptions are unwarranted in a large number of practical situations. As a result, a number of vehicle routing formulations designed to find robust solutions have been developed. Three types of such formulations, commonly referred to as stochastic vehicle routing problems can be distinguished: problems with stochastic travel times (e.g. [START_REF] Lambert | Designing collection routes through bank branches[END_REF]), problems with stochastic demands (e.g. [START_REF] Stewart | Stochastic vehicle routing: A comprehensive approach[END_REF]), and problems with stochastic customers, i.e. in which each customer has a certain probability of requiring service (e.g. [START_REF] Bertsimas | A vehicle routing problem with stochastic demand[END_REF]). For a survey, we refer to [START_REF] Gendreau | Stochastic vehicle routing[END_REF].

For stochastic problems, robust or flexible solutions are required. Consistent with the definitions found in the literature (see e.g. [START_REF] Kouvelis | Robust Discrete Optimisation and its Applications, volume 14 of Nonconvex Optimization and its Applications[END_REF], [START_REF] Mulvey | Robust optimization of large-scale systems[END_REF], [START_REF] Roy | A missing link in OR-DA: Robustness analysis[END_REF], [START_REF] Vincke | Robust solutions and methods in decision aid[END_REF], [START_REF] Dias | On computing ELECTRE's credibility indices under partial information[END_REF], [START_REF] Branke | Evolutionary Optimization in Dynamic Environments[END_REF]), we call a solution to a stochastic optimization problem robust if it has a high quality regardless of the actual realization of the stochastic parameters. In other words, a solution that consistently has a high performance across all possible outcomes of the stochastic parameters, is called robust. Usually, stochastic optimization algorithms attempt to find the solution that has either the best expected quality or the best worst-case quality over all potential outcomes of the stochastic parameters. In many cases however, a solution can be (partially) changed once the actual values of the stochastic parameters become known. In this case, a solution is preferred that can be successfully adapted to any realization of the stochastic parameters. Such a solution is called flexible. In other words, a flexible solution is one that has a high quality after adaptation to the outcomes of the stochastic parameters, whereas a robust solution is one that has a high quality without adaptation to the stochastic parameters. In the remainder, we will call the procedure that is used to adapt the solution to the actual values of the stochastic parameters a repair procedure.

Several problems arise with these methods, that inhibit their use in real-life vehicle routing applications. First, most of the methods are exact methods, implying that they are unable to solve problems of realistic size. Indeed, as indicated by [START_REF] Birge | Stochastic programming computation and applications[END_REF], the complexity of stochastic programs grows proportionally to the number of possible realisations of the stochastic parameters, which in turn grows exponentially with the number of stochastic parameters. Because in realistic cases this number is usually very large or even infinite (in the case of continuous distributions), only very small problems can be solved to "optimality". To partially overcome this problem, a part of the stochastic programming literature has fused on methods that use some form of Monte Carlo sampling of the stochastic parameters. For stochastic linear programming problems, stochastic decomposition [START_REF] Higle | Stochastic decomposition: An algorithm for two stage stochastic linear programs with recourse[END_REF] and importance sampling [START_REF] Infanger | Planning Under Uncertainty: Solving Large Scale Stochastic Linear Programs[END_REF] have been developed. For discrete stochastic problems, Norkin et al. [1998a,b] develop a sampling branch-and-bound. Whereas these methods use sampling at different steps of the optimization method to estimate e.g. function values, the sample average approximation method [START_REF] Kleywegt | The sample average approximation method for stochastic discrete optimization[END_REF] uses sampling to generate a set of sample scenarios and then attempts to optimize the corresponding deterministic expected-value problem. This method has been successfully applied to supply chain design problems [START_REF] Goetschalckx | Designing flexible and robust supply chains[END_REF] and routing problems [START_REF] Verweij | The sample average approximation method applied to stochastic routing problems: A computational study[END_REF], among others.

A second drawback to traditional stochastic programming methods is the fact that they make heavy use of the specific stochastic structure of the problem they consider and are very difficult to adapt to other problems. Thirdly, almost all algorithms developed for stochastic routing attempt to find the solution with the optimal expected value. Although the average performance of a solution is certainly an important indication of its robustness, it does not take into account the risk-preference of the decision maker.

In this paper, we develop a method that is able to overcome these problems. We adopt a pragmatic approach to stochastic optimization and approximate the stochastically optimal solution on two levels. First, we use a sampling-based approach to quickly estimate the robustness of a solution. This allows us to quickly and simultaneously evaluate both average and worst-case performance of a solution, but more complicated measures of robustness can be incorporated. Our sampling-based approach also allows us to evaluate solutions of problems with many stochastic parameters. Secondly, we use a metaheuristic to find the solution with the best approximate robustness characteristic, which allows us to solve large problems. Our approach adapts a metaheuristic for a deterministic problem by replacing its evaluation function by a so-called robust evaluation function.

Problem description and deterministic solution method

The CVRP is defined on an undirected graph G = (V, E) with a set of nodes V = {0, 1, . . . , n}. Node 0 corresponds to the depot, that has a set of identical vehicles of capacity Q and maximal travel cost C. Nodes 1 to n represent a set of spatially distributed customers, the demand of which is given by q i . The travel cost between customer i and customer j is given by c ij , the weight of the edge between node i and node j. The objective of the deterministic VRP is to find a set of minimum total cost routes that have the following properties: (1) each route begins and ends at the depot, (2) each customer is visited exactly once, (3) the capacity and maximal travel cost of the vehicles is not exceeded.

To solve this problem, we have developed a MA|PM, or menetic algorithm with population management. A MA|PM is a memetic algorithm (a GA hybridised with local search) that uses distances to measure and control the diversity of a small population. This allows to maintain diversity in the population while keeping the size of the population small. For a more elaborate description of MA|PM, we refer to [START_REF] Sörensen | MA|PM: memetic algorithms with population management[END_REF].

In our MA|PM, solutions are represented as a string of customers, without trip delimiters S = S 1 S 2 . . . S n . When necessary (e.g. to calculate their objective function value) they are decoded optimally using the split decoding procedure by [START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF]. This pro-cedure finds the optimal points at which to split the solution into tours so that the resulting decoded solution is feasible and the total travel cost is minimized. It makes use of an auxiliary directed weighted graph H = (X, F, W ). The vertex set X contains n + 1 nodes (where n is the number of customers served), labelled from 0 (the depot) to n.

The edge set F is constructed as follows. An edge (i, j) from vertex i to vertex j (i < j) is added when a trip containing customers S i+1 to S j is feasible, i.e. if

j k=i+1 q S k ≤ Q (1)
and

w ij = c 0S i+1 + j-1 k=i+1 c S k S k+1 + c S j 0 + j k=i+1 e S k ≤ C. (2)
The weight w ij of the arc (i, j) is equal to the sum of the travel distances covered in the trip 0 → S i+1 → S i+2 → . . . → S j → 0, plus the sum of drop costs in this trip. The minimum-cost path from vertex 0 to vertex n in H gives the shortest set of trips corresponding to this encoding. If more than one min-cost path exists, each of these paths will give a set of trips with equal total length. Since the graph H contains only positive weights and no circuits, the shortest path can be efficiently calculated using Bellman's algorithm [START_REF] Bellman | On a routing problem[END_REF].

A second procedure used heavily by our MA|PM for the CVRP is a distance measure based on the edit distance. Given three possible edit operations (add character, remove character and substitute character), the edit distance d(s, t) is defined as the minimal number of edit operations required to transform a string s into another string t. The edit distance is heavily studied in the literature. A simple dynamic programming algorithm [START_REF] Wagner | The string-to-string correction problem[END_REF] calculates this measure in O(n 2 ) where n is the length of both strings. Other, more efficient algorithms have been developed [START_REF] Sörensen | A framework for robust and flexible optimisation using metaheuristics with applications in supply chain design[END_REF].

Vehicle routing solutions can be regarded as sets of strings, each string representing a tour. As tours have no direction, each string may be reversed. The distance measure for the CVRP developed in [START_REF] Sörensen | A framework for robust and flexible optimisation using metaheuristics with applications in supply chain design[END_REF] calculates the minimal number of edit operations required to transform a VRP solution into another one, keeping into account that both solutions consist of a set of reversable strings. This is done by assigning tours from the first solution to the second one so that the total number of required edit operations is minimized.

Given two solutions s and t having m and n trips respectively, the distance measure constructs a square matrix M of size max(m, n). Each column and each row represents a trip of s and t respectively. Empty trips are added to the solution with the smallest number of trips to give both solutions an equal number of trips. Cell (i, j) of M contains the minimum of d(s i , t j ) and d(s i , tj ), where s i is the i-th trip of solution s and tj is the reversed j-th trip of solution t. A linear assignment problem is then calculated to find the matching of trips in solution s to trips in solution t. The cost of the assignment problem is the distance between the two solutions.

In the MA|PM, the distance measure is used for population management in order to maintain the diversity of a small population of high-quality individuals and overcome problems of premature convergence. Before a (candidate) solution is added to the population, the algorithm checks to see whether it satisfies the population diversity criterion, i.e. whether it is sufficiently different from all other members of the population. This is done by measuring the so-called distance to the population, i.e. the minimum distance of the candidate solution to any solution already in the population. The solution cannot be added until its distance to the population is greater than or equal to the value of the population diversity parameter ∆. This value can be used to control the diversity of the population as an increase of ∆ will tend to increase population diversity, while a small value will tend to decrease it. Intensification and diversification phases can be alternated by setting the value of ∆ to appropriate levels [START_REF] Sörensen | MA|PM: memetic algorithms with population management[END_REF].

The MA|PM also uses a fast local optimization algorithm, that consists of two simple tabu search procedures that are used alternatively until neither of them finds any more improvements. The insert tabu search procedure attempts to locate any customer at any other potential location in the solution (i.e. any other location in any other tour) and takes the move that decreases the total distance most (or increases it least when no improving moves can be found). An insertion of a customer into another location is only allowed when no constraints are violated by this move and when the customer does not appear on the tabu list. This insert procedure is repeated until the best-found solution during this tabu search phase has not been improved for a fixed number of moves. A fixed-length tabu list is maintained to avoid cycling. The tabu list contains the customers that were most recently moved. The swap tabu search procedure attempts to swap every pair of customers and swaps the pair that yields the largest improvement in objective function value. The swap of a pair of customers is only allowed when the resulting solution does not violate any constraints and when the pair of customers does not appear on the tabu list. This procedure is repeated until a fixed number of non-improving moves. In this case, the tabu list contains pairs of of items that were most recently swapped. Pairs that appear on the tabu list are excluded from swapping.

The crossover operator used in our MA|PM is the linear ordered crossover (LOX). This operator starts by randomly dividing both parents into three parts. The first offspring solution is found by copying the middle part from parent one and completing the solution by circularly scanning parent two, starting from the third part and wrapping to the first part. The second offspring is formed by reversing the roles of the parents.

Our MA|PM for the CVRP starts by generating a set of random permutations of all customers. These are then decoded (using the split procedure) and subjected to local search (alternating the two tabu search operators until no more improvement is found). The initial population consists of a small number (usually 10) solutions. At each generation, two solutions are drawn from the population using a binary tournament selection operator (i.e. we retain the best of two randomly chosen solutions). These two solutions are subjected to the LOX operator. Both offspring solutions are then decoded and subjected to the tabu search local searches. When the solution satisfies the diversity criterion, it is added to the population. If not, it is mutated until it does satisfy the diversity criterion and added to the population. The mutation operator randomly swaps two customers and repeats this step until the resulting solution satisfies the diversity criterion.

We use an adaptive population management scheme, in which the value of ∆ is slowly decreased to allow it to intensify and increased again when the search appears to be stuck in a local optimum, i.e. no more improving solutions are found for a fixed number of generations.

Some experiments show that our MA|PM is competitive with the best-known approaches in the literature. For detailed results, we refer to [START_REF] Sörensen | A framework for robust and flexible optimisation using metaheuristics with applications in supply chain design[END_REF].

A sampling-based approach for robust and flexible vehicle routing

In this section, we modify the MA|PM developed in the previous section, so that it is able to deal with stochastic vehicle routing problems. This is done by replacing the objective function by a so-called robust evaluation function, that measures the robustness or flexibility of the solution. A repair function is used if the solution can be adapted when the actual values of the stochastic parameters become known.

In a nutshell, our approach is the following. We use a metaheuristic, such as our MA|PM, that generates a set of solutions that are both diverse and have a high quality. For each solution we encounter, we calculate one or more measures of robustness or flexibility, referred to as a robust evaluation function value. Finally, we choose the solution that has the best robust evaluation function value. When two or more measures of robustness or flexibility are calculated, we choose the solution that has the best combined value using a multi-objective decision making process.

A robust evaluation function value for a given solution is calculated by repeatedly applying the solution to a sampling of the stochastic parameters and calculating the corresponding (deterministic) objective function value. These objective function values are then combined into one or more measures of robustness. If a repair procedure is available, the solution is repaired before it is evaluated. The concept of repair function is closely related to the stochastic programming concept of recourse, but differs from it that any procedure, even a heuristic one, may be used.

A typical robust evaluation function f * is of the form

f * (x) = 1 n e ne i=1 f (x, S i (π)),
where f is the (deterministic) objective function, π is the set of stochastic parameters, and S is the sampling function. S i (π) represents the i-th sampling of the stochastic parameters and f (x, S i (π)) is the objective function value of solution x when applied to the i-th sampling. n e is the number of deterministic evalutions per robust evaluations, i.e. the number of evaluations of the solution on a random sampling of the stochastic parameters.

If the solution can be repaired, the repair function R is invoked before evaluation the solution and the robust evaluation function becomes

f * (x) = 1 n e ne i=1 f (R(x, S i (π))).
Other useful robust evaluation functions include the worst-case performance, given (for a minimization problem) by

f * wc (x) = max i f (x, S i (π)),
or the standard deviation

σ * (x) = 1 n e -1 ne i=1 [f (x; S i (π)) -f * (x)] 2 .
These measures explicitly incorporate the risk preference of the decision maker into the process. E.g. the solution that has the best worst-case performance will generally be more "conservative" or "risk-averse" than the one that has the best average-case performance. The decision maker may decide to find the solution that minimizes f * (x)+λσ * (x), where λ is a parameter expressing the risk-averseness of the decision maker. Other, even more complex measures of robustness may be used, such as the probability that the quality of the solution falls below a certain threshold. Of course, since they are obtained using Monte Carlo simulation, the robust evaluation function values are only estimates of the true robustness measures. However, statistical theory can be used to estimate the reliability of the estimates. For a sufficiently large number of evaluations (e.g. n e > 30), the central limit theorem states that f * (x) is approximately normally distributed. From this, it follows that an 100(1α) confidence interval for the real average performance of a solution is given by

f * (x) ± z 1-α 2 (σ * (x)) 2 n e .
For a more elaborate discussion, we refer to [START_REF] Law | Simulation Modeling and Analysis[END_REF].

One possible extension of our method is the use of penalty functions to penalize violation of constraints. For some samples of the stochastic parameters, a solution might become infeasible. A penalty function P may be used to express the severity of the constraint violations and added to the robust evaluation function value as follows:

f * (x) = 1 n e ne i=1 [f (x; S i (π)) + P(x; S i (π))].
In our method, stochastic information can be expressed both as independent probability distributions on the parameters of the problem, or as a set of scenarios in which each of the parameters has a fixed value. In the latter case, the value of n e is the number of scenarios, and the sampling function S is used to evaluate the performance of a solution under the different scenarios.

We believe that our approach has several advantages over traditional methods based on stochastic programming, most importantly the fact that it is considerably more flexible, considerably easier to implement and applicable to much larger and much more complex problems. Moreover, it allows the decision maker to evaluate complex robustness characteristics of a solution, using any type of stochastic information that is available and incorporating his risk preference.

Experiments

In this section, we perform some experiments on different stochastic versions of the CVRP. Unless otherwise indicated, our MA|PM uses a population of size 10. The population diversity parameter ∆ is set to 20 initially and reduced by one unit for each 3 generations. ∆ is set to its initial level after 60 unproductive generations. Both the swap tabu search and the move tabu search procedures use a tabu tenure of 30 and repeat until 10 -improving moves.

Data sets for the stochastic versions are derived from the [START_REF] Christofides | Combinatorial optimization[END_REF] instances for the deterministic CVRP, available from the OR Library1 . We generally assume that the distributions of the stochastic parameters are independent and have their mean equal to the corresponding value in the deterministic data set. Unless otherwise indicated, the robust evaluation function value is based on n e = 1000 evaluations.

Vehicle routing with stochastic demand and cost

In this set of experiments, customer demand q i is assumed to be uniformly distributed between 0.75q i and 1.25q i . Travel cost between customers i and j are uniformly distributed between 0.8c ij and 1.2c ij . qi and cij are the averages of the demand and travel cost distributions respectively and equal to the values found in the deterministic data set. Note that the deterministic problem is the mean value formulation of the stochastic problem, i.e. the deterministic problem can be obtained by setting all stochastic parameters to their expected value.

We assume that the routes have to be determined before the actual values of customers demand and travel cost are known and cannot be changed afterwards. We therefore attempt to find robust solutions. For some samples of the stochastic demand and cost, a solution might violate the maximum travel cost constraints per route or the capacity constraints of the vehicles. We therefore introduce two penalty functions. Let c ijk and q ik represent the k-th random numbers taken from the distributions of c ij and q i respectively. The deterministic objective function value of solution x, evaluated on the k-th sampling of the stochastic parameters can be calculated as follows.

f (x; S k (π)) = l e ij ∈E l (x) c ijk ,
where E l (x) represents the set of edges included in the l-th tour of solution x, e ij the edge between customers i and j. If the total cost in a given route is larger than C (the maximum cost), then a fixed penalty per unit of exceeded cost is added.

P cost (x; S k (π)) = α 1 l max(0, e ij ∈E l (x) c ijk -C).
Similarly, if the total demand served in a given route exceeds the maximum demand Q, a fixed penalty cost per unit of exceeded demand is added.

P capacity (x; S k (π)) = α 2 l max(0, v i ∈V l (x) q ik -Q), (3)
where V l (x) represents the set of vertices (customers) in tour l of solution x. α 1 is the penalty cost per unit of exceeded travel cost. α 2 is the penalty per unit of exceeded capacity. We set α 1 = 100 and α 2 = 500.

We look for solutions that have both a good average-case and a good worst-case performance. To this end, we define two robust evaluation functions: f * (x) and f * wc . To calculate f * for solution x, we calculate the total travel cost increased with the values of both penalty functions for n e random samples of the stochastic parameters and take the average. For f * wc , we take the maximum value. For each solution encountered, we also measure the standard deviation σ * (x) and the value f (x), the deterministic objective function value.

The MA|PM is applied to the 14 vehicle routing problems. Robust evaluation function f * is used in the binary tournament selection procedure to select solutions from the population for crossover. This guides the search towards more robust solutions. We stop the MA|PM after 200 generations, yielding 400 solutions. Detailed results are in appendix (table3). Table1 is an extract from this table holding the results of a single experiment. For each data set, three rows are presented. The first row contains the characteristics of the solution that minimizes f , i.e. the solution that would be found by our method in the deterministic case, if we would ignore all stochastic information. The second and third rows contain the characteristics of the solution that minimizes f * and f * wc respectively. It can be seen that the solutions in the second and third rows score considerably better with respect to all robustness indicators (f * , f * wc and σ * ). These two solutions are therefore considerably more robust than the solution in the first row, proving that their is a strong need to take the stochastic information into account. This is further illustrated in figure1, where we plot the values of f versus f * for all 400 solutions. From this figure, it is clear that all robust solutions (low f * ) have a good deterministic objective function value, but that the reverse is not necessarily true. The decision maker can choose on of the two robust solutions presented in table3 or alternatively one of the other non-dominated solutions (solutions for which there does not exist a solution which scores better with respect to all robust evaluation functions) generated by the MA|PM. These are not shown because of space restrictions.

Data file n

Vehicle routing with stochastic customers and a repair function

In some situations, the vehicle dispatcher does not know in advance which customers will require service and which not. A possible strategy is to route all customers and repair the solution once the list of customers that do require service becomes known. This situation may occur when customers need to book a service in advance, but may cancel at short notice. A solution that can be effectively repaired is called flexible.

In an experiment, each customer has a 50% probability of requiring service. A 100% service level is required and the company therefore decides to design conservative solutions that are still feasible if all customers in a route require service. The repair procedure used by the company is to remove all customers that do not require service from a given route. If the customer in i-th position in a given route is removed, the vehicle travels directly from the customer in position i -1 to the customer in position i + 1. Total travel cost is measured after the customers have been removed from the routes.

The MA|PM is run for 200 generations. Each of the 400 generated solutions is evaluated 1000 times on a random list of customers. As in the previous experiment, we record the average f * , the maximum f * wc and the standard deviation σ * over all n e evaluations. Results can be found in table4.

From this table, it is clear that not much profit can be gained by using the robust optimization approach, as the solution that has the best value of f has approximately the same robustness characteristics as the solution that optimizes one of the robust evaluation function values. The main reason for this is the requirement that the solution should be feasible even when all customers are present. This results in very conservative solutions that-when the actual customers to serve are observed-use only a fraction of the available resources (capacity and maximum cost). The best solution of the deterministic problem, utilising as much of these resources as possible, is therefore likely to be a better solution for the reduced set of customers, than other solutions that do not use the resources as fully. As a result, this solution is very likely to be flexible.

Use of the robust evaluation function and computation times

The use of our approach for robust and flexible vehicle routing requires some design issues to be tackled. Especially the question when to use the robust evaluation function, should be taken with some care. Obviously, the robust evaluation function is used to evaluate the robustness properties of a solution. However, it can also be used to guide the search towards more robust solutions. In our MA|PM, the robust evaluation function is used to guide the search in the binary tournament selection process. However, it is not used for the move selection of the tabu search local optimization procedures. In both tabu search procedures, the (deterministic) objective function value is used to select the next move. The reason for this is that the use of the robust evaluation function for move selection in the tabu search procedures would have resulted in prohibitive computing time increases. The bulk of computing time is spent on improving solutions using the tabu search procedures. While this is necessary to obtain high quality solutions, the tabu search procedures require a large number of evaluations. Moreover, the calculation of the cost of a move can easily be short-circuited when a deterministic objective function is used, but this no longer holds when a robust evaluation function is used. Indeed, it can be easily seen that that when a customer k is inserted between customers i and j, the total cost of this route changes by c ik +c kj -c ij (provided that the route remains feasible). Similar reasoning applies for removing a customer from a route or swapping two customers. These shortcut calculations allow a move to be evaluated without re-evaluating the entire solution, but they are no longer valid for a robust evaluation function.

By using the robust evaluation function for selection of solutions from the population, it only has to be calculated twice for each generation (when both solutions are added to the population). This results in relatively small computing time increases as can be seen in table 2. This table compares the computing times for 200 generations for the MA|PM in the deterministic and the stochastic case (experiment with stochastic demand and cost). Computation time never increases by more than a factor of 3.08.

One of the design choices that influence the computation time when using our robust optimization approach is the value of n e , the number of objective function evaluations to perform to calculate the value of the robust evaluation function. In the experiments performed in this paper, we put the value of n e at 1000 and still obtained relatively small computing time increases, but this may not always be the case. Especially in situations where the objective function is time-intensive to calculate, the value of n e must be carefully chosen. As we mentioned before, a confidence interval can be calculated around the value of f * . The value of n e should be chosen in such a way that the confidence interval is small enough. How small exactly the confidence interval should be, depends on the specific application, the level of confidence that is required by the decision maker and the computing time available.

Conclusions

In this paper, we have developed an approach to find robust and flexible solutions of several stochastic versions of the capacitated vehicle routing problem. Our approach combines the optimization power of metaheuristics with Monte Carlo simulation based estimation of the robustness or flexibility of a solution. A repair procedure is introduced when flexible solutions are needed and penalty functions may be introduced to penalize violation of constraints. Our approach recognizes the need for more complex expressions of robustness or flexibility than the often-used average performance to be entered into the decision-making process.

We have tested our approach on vehicle routing problems with stochastic demand and cost and with stochastic customers. We have also discussed some design issues, specifically how the search can be guided towards robust or flexible solutions and how the computing time increase can be kept within limits.

We believe that our approach offers considerable advantages over more traditional methods based on stochastic programming, especially when applied to large-scale, reallife stochastic vehicle routing applications.

Figure 1 :

 1 Figure 1: Ordinary evaluation and robust evaluation for vrpnc01

Table 1 :

 1 Vehicle routing with stochastic demand and costs, example result

		Criterion	f	f *	σ *	f * wc
	vrpnc01 50	f	549.76 3211.60 2891.00 17779.56
		f * 1	604.76 605.01	13.21	876.66
		f * 2	629.12 629.28	10.61	661.48

Table 2 :

 2 Computing times (200 generations)

	Data file	n	avg. determ. (s) avg. stoch. (s) avg. stoch. avg. determ.
	vrpnc01	50	68.84	99.30	1.44
	vrpnc02	75	205.59	408.59	1.99
	vrpnc03 100	363.23	573.17	1.58
	vrpnc04 150	932.70	1856.12	1.99
	vrpnc05 199	1969.29	4266.21	2.17
	vrpnc06	50	92.57	97.52	1.05
	vrpnc07	75	205.43	423.80	2.06
	vrpnc08 100	369.90	655.88	1.77
	vrpnc09 150	967.54	2329.98	2.41
	vrpnc10 199	1657.46	5102.75	3.08
	vrpnc11 120	667.53	1215.21	1.82
	vrpnc12 100	292.09	698.23	2.39
	vrpnc13 120	772.81	1008.80	1.31
	vrpnc14 100	526.89	722.35	1.37
	All	N/A	649.42	1389.85	2.14

Table 3 :

 3 Vehicle routing with stochastic demand and costs

	Data file	n	Criterion	f	f *	σ *	f * wc	CPU (s)
	vrpnc01	50	f	549.76	3211.60 2891.00 17779.56	101
			f * 1 f * 2	604.76 629.12	605.01 629.28	13.21 10.61	876.66 661.48
	vrpnc02	75	f	891.86	8335.36 4986.75 30836.52	349
			f * 1 f * 2	956.57 956.57	1134.44 1134.44	617.56 617.56	6509.91 6509.91
	vrpnc03	100	f	887.97	2904.39 2753.68 16150.60	541
			f * 1 f * 2	971.38 971.38	977.51 977.51	72.99 72.99	2236.30 2236.30
	vrpnc04	150	f	1149.53	5313.22 3869.58 23450.44	1893
			f * 1 f * 2	1208.36 1208.36	1596.03 1063.23 1596.03 1063.23	9598.45 9598.45
	vrpnc05	199	f	1512.77 10112.13 5737.19 33182.36	4153
			f * 1 f * 2	1687.18 1687.18	3040.09 2118.64 15630.12 3040.09 2118.64 15630.12
	vrpnc06	50	f	559.87	567.68	122.64	3350.03	98
			f * 1 f * 2	559.87 579.51	567.68 579.42	122.64 9.58	3350.03 605.79
	vrpnc07	75	f	995.55	2583.40 2285.52 11686.45	415
			f * 1 f * 2	1070.72 1070.72	1077.22 1077.22	79.39 79.39	2512.38 2512.38
	vrpnc08	100	f	945.33	946.03	38.49	2106.94	648
			f * 1 f * 2	945.33 965.35	946.03 965.21	38.49 11.81	2106.94 999.43
	vrpnc09	150	f	1352.35	1506.98	659.95	6979.06	2265
			f * 1 f * 2	1404.10 1404.10	1404.68 1404.68	16.24 16.24	1460.92 1460.92
	vrpnc10	199	f	1597.51	3683.01 2822.85 18401.49	5117
			f * 1 f * 2	1697.23 1784.14	1698.18 1783.08	31.61 18.10	2415.77 1837.73
	vrpnc11	120	f	1260.55	4775.36 3331.22 21620.35	1229
			f * 1 f * 2	1492.17 1507.27	1504.12 1507.52	212.15 26.71	6115.15 1601.84
	vrpnc12	100	f	883.55	8004.79 5372.87 30825.31	698
			f * 1 f * 2	1138.66 1138.66	1638.14 1359.21 1638.14 1359.21	9534.97 9534.97
	vrpnc13	120	f	1368.33	4629.15 3227.07 18598.48	1079
			f * 1 f * 2	1565.28 1727.20	1611.20 1734.02	315.01 73.03	6819.88 3117.71
	vrpnc14	100	f	952.46	7706.49 5126.85 29321.32	719
			f * 1 f *	1033.58	1775.53 1707.77 13073.67

Table 4 :

 4 Vehicle routing with stochastic customers and repair procedure

	Data file	n	Criterion	f	f *	σ *	f * wc	CPU (s)
	vrpnc01	50	f	527.46	424.39 28.30	494.89 99.30
			f * 1 f * 2	528.22 528.22	423.92 28.85 423.92 28.85	488.88 488.88
	vrpnc02	75	f	905.22	736.65 45.95	854.83 408.59
			f * 1 f * 2	942.20 927.41	721.21 43.89 728.53 47.00	846.89 836.21
	vrpnc03	100	f	852.05	703.72 34.09	788.00 573.17
			f * 1 f * 2	852.14 852.14	696.88 32.00 696.88 32.00	776.41 776.41
	vrpnc04	150	f	1130.42	977.25 32.34 1067.67 1856.12
			f * 1 f * 2	1130.42 1134.48	977.25 32.34 1067.67 979.88 33.22 1065.78
	vrpnc05	199	f	1466.17 1286.28 32.03 1366.59 4266.21
			f * 1 f * 2	1483.99 1282.70 39.68 1384.81 1466.17 1286.28 32.03 1366.59
	vrpnc06	50	f	567.92	457.43 32.55	533.51 97.52
			f * 1 f * 2	577.01 569.51	456.09 34.14 459.30 32.16	553.61 528.00
	vrpnc07	75	f	990.71	809.90 51.78	919.41 423.80
			f * 1 f * 2	1000.02 990.71	806.21 56.00 809.90 51.78	945.15 919.41
	vrpnc08	100	f	933.72	786.23 36.69	872.24 655.88
			f * 1 f * 2	942.00 933.72	779.18 40.24 786.23 36.69	890.61 872.24
	vrpnc09	150	f	1351.54 1152.77 40.75 1257.38 2329.98
			f * 1 f * 2	1351.54 1152.77 40.75 1257.38 1351.54 1152.77 40.75 1257.38
	vrpnc10	199	f	1574.00 1379.91 41.69 1497.29 5102.75
			f * 1 f * 2	1577.60 1377.88 42.73 1483.68 1577.60 1377.88 42.73 1483.68
	vrpnc11	120	f	1122.20 1019.96 21.55 1083.10 1215.21
			f * 1 f * 2	1123.32 1019.26 22.61 1081.70 1125.26 1020.72 21.19 1074.73
	vrpnc12	100	f	861.27	774.77 24.27	839.91 698.23
			f * 1 f * 2	862.98 862.98	771.46 24.77 771.46 24.77	820.59 820.59
	vrpnc13	120	f	1239.06 1121.11 22.49 1177.47 1008.80
			f * 1 f * 2	1250.98 1117.50 23.57 1177.95 1239.06 1121.11 22.49 1177.47
	vrpnc14	100	f	878.46	768.03 23.92	833.47 722.35
			f * 1 f * 2	878.46 878.46	768.03 23.92 768.03 23.92	833.47 833.47

http://people.brunel.ac.uk/ ∼ mastjjb/jeb/orlib/vrpinfo.html

1125.71 2192.62 1891.78 12757.21


Md Umar Hashmi
email: md-umar.hashmi@inria.fr

Ana Bušić

Limiting Energy Storage Cycles of Operation

 to efficiently control the number of cycles of operation of a battery. The cycles of operation have to be tuned based on price volatility to maximize the battery life and arbitrage gains. We propose a mechanism to distinguish arbitrage returns.

I. INTRODUCTION

With increased share of power coming from Renewable Energy Sources (RES), the uncertainty in power networks is increasing. This variability makes the power, frequency and voltage regulation more challenging. Inverter based RES has small inertia compared to conventional synchronous power generators, reducing the capability of the future power grid to withstand perturbation. The future power network will need fast dynamic ancillary services to avoid additional infrastructure development in installing fast ramping power generators. Distributed energy storage devices such as batteries are suitable for assisting the grid because of its fast ramping capability. Authors in [START_REF] Divya | Battery energy storage technology for power systemsan overview[END_REF] highlight the diverse roles energy storage technologies can play in future power networks. Energy storage not only improves the reliability of power network but also facilitates arbitrage and increases the value of renewable energy sources in the energy markets. Authors in [START_REF] Xu | Factoring the cycle aging cost of batteries participating in electricity markets[END_REF] point that battery owners participating in electricity markets should consider cycle aging of batteries. This is due to finite battery cycle life, which is affected by charge-discharge operational cycles, temperature etc. Using a case study on ISO New England energy and reserve markets, authors demonstrate that participants maximize their operational profit while considering cycle aging cost.

The variability in the price is due to the variable generation cost and demand and supply mismatch components caused due to scheduling delays and errors in forecast of RES. This variability in the price can be used by batteries for performing arbitrage. However, sometimes charge and discharge cycles generate very small revenue and it would be beneficial for the battery to be idle. Maximizing idle time would lower stress on the battery and maximize the cycle life of the battery [START_REF] Bušić | Distributed control of a fleet of batteries[END_REF].

The key contributions are:

• Limiting cycles of operation: In this work we highlight how at the cost of small portion of arbitrage gains, the battery cycle life can be significantly improved.

• By introducing friction coefficient, η friction in energy arbitrage control we can efficiently control the cycles of operation. If the battery is operating more number of cycles per day then η friction is decreased so that the battery on average operates for lower number of cycles. Note the selection of η friction will governed by variability of price, number of cycles of operation in absence of any friction, i.e., η friction = 1.

• Distinguishing arbitrage returns: We present a mechanism to eliminate low returning transactions to maximize battery life.

The rest of the paper is organized as follows. In Section II we describe the cycle and calendar life of battery and list the typical values for power grid batteries. Section III summarizes the optimal arbitrage algorithm proposed in [START_REF] Hashmi | Optimal control of storage under time varying electricity prices[END_REF]. Section IV presents the mechanism to limit cycles of operation of a battery. Section V presents the numerical results. Section VI concludes the work.

II. BATTERY LIFE

There are two distinct, but inter-related ways that are used to measure the life of Li-ion batteries deployed in energy storage system applications for smart grid -calendar life and cycle life.

A. Cycle Life

One of the indexes which decides the life of operation is the rated cycles of operation of a battery. This information is often provided by manufacturers in their data-sheets. In this section we will enlist the values of the cycle of operation provided by the manufacturers under controlled environment. This value of cycles will with high probability lead to the end-of-life (EoL) of the battery. EoL is frequently defined as a state of the battery when the maximum capacity of the battery reduces to 80% of its rated initial capacity. In order to understand the typical cycle life of Li-ion based grid batteries, we summarize the cycle life information provided by the manufacturers in their data-sheets in Table I. [START_REF] Cdtechno | Featherweight lithium-ion battery li tel 48-170 c[END_REF] 4,300 cycles 20 o C, 80%DoD AES Advancion [START_REF] Aes | Advancion energy storage for the grid[END_REF] 2000 to 3000 at 80%DOD LG Chem 6.4 [START_REF] Lg-Chem | [END_REF] > 6000 90% DoD, 25 o C QINOUS Li-ion [START_REF] Qinous | Qinous ess compact, medium, large[END_REF] 4000 cycles 1C/1C, 80%DoD, 25 o C Tesla Powerwall 1 [START_REF]Tesla-powerwall[END_REF] 5000 cycles daily cycle application Tesvolt 8000 cycles for 70% DoD Li-phosphate [START_REF] Tesvolt | Tps energy storage[END_REF] 5000 cycles for 90% DoD 6 000, no aging 25 deg.C, 70% DOD Forsee Power 5 000, with aging 25 o C , 2 cycles/day HE48 [START_REF] Forseepower | Lithium-ion he 48[END_REF] 4400, with aging 25 o C , 1 cycles/day aentron Li-ion [START_REF] Aentron | Off grid battery system[END_REF] 3000 50% DoD, 20 o C Bosch Storage [START_REF] Bosch | Bosch energy storage solutions[END_REF] 4500 -REFU [START_REF] Refu-Elektronik-Gmbh | 003-001-0300[END_REF] 6000 80% DoD, 20 o C , max discharge=0.5C

The factors affecting the rated cycles of operation are:

• Battery type: Chemical composition of the battery,

• Depth of discharge (DoD): Operational cycles,

• Temperature should be maintained,

• Ramp rate: affects the stress on the battery,

• Cycles per day. Fig. 1 shows the effect of number of cycles of operation with variation in DoD. It should be noted that 10 cycles of 10% DoD will be equal to 1 cycle of 100% DoD, in terms of energy exchange.

B. Calendar Life

Calendar life refers to the number of years the battery is expected to last till the battery will reach EoL. It is independent of how much the battery is charged and discharged. However, calendar life is dependent on the state of charge of the battery and the temperature.

TABLE II: Calendar Life as listed in datasheets

Manufacturer

Calendar life and Remarks AES Advancion [START_REF] Aes | Advancion energy storage for the grid[END_REF] 7-8 years Lifespan Evolion Li-ion battery [START_REF] Saft-Batteries | Evolion, proven ultra-compact solution for on and off-grid installations[END_REF] at 20 deg. C is 20 years at 40 deg. C is > 10 years LG Chems RESU 6.4 EX [START_REF] Lg-Chem | [END_REF] at 25 deg. C is > 10 years Tesla Powerwall 2 [START_REF]Tesla-powerwall[END_REF] warranty of 10 years Tesvolt Li-phosphate [START_REF] Tesvolt | Tps energy storage[END_REF] 10 years Forsee Power Li-ion HE48 [START_REF] Forseepower | Lithium-ion he 48[END_REF] 10 years aentron Off grid Li-ion [START_REF] Aentron | Off grid battery system[END_REF] 10 yr. at 20 deg. C, Warranty 5 yr. Bosch Storage Solutions [START_REF] Bosch | Bosch energy storage solutions[END_REF] > 10 years BYD [START_REF] Byd | Byd home based energy storage solution[END_REF] Warranty of 5 years

Table II indicate that the typical value of calendar life is around 10 years. The calendar life deteriorates drastically with temperature as shown in Table II for Evolion Li-ion battery.

The life of a battery will be a function of cycle life and calendar life. If the battery operates more cycles per unit of time then the battery EoL will be achieved due to cycle life limitations and vice versa. Data-sheet [START_REF] Saft | Li-ion battery life[END_REF] shows how the calendar life of the battery is a function of SoC level and temperature. An interesting observation is that a high SoC level affects the calendar life drastically at high temperatures. It can be observed that below 35 o C expected life is independent of SoC. In this work we assume temperature is controlled and SoC is regulated in the optimal band of operation.

III. OPTIMAL STORAGE CONTROL

In a commodity market, arbitrage can be performed if the price of that commodity varies in time. Often, the buy and sell decisions are coupled in time because of the finite size of the buffer/capacity/inventory, where the commodity is stored. Similarly, when energy storage devices perform arbitrage, the buy and sell decisions are coupled because of the finite storage size; a decision to buy or sell in the present time instant will influence the potential to charge or discharge in future.

The total duration, T , of operation is divided into N equal steps, where in each step i ∈ {1, ..., N } the price of electricity is denoted as p i . The sampling time is denoted as h. The efficiency of charging and discharging of the battery are denoted by η ch ∈ (0, 1] and η dis ∈ (0, 1], respectively. We denote the change in the energy level of the battery at i th instant by x i = hδ i , where δ i denotes the storage ramp rate at i th instant; δ i > 0 implies charging and δ i < 0 implies discharging. Amount of energy that comes from the battery at i th instant is given by s

i = 1 ηch [x i ] + -η dis [x i ] -. Storage Constraints:
The battery has limited capacity. We incorporate the battery capacity constraints by imposing b i ∈ [b min , b max ], ∀i, where b min , b max is the minimum and the maximum battery charge level. In operating Li-ion battery it should be ensured that it is never over-charged or overdischarged. The instantaneous battery charge level is given as

b i = b i-1 + x i .
The ramp rate constraint is given as δ i ∈ [δ min , δ max ], ∀i, where δ min ≤ 0, δ max ≥ 0 are the minimum and the maximum ramp rates (kW). The change in battery charge level x i ∈ [X min , X max ], where X min = δ min h and X max = δ max h.

An optimal control problem for an energy storage device such as a battery is proposed in [START_REF] Cruise | Optimal control of storage for arbitrage, with applications to energy systems[END_REF], [START_REF] Cruise | Optimal control of storage incorporating market impact and with energy applications[END_REF] using convexity property of the optimization function and saddle point inequality. Authors in [START_REF] Hashmi | Optimal control of storage under time varying electricity prices[END_REF] propose optimal energy arbitrage algorithm under time varying electricity prices. In this work it is shown that a time horizon of optimization can be subdivided into sub-horizons. In each of these sub-horizons, the shadow price is a function of price levels in a sub-horizon and is independent of all past and future sub-horizons. This value of the shadow price is altered only when the battery capacity reaches its maximum or minimum permissible charge levels. Based on the value of shadow price in a sub-horizon, the optimal control decisions are selected depending on the level of price of electricity at that instant.

The threshold based structure of the optimal solution is presented in Remark 1.

Remark 1. The optimal control decision x * i in the ith instant minimizes the function

C (i) storage (x)-µ * i x for x ∈ [X min , X max ].
The optimal decision x * i (µ) is

x * i (µ) =                X min , if µ < p dis (i), [X min , 0] , if µ = p dis (i), 0, if p ch (i) > µ > p dis (i), [0, X max ] , if µ = p ch (i), X max , if µ > p ch (i), (1)
where

p ch (i) = p i /η ch , p dis (i) = p i η dis , C (i)
storage (x i ) = s i p i and µ represents the shadow price of the transaction.

Note for µ = p ch (i) or µ = p dis (i), x * i (µ) takes an envelope of values and for any other value of µ it is a singleton set.

In order to find optimal decisions among an envelope of possible solutions based on the price variations, the Backward

Step algorithm is used one time. The details of the algorithm can be found in prior work [START_REF] Hashmi | Optimal control of storage under time varying electricity prices[END_REF].

IV. LIMITING CYCLES OF OPERATION

In Section II we discussed that the cycles of operation is one of the primary indicators of battery life. In other words, if the battery is used a lot then the battery will last less. To maximize the operational life, the average number of cycles of operation per unit of time (I op ) should be set to

I op = Cycle Life Calendar Life . (2)
EoL of a battery operating more than I op will be caused by the limitation on cycle life of the battery, implying the battery is over-used. If the battery is under-used, EoL will be caused by the calendar life. We assume that if cycle life equals calendar life of the battery than the operational life the battery is maximized. This criterion will be met if battery operates I op cycles per unit time on average, shown in Eq. 2.

Note the volatility across all days in a year can vary drastically, implying battery should perform more cycles when the volatility is high. The friction coefficient ideally should be selected adaptively ensuring maximizing gains and battery's operational life.

We propose limiting the operation of the battery adding friction in mode changes for the battery. Adding friction will ensure that the battery does not operate for lower returns. This idea of creating dead-band of no operation is motivated by Eq. 1. The threshold based structure of the optimal solution indicates that the optimal decision for the battery when the condition p ch (i) > µ > p dis (i) is true, is to do nothing. This band signifies the additional profit the charge discharge cycle of a battery should make in order to compensate the losses incurred due losses in charging and discharging. Increase of this band will indicate eliminating low returns transactions in arbitrage. We define modified charging and discharging cost as a function of friction coefficient denoted as

p L ch (i) = p ch (i)/η friction = p i η ch η friction p L dis (i) = p dis (i)η friction = p i η dis η friction
The threshold based structure of the solution is modified as

x * i (µ) =                X min , if µ < p L dis (i), [X min , 0] , if µ = p L dis (i), 0, if p L ch (i) > µ > p L dis (i), [0, X max ] , if µ = p L ch (i), X max , if µ > p L ch (i), (3)
The selection of η friction will be governed by: volatility of the price, maximizing operational life of the battery and arbitrage gains. Note the value of η friction ∈ [0, 1]. The control of η friction should only be considered if the cycles of operation per day is higher than I op .

V. NUMERICAL RESULTS

In this section we present a numerical example demonstrating the effect of friction coefficient on arbitrage gains, number of cycles of operation and the shadow price of transaction. The battery parameters used for the numerical evaluation is listed in Table III. Consider the battery has a cycle life of 6000 cycles and arbitrage gains with friction coefficient. When friction coefficient is 1 implies battery operates more and the arbitrage gains are higher. As the coefficient is reduced the decrease in gains and cycles of operation can be observed. For I op cycles per day, the optimal value of friction coefficient for this numerical example is η frictionopt = 0.855 as shown in Fig. 3. Note that for achieving I op cycles the reduction in arbitrage gains is ≈ 10% and the cycles of operation decreases by ≈ 36%. A risk averse user may prefer to maximize it's current gains, for example for η friction = 0.928, the reduction in gains is 1.47% and the reduction in cycles of operation will be 11.7%. Fig. 4 shows the variation of shadow price

Fig. 1 :

 1 Fig.1: Li-ion cycle of operation with DoD[START_REF] Forseepower | Lithium-ion he 48[END_REF]

Fig. 2 :

 2 Fig.2: Price of Electricity NYISO[START_REF] Nyiso | Real time locational marginal price[END_REF]

Fig. 5 :

 5 Fig. 5: Variation of ramp rate of the battery with η friction

TABLE I :

 I Cycles of operation as listed in datasheets

	Manufacturer	Cycles of operation	Remarks
	C&D Techno

TABLE III :

 III Battery Parameters

		Parameter	Value	Parameter	Value
		δmax		0.5 kW	δ min		-0.5 kW
		bmax		1 kWh	b min		0.1 kWh
		η ch		0.95	η dis		0.95
		h		0.25 hour	b 0		0.1 kWh
	and calendar life of 10 years. The optimal number of cycles
	the battery should operate per day should be I op = 1.644
	per day, using Eq. 2. Real-time 15-minute locational marginal
	price data for 21st December, 2016 from NYISO [20] is
	used for simulation analysis. The electricity price is plotted
	in Fig. 2. Fig. 3 shows the variation of cycles of operation
		0.14
							Price of Electricity
		0.12
		0.1
	$/kWh	0.08
		0.06
		0.04
		0	4	8	12	16	20	24
					Hour of Day

Research supported by the French National Research Agency grant ANR-16-CE05-0008, and by PGMO (Gaspard Monge Program for Optimization and operations research).

M.U.H. and A.

of ramp rate of the battery with η friction . On comparing Fig. 5 with respect to price variations in Fig. 2, it is evident that battery operates only when the gains are significantly higher for lower values of η friction . As the η friction is reduced the battery operates only around midnight to charge in order to discharge during morning peak and during evening peak. The friction coefficient distinguishes arbitrage gains based on the returns. Inclusion of η friction ensures battery operates for only when the gains are high enough.

VI. CONCLUSION

We propose a mechanism to compare arbitrage gains and cycles of operation for a battery. The users can tune the coefficient of friction according the price volatility over a period to maximize operational life of a battery and its arbitrage gains. The inclusion of friction coefficient ensures that low returning transactions are eliminated and battery operates only when returns are higher.