K V Avramov

USING NONLINEAR NORMAL MODES TO ANALYZE FORCED VIBRATIONS

Keywords: nonlinear forced vibrations, discrete systems, Rauscher's method, Pierre-Shaw nonlinear modes, shallow arch

The paper proposes a method to analyze forced vibrations in nonlinear systems. The procedure combines Rauscher's method and Pierre-Shaw nonlinear modes. Results from an analysis of the forced vibrations of a shallow arch are presented as an example

&&

( , , ) cos( ), , q q R q q h t i m

i i i i m i + + = = w 2 1 1 K W , (1)
undergong forced vibrations. Assume that the vibrations are not damped. The vector function R q q i m ( ,... , )

1 describes nonlinear terms that contain the generalized coordinates raised to second, third, and other powers. It is assumed that the natural frequencies w i do not satisfy the internal-resonance condition.

Assume that the excitation frequency W is close to the natural frequency w l . Since there is no internal resonance, we can assume that the amplitude of the generalized coordinate q l is much higher than all the other generalized coordinates in the neighborhood of the principal resonance. For the purpose of vibration analysis, we will use an iterative procedure. At the first iteration, we set q l ¹ 0 and q m = 0, m =1, m, m ¹ l. Then an oscillator with one degree of freedom follows from system (1): && ~( ) cos( ) q q R q h t l l l l l l

+ + = w 2 W , (2)
where ~( ) ( ,... , , , ,... , ) R q R q l l l l = 0 0 0 0 . Following Rauscher's method, we represent the solution of system (2) as cos ( ) Wt r q q q l l l

= = + + + a a a 0 1 2 2 K . ( 3
)
The idea of combining nonlinear normal modes and Rauscher's method is as follows. Substituting (3) into system (1), we obtain a so-called pseudoautonomous dynamic system [START_REF] Avramov | Snap-through truss as an absorber of forced oscillations[END_REF][START_REF] Vakakis | Normal Modes and Localization in Nonlinear Systems[END_REF]. Then nonlinear normal modes are analyzed in the pseudoautonomous dynamic system [START_REF] Shaw | Normal modes for nonlinear vibratory systems[END_REF]. Using nonlinear normal modes, we reduce the system of m equations (1) to one equation analyzed approximately.

Let us now obtain the solution of Eq. ( 2) in the form (3). To this end, we apply the harmonic balance method to system (2). The solution is represented as

q A A t A t A t l = + + + + 0 1 2 3 2 3 cos cos cos W W W K . ( 4
)
Restricting the analysis to three Fourier harmonics, we substitute solution (4) into (2) and balance the harmonics. Doing so gives a system of four nonlinear algebraic equations for ( , , , , )

W A A A A 0 1 2 3 . Its general form is F W n n ( , , , , ) , , A A A A 0 1 2 3 0 14 = = . (5)
The objective of the analysis is to plot the amplitude-frequency characteristic at the principal resonance W » w l . This is why the value of A 1 is incremented, and the system of the nonlinear algebraic equations ( 5) is solved for each value of A 1 . As a result, we obtain parameters ( , , , )

W A A A 0 2 3 .
Then, the values of W, , A A 0 2 , and A 3 are used to determine the coefficients of series (3): a 0 , a 1 , a 2 , ….

Let us outline a procedure of determining these coefficients. To this end, we represent the solution of Eq. (4) as

q A A z A A z A z A l = - + - + + + 0 2 1 3 2 2 3 3 3 2 4 ( ) K , z t = cos W . ( 6
)
Next, series (3) is substituted into (6) and the coefficients of like powers of q l are equated. As a result, we obtain a cubic equation for a 0 :

4 2 3 0 3 0 3 2 0 2 1 3 0 0 2 A A A A A A a a a + + - + - = ( ) . (7)
The parameters a 1 , a 2 , and a 3 are defined by a a a

1 1 3 2 0 3 0 2 1 3 4 12 = - + + A A A A , a a a a a 2 1 2 3 0 2 1 3 2 0 3 0 2 12 2 3 4 12 = - + - + + ( ) A A A A A A , a a a a a a a a 3 1 2 2 3 1 2 3 0 2 1 3 2 0 3 0 2 4 6 3 4 12 = - + + - + + ( ) A A A A A A A . ( 8
)
Equation ( 3) is substituted into the nonautonomous dynamic system (2) to obtain a pseudoautonomous dynamic system: && ( ,... , ) ( ), , q q R q q h q q i m i i i i m i l l

+ + = + + + = w a a a 2 1 0 1 2 2 1 K . ( 9
)
For the purpose of further analysis, we change variables:

q q t i i i = + 0, ( ) h , (10)
where q i m i 0 1 , ( , ) = are the coordinates of the fixed point of dynamic system (9) that are determined from a system of m nonlinear algebraic equations:

w a a a i i i m i l l q R q q h q q 2 0 01 0 0 1 0 2 0 2 , , , , , ( ,... , ) ( ),
+ = + + + K i m = 1, . (11)
The dynamic system (9) takes the following form for the variables h h 1 , ,

K m : && ( ) , , h a h h h i i j j i m j m F i m + + , , = = = å 1 1 0 1 K , ( 12
)
where F i ( ,... h 1 , h m ) are nonlinear functions.

The matrix

L = = = || || , , a ij j m i m 1 1 is represented as L = - UAU 1 , ( 13
)
where U is a matrix consisting of the eigenvectors of the matrix L;

A m = diag ( , , ) n n 1 2 2

K

. We introduce new variables:

x h = - U 1 , ( [ ], [ ] ) x x x h h h = , , = , , 1 1 K
K m m T . (14)
Then the dynamic system (12) takes the form

&& ( , , ) , , x n x x x i i i i m L i m + + = = 2 1 0 1 K , ( 15
)
where L i ( ,... x 1 , x m ), i m = 1, , is the nonlinear part of the dynamic system.

To analyze the dynamic system (15), we will use Pierre-Shaw nonlinear normal modes [START_REF] Shaw | Normal modes for nonlinear vibratory systems[END_REF]. Then, all the generalized coordinates x i and velocities v i i = & x are expressed in terms of one coordinate x l and velocity v l :

x x x x x k k l l l l l l v F F v F v = = + + + ( , ) ( ) ( ) ( ) 1 1 2 2 1 3 1 2 K , v v v F F v F v k k l l l l l l = = + + + ( , ) ( ) ( ) ( ) x x x 1 2 2 2 2 3 2 2 K , k m = 1, , k l ¹ , (16)
where

F F 1 1 2 1 ( ) ( ) ,
,K are unknown coefficients to be determined.

Let us differentiate Eq. ( 16) with respect to time and use Eq. [START_REF] Nayfeh | Nonlinear Oscillations[END_REF]. Doing so gives the following system of partial differential equations:

v v v v L v k l l k l l k l l l l l l l ( , ) { [ ( , ) x x x x n x x x x = ¶ ¶ - ¶ ¶ + , , , 2 1 K
K, x x m l l v ( , )]}, [ ] n x x x x x x x k k l l k l l l m l l v L v v 2 1 ( , ) ( , ) ( , ) + , , , , K K = - ¶ ¶ + ¶ ¶ + , , , , v v v v L v k l l k l l l l l l l m l x n x x x x x x { [ ( , ) ( , 2 1 K K v l )]}, k m = 1, , k l ¹ . ( 17
)
Following [START_REF] Bespalova | Vibrations of polygonal plates with various boundary conditions[END_REF], we substitute Eqs. ( 16) into ( 17) and collect the coefficients of like terms x l j l j v 1 2 , j 1 0 1 2 = , , ,... , j 2 0 1 2 = , , ,K . As a result, we obtain a system of linear algebraic equations for F F

1 1 2 1 ( ) ( )
, ,K . Solving this system yields functions [START_REF] Rauscher | Steady oscillations of system with nonlinear and unsymmetrical elasticity[END_REF]. Using Eq. ( 14), we write functions ( 16) for the variables h h 1 , ,

K m as h h h h h h k k l l l l l l l s C C s C C s C = = + + + + ( , ) ( ) ( ) ( ) ( ) 1 1 2 1 3 1 2 4 1 5
( ) 1 2 s l +K , s s s C C s C C s C k k l l l l l l l = = + + + + ( , ) ( ) ( ) ( ) ( ) h h h h 1 2 2 2 3 2 2 4 2 5
( ) 2 2 s l +K , ( 18
)
where

s k k = & h , C C 1 1 2 1 ( ) ( )
, ,K are known parameters.

We substitute Eq. ( 18) into [START_REF] Avramov | Snap-through truss as an absorber of forced oscillations[END_REF] and the resulting functions into the equation l of system [START_REF] Gulyaev | Applied Problems in the Theory of Nonlinear Vibrations of Mechanical Systems[END_REF]. Doing so gives a dynamic system with one degree of freedom, which can be represented in the form (2). To solve this equation, we use the harmonic-balance method. The solution is presented in the form (4). Then, we obtain a system of four nonlinear algebraic equations for W, , A A 0 2 , and A 3 . In the general case, this system can be represented in the form [START_REF] Yu | Dynamic analysis of a two-mass system with essentially nonlinear vibration damping[END_REF]. Let us now generalize the iterative procedure outlined above. At the first iteration, one equation ( 2) is solved and the contribution of the coordinates q q q q q l l m 1 2

-+ is neglected. At the second iteration, the nonlinear vibration mode (18) is used to take into account the contribution of all the coordinates to the dynamics of the system. The third, fourth, etc. iterations are set up similarly. If the values of W, , A A 0 2 , and A 3 obtained in two successive iterations are close, the iterative process is terminated. Then new value of A 1 is specified and a new point of the amplitude-frequency characteristic starts to be calculated.

2. Nonlinear Vibrations of a Shallow Arch. The above method can be used to analyze the nonlinear vibrations of a shallow arch (Fig. 1). The nonlinear dynamics of the arch is described by a system of integro-partial differential equations [START_REF] Hsu | Stability of shallow arches against snap-through under timewise step loads[END_REF]:

A y EI y y Hy x L n F t tt xxxx xx r d + - - = - ae è ç ö ø ÷ ( ) cos( ) 0 W , H EA L y y dx x x L = - ò 2 2 0 2 0 ( ) , , (19)
where A and I are the cross-sectional area and moment of inertia; y x t ( , )is the current deflection of the arch; E and r are Young's modulus and density of the material; y x 0 ( )is the initial deflection of the arch; d( / ) x L n is the delta-function; n is an arbitrary number; and F t cos( ) W is a concentrated force. We introduce the following dimensionless variables and parameters:

x x L * = , y y r * = , t t EI Al * = r 4
, y ry

0 0 = * , d d * * , x n L x L n f FL EIr - ae è ç ö ø ÷ = - ae è ç ö ø ÷ = 1 3 , ( 20
)
where r is the cross-sectional radius of inertia.

Omittng the asterisks in [START_REF] Shaw | Modal analysis-based reduced-order models for nonlinear structures-an invariant manifolds approach[END_REF], we represent the dynamic system [START_REF] Shaw | Normal modes for nonlinear vibratory systems[END_REF] in the form

+ - - - = - ae è ç ö ø ÷ ò ( ) ( ) co , 0 2 0 2 0 1 2 1 d s( ) Wt . (21)
The initial deflection of the arch can be represented as

y x x 0 1 2 2 = + l p l p sin( ) sin( ). (22)
Then the vibrations of the arch can be described by

y x t q t x q t x ( , ) ( ( ))sin( ) ( ( ))sin( ) = + + + l p l p 1 1 2 2 2 . ( 23
)
Next we set n = 2, substitute Eqs. ( 22) and ( 23) into [START_REF] Vakakis | Normal Modes and Localization in Nonlinear Systems[END_REF], and apply the Bubnov-Galerkin method. Doing so gives the following dynamic system:

&& ( ) ( , ) c o s ( ) q q q q q f t 1 4 1 4 1 1 1 2 4 2 + + + = p p l m W , && ( ) ( , ) q q q q q 2 4 2 4 2 2 1 2 16 0 + + + = p p l m , m l l ( , ) ( ) ( ) q q q q q q 1 2 1 1 1 2 2 2 2 4 2 = + + + . ( 24
)
Now we introduce normal coordinates ( , )

x x 1 2 of the linear part of system (24), which are determined from the linear relations q q r r x x ] / .

1 2 1 2 1 2 1 1 é ë ê ù û ú = é ë ê ù û ú é ë ê ù û ú , (25)
The dynamic system (24) takes the following form for the variables ( , )

x x 1 2 :

&& x p x i i i + 2 = + + + + s q i i i i i f t A x A x A x x A x 1 1 1 1 2 2 2 2 3 1 2 4 1 cos( ) ( ) ( ) ( ) ( ) W [ ] 3 5 1 2 2 6 1 2 2 7 2 3 1 2 + + + = A x x A x x A x i i i i ( ) ( ) ( ) , , , (26)
where

q 1 1 = , q 2 1 = -, s l l = - 2 1 2 1 2 2 2 p p
, and f f

1 4 2 = p . The parameters A A 1 1 2 1 ( ) ( )
, ,K are omitted for brevity.

Let us analyze forced vibrations at the principal resonance. In this case, the excitation frequency W is close to p 1 . We set x 2 0 = at the first iteration. Then we obtain the following dynamic system:

[ ] && cos( ) ( ) ( ) x p x f t A x A x 1 1 2 1 1 1 1 1 1 2 4 1 1 3 + = + + s W . ( 27
)
The solution of this equation can be represented in the form (4). Applying the harmonic-balance method, we arrive at a system of nonlinear algebraic equations, which can be written in the form (5). If system (5) has been solved numerically, the solution of Eq. ( 27) could be represented in the form (3). This solution is substituted into (26). As a result, we obtain a pseudoautonomous dynamic system:

&& x p x i i i + 2 ( ) = + + + + + s q i i i i i f r x A x A x A x x A x A 1 1 1 1 2 2 2 2 3 1 2 4 1 3 5 ( ) ( ) ( ) ( ) ( [ ]
i i i x x A x x A x i ) () () , 1 2 2 6 1 2 2 7 2 3 + = 1, 2. (28)
Let us change variables as in [START_REF] Avramov | Snap-through truss as an absorber of forced oscillations[END_REF]. For (28), we have

x x t i i i = + 0, ( ) h , i = 1, 2. (29)
Then the parameters x 0 1 , and x 0 2 , can be found from a system of nonlinear algebraic equations in the form [START_REF] Avramov | Nonlinear equations of flexural-flexural-torsional oscillations of rotating beams with arbitrary cross-section[END_REF]. System (28) can be represented for the variables ( , ) h h

1 2 as [ ] && ( , ) ( , ) ( ) ( )
h h s h h h h h h i i i i i i i p c c P P + = + + + 2 1 1 2 2 2 1 2 3 1 2 , P B B B i i i i ( ) ( ) ( ) ( ) ( , ) 2 1 2 1 1 2 2 2 2 3 2 1 h h h h h h = + + , P B B B B i i i i i ( ) ( ) ( ) ( ) ( ) ( , ) 3 1 2 4 1 3 5 1 2 2 6 2 2 1 7 h h h h h h h h = + + + 2 3 1 2 , , i = , (30)
where the parameters B B

, ,K are omitted for brevity.

We rearrange the linear part of system (30) for the normal coordinates using the equation

h h x x 1 2 1 2 1 2 1 1 é ë ê ù û ú = V V é ë ê ù û ú é ë ê ù û ú , V = -- i i c p c s n s 12 1 2 2 11 , i = 1 2 , . (31)
The frequencies n 1 and n 2 are found as the roots of the biquadratic equation n n

4 2 0 1 0 + + = R R , (32)
where R c c p p System (30) takes the following form for the variables ( , )

0 2 2 1 1 2 2 1 2 = + - - s( ) , R p p c p c p c c c c
x x 1 2 : && [ ( , ) ( , )] ( ) ( )
x n x a x x x x The nonlinear mode (invariant manifold) is represented in the form [START_REF] Rauscher | Steady oscillations of system with nonlinear and unsymmetrical elasticity[END_REF] when k = 2 and l = 1. Using the approach outlined in Sec. 1, we express the coefficients F F F

i i i i i R R + = + 2 2 1 2 3 1 2 , R D D D i i i i ( ) ( ) ( ) ( ) 2 1 1 2 2 2 2 3 1 2 = + + x x xx , R D D D D i i i i i ( ) ( ) ( ) ( ) ( ) 3
1 1 2 1 3 2 ( ) ( ) ( )
, ,... , as

F D F D 1 1 1 2 1 2 2 2 2 2 1 2 2 2 2 2 1 2 2 4 2 ( ) ( ) ( ) ( ) ( ) ( )
, = - - = - a n n n n n a 4 1 2 2 2 n n - , F D F F F 3 1 1 2 2 2 1 2 2 2 2 1 3 2 1 2 2 4 0 ( ) ( ) ( ) ( ) ( ) ( ) , = - = = = a n n n . ( 34
)
The invariant manifolds are represented for the variables ( , , , ) h h

1 2 1 2 s s as h h d h 2 1 1 1 1 1 1 2 3 1 1 2 = V + + + - * ( ) ( ) ( ) F F s K , s s F s 2 1 1 1 2 2 1 1 = V + + - ( ) * d h K , (35)
where d * ( )

( ) = -V V V V -V - 1 2 1 1 2 1 1 2 .
The motion in a nonlinear normal mode is described by a nonlinear ordinary differential equation:

&& [ cos ( ) ( ) ( ) ( ) x p x f t G G x G x G 1 1 2 1 1 1 0 2 1 2 1 2 2 1 2 3 2 + = + + + + s W & &] ( ) ( ) x G x G x x 1 2 4 2 1 3 5 2 1 1 2 + + . ( 36
)
The parameters of system (36) are omitted for brevity. The periodic motions of system (36) are represented in the form (4). Then the parameters ( , , , )

W A A A 0 2
3 follow from the system of nonlinear algebraic equations ~( , , ,

n A A A 0 2 3 0 14 = = ) , , F W n
3. Numerical Vibration Analysis. Let us numerically analyze the forced vibrations of the shallow arch [START_REF] Shaw | Normal modes for nonlinear vibratory systems[END_REF]. We will use the following parameters of the system: l 1 5

= , l 2 =0.8, f =0.05. Let the vibrations occur with a moderate amplitudes about the equilibrium state q q 1 2 0 = = . Snap-through buckling is not considered here. The natural frequencies of linear vibrations Conclusions. The paper proposes an iterative procedure to analyze the forced vibrations of discrete nonlinear systems. This method combines Rauscher's method and nonlinear normal modes. Rauscher's method reduces a nonautonomous dynamic system to an iterative process of finding the normal mode in a pseudoautonomous system.

Nonlinear normal modes are an efficient approach to the analysis of dynamic systems with a great number of degrees of freedom [START_REF] Avramov | Nonlinear equations of flexural-flexural-torsional oscillations of rotating beams with arbitrary cross-section[END_REF][START_REF] Shaw | Modal analysis-based reduced-order models for nonlinear structures-an invariant manifolds approach[END_REF]. Nonlinear normal modes can be applied to a dynamic system of high dimension resulting from finite-element discretization of an elastic medium.

Fig. 1 x

 1 Fig. 1

  omitted for brevity.

about the static equilibrium position are p 1 = 2 .

 12 9067 and p 2 = 5.1158. The procedure presented in Sec. 2 is used to analyze forced vibrations. The amplitude-frequency characteristic of forced vibrations (the modulus of the first Fourier harmonic | | A 1 versus the excitation frequency W 1 ) is shown in Fig. 2. The dashed line represents the skeletal curve, while the solid line the curve of forced vibrations. The right branch of the amplitude-frequency characteristic (relative to the skeletal curve) is determined by the parameter A 1 0 > , and the left branch by A 1 0 < . The surface of the nonlinear normal mode for A 1 = -1.2 is shown in Fig. 3 as an example.

  Fig.2Fig.3


Jokin Munoa
email: jmunoa@ideko.es

Zoltan Dombovari

Alexander Iglesias

Gabor Stepan

Tuneable mass dampers with variable stiffness for chatter suppression

Keywords: stability, chatter, milling, vibration absorber

come

Tuneable mass dampers with variable stiffness for chatter suppression

INTRODUCTION

New self-acting and self-characterizing functions of a future machine tool require integrated equipments to improve precision, productivity and remote maintenance capabilities [Monostori, 2014]. This trend in the machine tool industry intends to decrease human interaction in large manufacturing lines and helps further optimizing machining processes taking into account the overall dynamic behavior of the machine tool/tool/workpiece system.

Embedded solutions are becoming cheaper and increasingly affordable for machine tool constructors, which results in the design of independent active [START_REF] Lu | Magnetic actuator for active damping of boring bars[END_REF], Munoa et al., 2013] and semi-active [START_REF] Aguirre | Self-tuning dynamic vibration absorber for machine tool chatter suppression[END_REF]] mechatronic devices to attenuate machine tool vibrations. Variable stiffness tuneable mass dampers (VSTMD) dampers offers promising advantages since they can be tuned during the machining process even remotely, and therefore they are excellent candidates to be included in embedded semi-active [START_REF] Aguirre | Self-tuning dynamic vibration absorber for machine tool chatter suppression[END_REF] and cyber-physical systems.

The original tuned mass dampers (TMD) provide a classical cheap and effective classical solution to increase damping in a mechanical structure. Exact analytical solution exists to optimize the parameters of TMD when the original structure is considered with a single undamped mode by [Den Hartog, 1934]. The above mentioned classic theory provides a sufficiently good approximation if the cutting performance of the machine tool is limited basically by a single dominant mode with low damping in the structure. If the machine has high damping, the exact tuning of the TMD requires semi-analytical or numerical solutions even for a single mode [Asami, et al. 2002]. Several variations of the primary idea have been proposed in the literature. In recent years, a single TMD with multiple degrees of freedom (MDOF) [Zuo, et al. 2006] and multiple TMDs (MTMD) [Li, et al. 2007] and even nonlinear ones [Habib, et al. 2015] have been developed to damp a single mode.

However, the tuning requirements of chatter suppression in machining processes are quite different. The repeating surface pattern on the workpiece causes regeneration that often leads to self-excited vibrations in a given cutting process. This effect can drive the cutting process to loss of its stability leading to large amplitude chatter vibration, which clearly must be prevented. The effect of this regeneration was recognized by pioneers in the field like Tlusty and Tobias [START_REF] Tlusty | Self-Excited Vibrations on Machine Tools[END_REF][START_REF] Tobias | [END_REF]. They demonstrated that the minimum stability of a certain process limited by a single mode is directly related to the real part of the frequency response function (FRF) of the system in the cutting point. Mathematically, regenerative chatter is described by delay differential equations DDE [Hale et al., 1977]. Chatter vibrations can grow when the dynamic stiffness of the machine/tool/part system is lower than the cutting stiffness of the cutting process [Merrit, 1965]. Therefore, the higher is the damping of the system is, the stronger stability against chatter vibrations is [START_REF] Munoa | Active suppression of structural chatter vibrations using machine drives and accelerometers[END_REF].

The application of TMD dampers to increase the damping of the machine tool was one of the first proposed solutions for chatter avoidance [START_REF] Koegnisberger | Machine Tool Structures[END_REF]. Some experimental studies were performed to find the ideal tuning of TMDs. For instance, [START_REF] Koegnisberger | Machine Tool Structures[END_REF] were proposing to build an electrical circuit equivalent to the mechanical system in order to optimize the real part of the FRF. [Rivin, et al. 1992] dealt with the improvement of stability in boring bars considering TMD. Also as an experimental work [Tarng, et al. 2000] and [Rashid, et al. 2008] tuned the natural frequency of their TMD to match with the natural frequency of the structure target mode of the structure. Finally, [Sims, 2007] found an analytical expression for the optimal tuning of a single TMD for chatter avoidance. The damper was tuned to maximize the negative real part of the main frequency response function (FRF) creating two equal peaks. The parameters of a viscoelastic cantilever beam was tuned in [START_REF] Saffury | Optimized chatter resistance of viscoelastic turning bars[END_REF]] by maximizing the most negative real part of the corresponding FRF. [Yang, et al. 2010] optimizes the MTMD for one dominant mode considering the real part of the FRF.

Recently, new self-tuneable TMD have been proposed based on variable stiffness concepts [START_REF] Aguirre | Self-tuning dynamic vibration absorber for machine tool chatter suppression[END_REF]. These VSTMD dampers offer the possibility to change the tuning of the damper depending on the cutting conditions. The aim of this study is to find the ideally best tuning to suppress chatter vibrations by taking the advantage of the variable values of the spring of a VSTMD concept (see Figure 1a). The optimal performance of the VSTMD is to be comparable with the results obtained with the best tuning of a single TMD.

In practice, the design of a VSTMD is not an obvious task. In fact, several concepts have been proposed in the literature [START_REF] Aguirre | Self-tuning dynamic vibration absorber for machine tool chatter suppression[END_REF]. For instance, [START_REF] Aguirre | Self-tuning dynamic vibration absorber for machine tool chatter suppression[END_REF] created this damper using a movable mass guided by means of a flexures and a rotary spring controlled by a stepped motor. Changing the angular position of the rotary spring, the natural frequency and therefore, the tuning of the damper, could be modified. In parallel, pure viscous damping was introduced in the system by means of the eddy current effect (see Figure 1).

The paper is divided into three parts, which include model development, the determination of the optimal tuning for a rotary spring self-tuning damper using the method of zeroth order approximation (ZOA) [START_REF] Altintas | Analytical prediction of stability lobes in milling[END_REF], Zatarain, et al. 2004] and the verification part. Finally, the concept is verified by using time domain simulations and by using semi-discretization method (SDM) [START_REF] Insperger | Semi-discretization for time-delay systems: stability and engineering applications[END_REF]. The model presented in the paper represents an idealistic situation, which needs further more development to fit a in real case scenario.

STABILITY MODEL FOR MILLING OPERATIONS WITH A VSTMD

A simple milling model is introduced here in order to describe the dynamic behavior of a VSTMD concept. The model is adequate for applying semi SDM and ZOA with and without taking time averaging in the tooth passing period of the milling operation. The ZOA permits to obtain a fast analytical parametric solution with a good accuracy for continuous milling processes. In the case of interrupted cutting, the presence of double period chatter and mode couplings reduces the precision of the ZOA [Munoa, et al. 2009].

The VSTMD is especially interesting for structural chatter cases in heavy duty operations. In these cases, tools with a high number of flutes and large engagements are used [Iglesias, et al. 2014]. [Den Hartog, 1934] (blue) and [Sims, 2007] (red).

Therefore, the periodic excitation is close to "flat" (almost constant) and the system can be roughly considered roughly as a time independent autonomous system. This way, the ZOA can also serve acceptable results, and simple assumptions can be made regarding the stability limits (see Figure 2a).

The first step to optimize the VSTMD is to develop a milling model that includes the effect of a passive damper. For simplicity, it is assumed that the damper is located at the point where the cutting force is acting. Two different points can also be considered with modal parameters [START_REF] Munoa | Chatter suppression in ram type travelling column milling machines using a biaxial inertial actuator[END_REF] but this increases the complexity of the main expressions and does not help to explain the main innovation in the present work. A single dominant mode parallel to the cutting plane has been considered, too. Finally, a lead angle of 90º has also been assumed. With this set of assumptions, a 2 DOF model has been formed to describe the stability of the milling process.

By assuming linear cutting force characteristics, the resultant cutting force can be derived in the following way for the simple straight fluted milling operation

)) ( ( )) ( ( )) ( ( ) ( c e 1 p t h t t g a t i i i Z i i i K K T F        , (1)
where a p and Z are the depth of cut and the number of teeth, respectively. The cutting force is taken into account by using K e = col(K e,t , K e,r ) edge coefficients and K c = K c,t κ c (κ c = col(1,  r )) cutting coefficients [Altintas, 2000]. A ratio  r has been defined

between the radial and the tangential force components.

The instantaneous chip thickness is a function of the chip thickness direction s i (t) = col(sin φ i (t), cos φ i (t)) of the i th tooth:

) ( cos )) ( ) ( ( ) ( sin )) ( ) ( ( ) ( t t y t y t t x t x f t h i i Z i            , (2)
where x(t) and y(t) are the relative spatial vibration of the center of the milling tool compared to its secondary motion described by the feed per tooth f Z . The instantaneous chip thickness is considered at the angle

Z i i t t π 2 ) 1 ( ) (      (3)
for each ith teeth. The milling tool rotates with Ω = 2 π n / 60 angular spindle frequency.

The cutting force is transformed back from local (tr) to tool (xy) system by T

i (φ i (t)) in (1) considering g i (φ i (t)
) screen function [Dombovari, et al. 2010b]. This takes into account the radial immersion by the entrance (φ en ) and exit angles (φ ex ).

The dominant mode direction is considered under the angle ψ relative to the y axis (see Figure 2a), which results in x(t) = q 1 (t) sin ψ and y(t) = q 1 (t) cos ψ. By means of the modal mass m 1 , damping c 1 and stiffness k 1 of the relevant mode, the one degree of freedom (DOF) milling system has governing equation in the form

)) ( ), ( ; ( ) ( ) ( ) ( 1 1 T 1 1 1 1 1 1      t q t q t t q k t q c t q m F p    with          cos sin p . (4)
Assuming that the passive damper has the same vibration mode defined by angle ψ, the equation of motion ( 4) is modified to

)) ( ), ( ; ( ) ( ) ( ) (      t t t t t t q q Q q K q C q M    , (5)
where

       2 1 0 0 m m M ,           2 2 2 2 1 c c c c c C ,           2 2 2 2 1 k k k k k K ,        0 ) ( ) ( T t t F p Q .
Linearizing around the periodic stationary solution (forced vibration) q p (t) = q p (t + T) (T = 2π/Ω/Z), the following linearized system can be derived for the perturbation u(t) = q(t)  q p (t):

)) ( ), ( ; ( ) ( ) ( ) (       t t t t t t u u Q u K u C u M    , (6)
where

                   0 )) ( ) ( ( ) ( ) ( ) ( ) ( 1 1 p c, 2 1  t u t u t B a K t Q t Q t t Q (7)
and

p s T p             ) ( )) ( ( )) ( ( ) ( T c 1 T t t t g t B i i i Z i i i    . ( 8
)
A time dependent factor B(t) has been obtained in the expression of the perturbed force. This coefficient is named time domain directional factor and it is a periodic function that collects the projections of the cutting force generated by the different flutes onto the direction of the single vibration mode and the projection of the generated vibration onto the chip thickness direction (s i ) [Zatarain, et al. 2010].

CALCULATION METHODS

This time-periodic directional factor can be averaged to obtain a zero order approximation (ZOA) that permits an analytical solution. This approximation is going to be used to investigate the optimal tuning of the different TMDs.

Since ZOA gives back the borders of stability, the SDM is used to calculate exact stability properties on the exact model of milling, which, although it is hardly usable for analytical calculations.

Non-smooth time domain simulation was used on the fly-over [Stepan, et al. 2011] DDE model of milling to determine the real chatter frequency which is defined as the dominant frequency of the large amplitude threshold vibration.

These three methods have been combined to explore the potential of the VSTMD.

Time averaged model for ZOA

In order to find an optimal tuning for the passive damper, analytical ZOA is used, which is based on the time averaged model of the milling operation. Due to the periodicity of the directional factor, it can be developed in Fourier series. ZOA is obtained when this periodic function is approximated by the average value. Milling operation can be approximated really well by ZOA if the operation is near to full immersion (not interrupted) milling. By the introduction of an equivalent cutting coefficient [START_REF] Opitz | Investigation and calculation of chatter behaviour of lathes and milling machines[END_REF] ) , , , ( π 2 :

ex en 0 0 r Z B       , ( 9
)
the periodicity is averaged out from the describing milling model. The averaged directional factor is

)))), ( 2 sin( )) ( 2 sin( ) ( 2 ( )) ( 2 cos( )) ( 2 (cos( 4 1 ) , , , ( en ex en ex en ex ex en 0                            r r (10)
The force variation ( 6) is in the direction of q 1 and can be written as

)) ( ) ( ( π 2 : ) ( 1 1 0 p c, 1       t u t u a K Z t Q t .
(11)

From [START_REF] Tlusty | Self-Excited Vibrations on Machine Tools[END_REF] to [Altintas, 2000], this leads to the following conditions for the border of asymptotic exponential stability

)) ( Re( π ) ( 11 0 c, lim    H K Z a t   and   ) ( 2 π ) 1 2 ( 1 ) (         k k , (12)
where k = 1, 2, … and H 11 (ω) is originated from the equation defined at (6) as

) ( ) ( ) (    Q H U   , (13) and )
] ( [ ) i ( ) ( 1 2     j i H       K C M H , if i, j = 1, 2, while ψ (ω)[  π, 0] is the phase of H 11 (ω).
This classical approach relates the stability of the cutting process to the value of real part of the FRF [START_REF] Mancisidor | Receptance coupling for tool point dynamic prediction by fixed boundaries approach[END_REF]. Therefore, equations at (12) open the way for using simple optimization based on the analytical representation of the FRF H 11 (ω) supposing c 1  0 [Sims, 2007]. The stability of regenerative cutting process requires a special tuning where the real part should be shrunk. The average directional factor β 0 can be positive or negative value depending on the engagement and the cutting characteristics. Therefore, if the directional factor is positive, the negative side of the real part should be maximized creating two equal peaks. On the other hand, if the directional factor is negative the positive side should be minimized.

  m c  1 2 :    f 2 2 2 m k   1 :    g 1 1 1 m k   Table 1; definition of dimensionless dynamic parameters.
The real part of a FRF dominated by a single mode can be given by means of several dimensionless parameters (see Table 1):

2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 11 ) 1 ( 4 ) ) )( 1 (( ) 1 ( 4 ) ) )( 1 )(( ( 1 )) ( Re( g g f g f g f g g g g f g f g f g g f k g H                     ,( 14
)
where f is the tuning of the absorber, μ is the mass ratio, g is the dimensionless frequency and χ = c 2 / c cr is the damping ratio relative to c cr = 2 m 2 ω 1 "critical" damping defined in [Den Hartog, 1934].

Linear stability of stationary cutting determined by SDM

The SDM is used in this work to obtain the linear stability borders. The complete description of this method can be found in the literature [START_REF] Insperger | Semi-discretization for time-delay systems: stability and engineering applications[END_REF]. It compiles the discrete representation of the solution operator (step matrix) of autonomous (time independent) DDEs. The step matrix B l is derived for one single arbitrary chosen time step Δt = Δθ = τ max /m. Time periodicity is taken into account by the use of Floquet theory [Farkas, 1994] over the tooth passing period T by simply multiplying the corresponding step matrices, that is for

l l r l r l r l z B B B z ... 2 1        
, where

) ) 1 ( ( col 1       k t l m k l x z . ( 15
)
In equation ( 15) r = T/Δt represents the resolution of the approximation, while θ[ τ max , 0] is the relative "delayed time" to describe a state of the delayed system considering the following solution definition x(t + θ) := col(u(t + θ), u  (t + θ)). The eigenvalues of the so-called Floquet multipliers determine if the resulting periodic stationary solution is stable: | μ k | < 1. The multiplications at (15) can be accelerated by using spare matrices or using special techniques like in [Henninger, et al. 2008] or determining efficiently the eigenvalue problem itself [START_REF] Zatarain | Stability analysis of milling with irregular pitch tools by the implicit subspace iteration method[END_REF]].

Time domain solution

The final verification of the optimal tuning has been performed using pure time domain simulation. The simulation of the milling process was performed by using the original milling model as a base, (5) is supplemented by the non-smooth behavior of the fly-over effect. Fly-over appears if the local instantaneous chip thickness

) ( cos )) ( ) ( ( cos ) ( sin ))) ( ) ( ( sin ( ) ; ( 1 1 1 1 t r t q t q t r t q t q f r r t h i i i i Z i i i              (16)
is less than zero considering all possible previous cuts. Considering all possible fly-overs and cuts, the effective chip thickness can be determined as a minimum h ie (t):= min l h i (t; l ). If h ie (t)  0 the corresponding ith tooth flies-over at the time instant t and the cutting force switches off. This results in the following force definition after (1)

)) ( ( )) ( ( )) ( ( ) ( c e 1 fo p t h t t g a t ie i i Z i i i K K T F        , ( 17
)
where the switching function g foi = 1, if the corresponding ith edge is in radial immersion and if that edge is not in fly-over state, that is, h ie (t) > 0. Formulating the first order form using the dynamics from (5) with the excitation at (17), the simulation was performed by standard dde23 introducing special initial function (IF). The IF is the corresponding stationary solution y 0 (θ

) = col( ) ( ), ( p p   q q 
) in order to avoid the effect of weak attraction zones appearing in the non-smooth fly-over system [Dombovari, et al. 2010].

VARIABLE STIFFNESS TUNEABLE MASS DAMPER (VSTMD)

In the literature, there are two main analytical results for the ideal tuning of passive dampers. On the one hand, the work of [Den Hartog, 1934] uses two invariant analytical points to achieve the decrease of the receptance magnitude in all frequencies. On the other hand, the work of [Sims, 2007] shows three invariant points to achieve the best possible real parts of the corresponding FRF considering a lim by ZOA (12) as an objective.

A case study described by [START_REF] Aguirre | Self-tuning dynamic vibration absorber for machine tool chatter suppression[END_REF] has been used to compare the performance of the different strategies. In this case, the moving mass was m 2 = 7 kg and the mass ratio is μ = 4.7 %. The optimal tuning f (14), the optimal stiffness k 2 and damping c 2 can be determined for both methods by considering initial process parameters presented in Table 2 and the optimal parameters described in Table 3. The results are presented in Figure 2b about a near full immersion milling process (see Table 2) without any passive damper, a passive damper tuned according to [Den Hartog, 1934], and a passive damper tuned according to [Sims, 2007].

However, all of these theories work with fixed and invariant stiffness of the damper (see k 2 at (5)) for the different spindle speeds, which is not the case for VSTMD. This tuning procedure and the design of the rotary spring (see Figure 1a) allow setting the stiffness value iteratively or even continuously between a minimum and a maximum stiffness values, which means different optimal tuning can be realized along the stability limits.

Frequency ratio

Damping ratio [Den Hartog, 1934] [Sims, 2007] (

   1 1 DH f 3 ) 1 ( 8 3      [Sims, 2007] (β 0 > 0) 2 2 n opt, ) 1 ( 2 2 2          f ) 1 ( 8 3 :     
β 0 < 0) 2 2 p opt, ) 1 ( 2 2 2          f
Table 3; Optimal tuning analytical expressions.

By the modification of ( 6), the equation of motion has the form

)) ( ), ( ; ( ) ( ) ( ) (       t t t t t t j u u Q u K u C u M    , ( 18
)
where

          j j j j k k k k k , 2 , 2 , 2 , 2 1 K , ] , [ max , 2 min , 2 , 2 k k k j  . ( 19
)
Here j represents the jth set stiffness of the iteration described in [START_REF] Aguirre | Self-tuning dynamic vibration absorber for machine tool chatter suppression[END_REF]. The iteration is only performed if unstable stationary cutting appears, that is, when the chatter frequency ω ch and its modulations raise up-to a previously defined level compared to the harmonics of the forced vibration. In this case, the embedded VSTMD system acts by rotating the specially designed spring. It tunes the system to the dominant chatter frequency if it is possible without the violations of the bounds (19). The idea is behind to set the "resonant" stable pockets (see Figure 2b) to the used spindle speed zone by changing the stiffness of the tunable damper to k 2, j = ω ch 2 m 2 .

Ideal tuning of VSTMD

Regardless of the tuning iteration procedure applied on the VSTMD, it can only be successful if the system can be stabilized at all. This means that there is a special stiffness value where the cutting process (stationary cutting) is stable. This fact is convenient because the ideal tuning can be found by means of linear theories presented at ( 12) and ( 15). The ZOA and the FRF at ( 14) can be used for the analytical derivation. In order to ease the notation the magnifying function

M 11 (g, f ) : = H 11 (g) k 1 is introduced.
According to the derivation of ZOA [Altintas, et al. 2008] the arising chatter frequency ω ch can be considered as sampling on the corresponding FRF. Consequently, considering non-interrupted case the FRF can be optimized at each dimensionless frequency g = ω ch / ω 1 . This means finding an extremum of Re(M 11 (g, f )) by taking dRe(M 11 (g, f ))/df, which leads to five roots for tuning including one at f 0 = 0 and four others symmetrically placed. Thus, these latter four can be considered as two extremums as f + and f -

2 2 2 3 2 2 4 6 2 ) 1 ) 1 ( ( ) 1 ) 1 ( )( 1 ( 2 ) 1 ( ) (               g g g g g g g f . (20)
This new formula shows the best tuning of a damper to have the maximum stability for a certain spindle speed represented indirectly (see ( 12)) by the dimensionless value (g). This formula can improve the result of Sims" proposal for a single spindle speed.

In Figure 3, Sims" locked frequencies are g p (f ), g n (f ) and g i (f ) = f and their limits are g plim and g nlim . Among the solutions presented at (20), it is always f + (g) that gives the maximum and f  (g) that gives the minimum for Re(M 11 (g, f )). Keeping in mind the results of the ZOA depending on the sign of the directional factor β 0 > 0 (β 0 < 0) at (10) Re(M 11 (g, f)) needs to be maximized (minimized). Namely, the best solution for the tuning is always given by f + (g) (f -(g)), except between the values of g plim and g nlim where the optimum values determined at (20) are extremes.

a) b)

Figure 4; Real parts of the optimally tuned magnifying function for positive directional factor (β 0 > 0) a), and the same for negative directional factor (β 0 < 0) compared to Sims' optimal tuning (red) [Sims, 2007].

In this region, depending on β 0 > 0 (β 0 < 0) a minimum f min (maximum f max ) tuning is limiting not realizable stiffness values on the rotary spring (see ( 21) and ( 22)).

In Figure 3, one can realize that Sims" constant optimal tuning f opt,p (f opt,n ) for β 0 > 0 (β 0 < 0) is only optimal at two dimensionless frequencies g, below and above g = 1. Although, along the lobes limited frequencies are possible, usually g>~1 (g<~1) if β 0 > 0 (β 0 < 0). This means one point along the lobe is optimal for a VSTMD in this simple case (see Figure 4).

Ideal stability behavior of VSTMD

Based on the calculated optimal tuning function f  (g), piecewise smooth definition can be defined, which take into account the physical limitation of the rotary spring. The limitations are only active close to the original natural frequency ω 1 of the structure.

When

β 0 > 0, } ), ( { max ) ( min p opt, f g f g f   , (21)
while

β 0 < 0, } ), ( min{ ) ( max n opt, f g f g f   . (22)
Using these piecewise smooth definitions, the stability boundaries can be depicted (see 5) by simply applying ZOA shown at (12). For both cases of negative and positive directional factors, the stability is improved by using VSTMD solution. In these ideal examples, one can realize there is a point in the stability limits where the stability is not improved with respect to the Sims" optimal tuning values (see Figure 3). It can be also noticed that the "double lobe" shape is disappears and in this ideal solution the lobes are simple shaped stability boundaries. Figure 5; ZOA stability limits based on the optimal FRF for positive directional factor (β 0 > 0) a), and the same for negative directional factor (β 0 < 0) b) compared with Sims' solutions (red).

Iterative algorithm for optimal tuning VSTMD

It has been demonstrated that if the tuning of the damper is changed for each spindle speed using a VSTMD, the stability of the process can be improved more than an ordinary TMD with the same mass tuned according to Sims" parameter. However, it is a complex procedure to implement it in practice. In this section, an iterative method to approach this optimal solution is going to be proposed.

The first step towards chatter suppression is to detect whether chatter is occurring or not during the machining process. The vibration measured by the accelerometer installed on the structure of the damper is processed in order to find its main frequency components. Here it is important to distinguish between forced vibrations, induced directly by the cutting forces, and chatter, which is due to an unstable regenerative process generated only under certain working conditions.

Forced vibrations appear at harmonics of the tooth passing frequency, but are stable, and thus are usually not a problem for machining, except in finishing operations where surface roughness needs to be improved. Chatter appears at other frequencies than tooth passing frequencies. It is an unstable cutting process, meaning that the cutting forces and vibrations increase with time, leading to unacceptable machining conditions, since they produce very bad surface quality and can lead to damage in the machine.

The chatter detection and suppression algorithm is presented in Figure 6 Figure 6; Chatter detection and suppression algorithm.

This algorithm is implemented on a real-time controller. It is running continuously during the machining process, calculating the spectrum of the measured vibration of the machine, as shown in Figure 7b. The algorithm detects the frequency of the maximum vibration peak, and compares it with the tooth passing frequency: if it is an integer multiple of the tooth passing frequency, it is considered to be a forced vibration, and no corrective action is taken. If it is not an integer multiple, it is considered to be chatter, and the angular position in the damper is modified in order to tune it to the chatter frequency.

It is very important to clearly distinguish clearly chatter from forced vibrations, so that the damper is only tuned to chatter frequencies. Otherwise, once the damper is tuned to the chatter frequency, the vibration level at this frequency will drop, and the algorithm will detect a forced vibration as main frequency. If the damper is tuned to this new frequency, chatter generation could start again, which needs to be avoided.

The detection algorithm first filters the acquired acceleration signal with a low pass filter. Then the FFT of the filtered signal is calculated, finding the frequency of the highest peak. When the highest peak"s frequency does not match any harmonic of the spindle speed, chatter detection is considered positive. In this section, real case simulations are presented by the real tuning procedure of the modelled VSTMD solution. Time domain simulations are only performed if the technological parameters set result in unstable cutting determined by SDM (15). The time domain simulations are started from the perturbed stationary solution q p to reach threshold fly-over [Stepan, et al. 2011] effect fast. Roughly 30-40 periods are simulated and 10-20 periods long of signal from the end of the simulation is taken for FFT analysis. With the removal of the DC component, the chatter frequency can be extracted easily since the stationary cutting solution q p is known and the harmonics of its amplitude can be used as a threshold for chatter. The maximum among the peaks overtaking the limiting amplitude is considered as chatter frequency and the VSTMD is tuned to this main frequency with k 2, j = ω ch 2 m 2 . This tuning is done subsequently until stability determined by SDM or a preset maximum number of iterations are reached.

Sims" and the ideal tuning solutions are presented in Figure 7a for milling operation calculated by SDM in this case. Time domain calculations were performed for selected spindle speeds (A, B and C) listed in Table 2 in order to see the tuning iterations of the VSTMD. The time domain simulations combined with the fly-over effect must be performed because there is no sufficiently mature method that can determine the dominant frequency of the threshold periodic, quasi-periodic or completely chaotic motion [Dombovari, et al. 2010].

The "originally stable" cases in Figure 7 and Table 4 are stable according to Sims", no tunings are necessary. The "tuned" cases are parameters, when the tuning was effective by using the simulated chatter frequencies. The "not tuneable" situations are the cases when even the VSTMD was not effective. =118.8, 113.5, 121.9, 111.5 Hz Stable in four tuning iterations C 4 : n=2500 rpm, a p,4 =70 mm ω ch,j =119.8,113.5,124.0,113.5,123… Hz Unstable Table 4;The detailed results of the time domain simulations of the selected points (see Figure 7a).

a) b)

A 1 (originally stable) B 1 (originally stable) C 1 (originally stable)

A 2 (tuned) B 2 (originally stable) C 3 (tuned)

A 3 (not tuneable) B 3 (not tuneable) C 4 (not tuneable) The specific parameters and tuning data of the simulated cases presented in Figure 7 can be followed in Table 4.

It is important to mention that the authors experienced small attraction zones at some parameter sets like C 3 when the non-smooth fly-over effect introduced weak attraction zone. In this case the accurate selection of IF was of great importance, although this situation hardly appears in realistic cases.

Summarizing, we can say that the time domain simulations confirm that VSTMD can improve the results obtained by a single VSTMD with the optimal design for chatter, and the proposed iterative algorithm can approach the optimal solution defined in the previous sections for a VSTMD.

CONCLUSION

In the presented work a variable stiffness tuneable mass damper concept was investigated. The real system is built as an embedded system that is capable to decouple the applied damping and stiffness.

The paper deals with the idealized tuning of this semi-active damper. Sims" idea to use the real part of the corresponding frequency response function was used when but in this case the optimization was performed at each frequency. This way, a new tuning method has been proposed to maximize the stability in a certain spindle speed. Using this concept, the frequency dependent tuning function was derived analytically and the best possible stability for a VSTMD has been obtained. Later the concept was confirmed by using the mathematical model of the tuned mass damper performing semidiscretization and time domain simulations.

An iterative method to approach the optimal solution has been proposed and verified by means of time domain simulations. It can be seen that the method together with the VSTMD solution is able to improve the stability behavior by simply retuning the damper to the best optimal stiffness value in each spindle speed. Time domain simulations sometimes showed weak attraction behavior which effect must be checked more closely in the future.

Figure 1 ;

 1 Figure 1;Rotary spring a), the magnets producing Eddy current and their guiding system b). The assembled passive damper block c).

  Figure 2; a) sketch of the simple one dimensional milling model. Panel b) shows the effect of mass dampers (MD) by [DenHartog, 1934] (blue) and[Sims, 2007] (red).

Figure 3 ;

 3 Figure 3; Optimal tuning w.r.t. the dimensionless frequency g.

Figure 7

 7 Figure 7; a) shows ideal tuning solutions for VSTMD (red) performed by SDM (β 0 > 0) compared to Sims' case (black). b) shows time domain solution at selected points.

  partially supported by EU FP7 HIPPOCAMP project (608800/ FP7-2013-NMP-ICT-FOF), the Hungarian Scientific Research Fundation OTKA Grant no. K108779 and by the European Research Council under the European Union"s Seventh Framework Programme (FP/2007-2013) / ERC Advanced Grant Agreement n. 340889.

Table 2 ;

 2 Milling process parameters.


Cherrelle Eid
email: c.eid@tudelft.nl

Paul Codani
email: paul.codani@supelec.fr

Yannick Perez
email: yannick.perez@u-psud.fr

Javier Reneses
email: javier.reneses@iit.upcomillas.es

Rudi Hakvoort
email: a.hakvoort@tudelft.nl

Javier Reneses Guillén

Managing Distributed Energy Resources in a Smart Grid Environment: A review for incentives, aggregation and market design

Keywords: Tariffs, distributed generation, aggregation, incentives, decentralized systems Highlights

Introduction

Traditionally low voltage grids have been designed to transport electricity towards residential users for consumption. However, due to the increased penetration of distributed energy resources (DER), low voltage grids are increasingly used as carriers of bi-directional electricity flows. The penetration of DER such as distributed generation (DG), electric storage, electric vehicles (EVs) and demand response significantly affect the operations of distribution grids (Arriaga & Bharatkumar, 2014;[START_REF] Pudjianto | Virtual power plant and system integration of distributed energy resources[END_REF]). In Germany for example, the growth of Solar Photovoltaics (PV) reached a level of 38 GW installed in 2013 and affected grid stability in some local areas [START_REF] Von Appen | Time in the Sun[END_REF]. Large numbers of PV installations are noticed in The United States (US) within California, Arizona and Hawaii (Greentech Media & Solar Energy Industries Association, 2013). Other examples of DER rises are a significant growth of EVs in Norwaywhere EVs stood for 12.5% of new car sales in 2014and Californiawith almost 130.000 plug-in vehicles on the roads by the end of 2014and CHP in Denmark (International Energy Agency, 2005).

On one hand, this DER development is positive due to reductions in CO 2 emissions with sustainable DG, decreased use of transmission lines, increased self-consumption and the increasing independence of customers from central grid power (Alanne & Saari, 2006). However, regardless of those, DER is potentially problematic for grid stability and reliability due to congestion and voltage issues [START_REF] Eftekharnejad | Impact of Increased Penetration of Photovoltaic Generation on Power Systems[END_REF][START_REF] Walling | Summary of distributed resources impact on power delivery systems[END_REF]. These concerns are mostly posing effects on the distribution net which is under supervision of the DSO (in Europe) or central service utility (in the USA in some places) and system operator for maintaining the overall grid balance. The German example shows that due to local electricity over production at the sunny moments of the day, reliability of supply is endangered in distribution grids [START_REF] Epia | Global Market Outlook for Photovoltaics[END_REF][START_REF] Epri | The Intergrated Grid: Realizing the Full value of Central and Distributed Energy Resources[END_REF][START_REF] Von Appen | Time in the Sun[END_REF][START_REF] Yan | Operation strategies in distribution systems with high level PV Penetration[END_REF]. In France, realistic forecasts count on 450.000 Plug-in Electric Vehicles on the roads by 2020 [START_REF] Rte | Bilan previsionnel de l'equilibre offre-demande d'electricite en France[END_REF]; if this objective is reached, simultaneous charging of these EVs could stand for between 5 to 20% of the annual peak load [START_REF] Rte | Bilan previsionnel de l'equilibre offre-demande d'electricite en France[END_REF]. Similar issues are found in … (add examples of challenges DER penetration; EVs, CHP, PV etc.).

Current research describes issues of DER penetration from both a technical perspective and market design perspective, including the solutions that could be proposed. For example, [START_REF] Eftekharnejad | Impact of Increased Penetration of Photovoltaic Generation on Power Systems[END_REF] and [START_REF] Dang | Transformer Operating conditions under introduction of PV and EVs in an Eco-district[END_REF] discuss the impact of PV penetration on grid stability and the improving effects that storage would provide. An holistic approach of DER management has been briefly discussed for the Norway sector [START_REF] Ottesen | How Demand Side Can Provide Operational Flexibility To Power System Through a Holistic Aggregator Concept[END_REF]. A description of DER from a DSO perspective is provided in the report of Perez-Arriaga et al, describing the effects of DER on the financial position of the DSO [START_REF] Pérez-Arriaga | From Distribution Networks to Smart Distribution Systems : Rethinking the Regulation of European Electricity DSOs[END_REF]. Cossent described the way in which DSO charges should be set up to incorporate different DER (Cossent, 2013). More recently, evolving roles for the DSO with increasing penetration of DER and smart metering has been presented as well [START_REF] Edso | Response to CEER public consultation on the future role of the DSO[END_REF][START_REF] Evolvdso | Development of methodologies and tools for new and evolving DSO roles for efficient DRES integration in distribution networks[END_REF].

However, an analysis of individual DER sources, their technical characteristics and the limitations that those technical issues cause for system use and economic trading of electric flexibility is lacking. Therefore, an interesting issue that remains is the ability of different DER sources to provide grid flexibility services and relate those to suited (price) incentives. Effective market design takes into account the technical complexity of DER management, the economic possibilities for designing markets for such services and eventual roles that should involve managing those transactions. Consequently, this paper presents a review of distributed energy resources from a technical perspective, the related system requirements for flexibility and the required incentives for DER management like (dynamic) tariffs and contractual arrangements. Depending on system status and policy objectives, some arrangements might better serve system purposes than others. Both in a liberalized and a vertical integrated sector this work is of interest for policy makers, electricity suppliers, network managers and emerging actors like aggregators and Energy Service Companies (ESCOs).

This paper starts with a description of general changes in the electricity system in Section 2. Section 3 presents a review of the most common distributed energy sources and their technical characteristics. Section 4 presents an overview of centralized markets for flexibility trading. Next, Section 5 reviews incentives for DER management like tariffs, contracts and direct control. Lastly, Section 6 presents the discussion and Section 7 the conclusions and policy recommendations.

From Traditional to Smart electricity Systems

Worldwide, different developments in electricity sectors challenge the traditional centralized management of electricity systems. The increased penetration of renewable energy sources (RES), the distribution of electricity production, the penetration of distributed energy resources and the move towards smart-metering call for a different approach on electricity consumption and production.

Last years, supportive Feed-in-Tariffs in for example Germany incentivized the installation of small solar panels in the residential and commercial sector. In 2014 Germany had 38 GW capacity of Solar PV installed, with a large part, (more than 60%) located at low voltage levels [START_REF] Epia | Global Market Outlook for Photovoltaics[END_REF]. Other examples of rapidly developing residential solar PV segment are found in Belgium (where 1 out of 13 households are equipped with a PV system), Denmark, Greece and the UK (EPIA, 2014). Likewise, large numbers of PV installations are noticed in The United States (US) in California, Arizona and Hawaii (Greentech Media & Solar Energy Industries Association, 2013). Electricity generation is thus increasingly placed at the distribution grid level as an alternative of transmission grid level. This affects the distributed nature of electricity generation (Alanne & Saari, 2006).

Demand response is a term that refers to "the changes in electric usage by end-use customers from their normal consumption patterns in response to changes in the price of electricity over time, or to incentive payments designed to induce lower electricity use at times of high wholesale market prices or when system reliability is jeopardized" (Aghaei & Alizadeh, 2013) meaning the electricity demand is being activated to respond on for example a price trigger. Regarding the terminology used, in this paper, we focus on price triggers and other incentives that could affect the "response" of the operation of different DER, including demand response itself.

Important system benefits related to the operation of renewable energy sources like wind and solar are decreased emissions and economic efficiency due to low marginal cost. However, for system operation, RES increase risks because of unpredictable production patterns due to their highly intermittent character. Therefore RES require flexibility services like back up generation to supply for balancing needs of the non-supplied demand. Next to those traditional methods of system balancing, demand response and storage can potentially supply the system with flexibility services. Storage units are potentially beneficial for electric energy time-shift, power supply capacity and transmission congestion relief [START_REF] Eyer | Energy Storage for the Electricity Grid : Benefits and Market Potential Assessment Guide[END_REF].

Next to the previous named developments regarding the variability of RES generation, the distributed nature of generation and the change of demand from static to responsive, other developments affect the way in which distribution grids require possible more decentralized management. The first one relates to the electrification of transport with the electric vehicle (EV). The EV development is important because EV charging significantly increases electricity consumption at distribution grids during peak periods, potentially jeopardizing security of supply due to congestion and voltage issues [START_REF] Clement-Nyns | The impact of vehicle-to-grid on the distribution grid[END_REF][START_REF] Green Emotion | D4.3 -B2: Grid Impact studies of electric vehicles ; Reinforcement Costs in Low-voltage Grids[END_REF].

Distributed Energy Resources as flexibility service providers within electricity systems

As described in the previous section, Distributed energy resources (DER) e.g. Electric Vehicles (EVs), combined heat and power (CHP) units, electric water heaters and storage units are potentially providers of flexibility services. Technically speaking, an electric flexibility product can be defined as a power adjustment sustained for a given duration in order to balance supply and demand at a given moment in time. Thus, a flexibility service is a service characterized by four attributes (see Figure 1): its direction (upward or downward) (a); its electrical composition in power (b); and its temporal characteristics defined by its starting time (c) and duration (d). Some DER may have a single direction potential contribution (for instance typical household loads, such as water heaters, dishwashers and electric heaters), while others have bidirectional capabilities and could both act as a consuming and producing unit (e.g. EVs and storage units).

Furthermore, the electrical composition is of importance in order to state for what system flexibility needs the DER could serve, which calls for a differentiation between power and energy resources. The former have a rather low energy/power ratio. Those DER can provide the electricity system with a high power value, but are not able to maintain this power level for a long period of time. The latter have a high energy/power ratio and are more appropriate to maintain a change in power level for a longer period of time. The power resources are consequently better suited for short term markets (e.g. on the ancillary service markets) while energy resources are better suited for long term markets like for instance balancing mechanisms and trading DR in the bulk electricity market.

In order to compare the different DER on this criterion, we define the max power temporal ratio t r (expressed in time) as the maximum duration a DER can sustain its maximum power variation with respect to its nominal power. For some DER types, this parameter can be computed by dividing the allowed energy range by the maximum power capacity (e.g. considering a stationary battery with a charging/discharging power equal to 10kW and an energy capacity equal to 50kWh, we find t r = 5h). For some other DERs, it may be related to physical characteristics (for instance for a water heater with thermic inertia, we may find t r = 30min). The lower this value, the more the DER can be considered as a capacity type DER, and vice-versa. This variable is intended to provide insights on differences between DER categories and is not a singular value for all DER in such category, simply because this is dependent on the specific technology.

What is more, the availability (in time) is a constraint that distinguishes the moment at which DER could provide services to the system. Some resources may only be available during specific periods of timefor instance EVs are most likely to be available from 6 PM to 6 AM. In order to compare the flexibility providing units on this criterion, we compute for each of the DER the ratio a r defined as the average number of hours during which the unit is available divided by the total number of hours in a week. Also here, depending on the specific source (one EV differs from the other for example), in reality similar DER may offer different availability times, however we aim to provide insight in expected values.

Furthermore, the activation time refers to the aspect that some resources may be able to adjust their power much quicker than other sources. Generally, almost all electric appliances have a fast activation time, ranging from the order of a second to one minute, except for CHP units which have longer ramping times [START_REF] Houwing | Demand Response With Micro-CHP Systems[END_REF]. Lastly, the geolocation of DER is of importance for the supplied nature of the required demand response. For example, geo-location based demand response could be of interest for local congestion management or distributed generation (DG) optimization. Table 1 provides an overview of the most common DER and their characteristics. The table is divided in different types of DER; consumption, bi-directional and generation.

Consuming DER: residential loads

New generation LED lightings could adapt their power consumption to required grid power variations [START_REF] Lee | A design methodology for smart LED lighting systems powered by weakly regulated renewable power grids[END_REF]. Future LED systems could undergo system power variations up to 35% while humans would only perceive a variation of 15% in light intensity. This method would be particularly interesting for public lighting. On the contrary, older lighting systems do not have the same abilities [START_REF] Lee | A design methodology for smart LED lighting systems powered by weakly regulated renewable power grids[END_REF][START_REF] Samarakoon | Investigation of domestic load control to provide primary frequency response using smart meters[END_REF]. Changing their power consumption would have serious impact on their luminous capability. LED lighting systems can maintain this power variation for a 1 Where DS stands for distribution grid and TS for transmission grid significant period of time and therefore can be considered as energy type flexibility resources. However, their potential power modulation is relatively low in absolute values. Their predictability is relatively good, while their availability highly depends on the usage considered. Typical lightings would be turned on from a few hours a day during peak hours to 12 hours a day, thus we have 0.2 < a r < 0.5. It is noticeable that this criterion is highly seasonal dependent.

DER

Residential appliances, such as water heaters, washing machines, electrical heaters, air conditioners etc., have rather low max power temporal ratios t r (many of those units are off or stand-by most of the time). The latter can range from a few seconds (e.g. for cookers) to a dozen of minutes (electric space heaters) [START_REF] Samarakoon | Investigation of domestic load control to provide primary frequency response using smart meters[END_REF], thus providing a maximum temporal ratio of 5s < t r < 15min.

Their availability depends a lot on the appliance considered: whereas electric space heaters have a good availability, therefore named continuous loads here (0.4 < a r < 1) while washing machines have a very limited one (a r < 0.1). Similar rationale applies for their predictability [START_REF] Tomiyama | Modeling air conditioner load for power system studies[END_REF][START_REF] Wong | Demand response potential of water heaters to mitigate minimum generation conditions[END_REF]. Heat pumps coupled with thermal energy storage stand out in this category: indeed, their max power temporal ratio can reach up to 3h without any inconvenience for end-users [START_REF] Arteconi | Domestic demand-side management (DSM): Role of heat pumps and thermal energy storage (TES) systems[END_REF], making those units suited for longer grid services such as peak shaving.

Bi-directional DER: Electrochemical storage and EVs

Storage units are potentially beneficial for electric energy time-shift, power supply capacity and transmission congestion relief [START_REF] Eyer | Energy Storage for the Electricity Grid : Benefits and Market Potential Assessment Guide[END_REF]. Electrochemical Energy Storage (EES) units have a perfect availability and predictability (a r ≈ 1). Whether they should be considered as energy type or power type resources depend on their power density and energy density characteristics, both parameters being much related to the type of battery technology (Li-ion, Ni-MH, Ni-Cd, etc.) [START_REF] Yang | Electrochemical energy storage for green grid[END_REF]. Thus, it is possible to find EES units for all kind of applications, from very-fast highpower responding units (such as supercapacitors, t r ≈ 4s) to energy type chemical batteries (such as Li-ion batteries, t r ≈ 10h).

Most Electric Vehicles that are on the roads today have a battery capacity around 20kWh2 . Their max power temporal ratio depends on the power of the charging station they are plugged in. Typical charging station powers range from 3kW to 50kW, leading to approximately 30 min < t r < 6 h. Because EVs are primarily used for transportation, capacity type services that would not empty the battery should be encouraged. Privately owned EVs are mainly available during nighttime and weekends (a r ≈ 0.5), but the availability could rise up to a r > 0.9 if charging points are installed at working places. Company fleets have slightly different usage patterns and could also be available in the afternoons (a r ≈ 0.8). Their predictability patterns are easily foreseeable [START_REF] Pearre | Electric vehicles: How much range is required for a day's driving?[END_REF], especially considering a fleet of EVs and not a single vehicle.

Producing DER: Micro CHP and PV units

Micro-CHP units are small heat and electricity generating entities. They have a large availability and predictability since they are dedicated to heat and electricity production (a r ≈ 1). It is more difficult to define a max power temporal ratio for micro-CHP units because they could produce electricity at maximum power continuously, as far as they are being supplied in primary energy source (mainly gas). Rather, their availability to maintain a change in their electricity production will be based on economics considerations. The control strategies of micro-CHP units are likely to take into account energy costs [START_REF] Houwing | Demand Response With Micro-CHP Systems[END_REF] in their economic balance. Therefore, micro-CHP units would fit in the energy type category.

PV units are different from the others, in the sense that their production output cannot be controlledonly, they can be curtailed in a few minutes in case of extreme network conditions. They are available in daytime, so between 6 hours and ten hours a day depending on their location. Good predictability is difficult to achieve, and is dependent on the method used and the desired spatial resolution. Generally speaking, we can say that acceptable forecasts can be achieved a few hours a-head (International Energy Agency, 2013).

4

Markets for electric flexibility trading

Traditional electricity systems are managed in a top-down manner, meaning that generally large generation units connected at high voltage levels feed in electricity for electricity consumers that are located at all other voltage levels. Flexible generation units (mostly hydro units, gas and coal fired plants) are besides providers of bulk electricity supply, also providers of electric flexibility by means of upward and downwards adjustments. Those adjustments could be incentivized by for example capacity contracts with the SO for automatic adjustments.

Beside generation, also consuming units might be suppliers of electric flexibility. France and the United Kingdom (UK) are important frontrunners regarding demand response developments in Europe (SEDC, 2014). In France, already before sector liberalization, demand response activity was triggered by the electricity utility EDF for industrial electricity customers. These units received dynamic tariffs that incentivized consumption shifting. Table 2 provides an overview of the most common traditional markets for electricity trading in the short and long term, based on the French system. The next sections describe in further detail the presented examples.

DER trading for ancillary services

Ancillary service markets are in place in order to manage transaction for system balancing in the short to very short term. These markets are organized very close to real-time and require automated load adjustment. In France ancillary service markets are organized shorter than 30 seconds before real-time for Frequency Containment reserves (FCR, also named primary reserve), below 15 minutes for frequency restauration reserves (FRR, also named secondary reserve) and lastly replacement reserves (RR, tertiary reserve) for system balancing between 13-2h before real time (see Table 2). In the United States (US) and UK different projects present examples of DER flexibility provision within ancillary service markets [START_REF] Dynamicdemand | A dynamically-controlled refrigerator[END_REF][START_REF] Kempton | A Test of Vehicle-to-Grid ( V2G ) for Energy Storage and Frequency Regulation in the PJM System[END_REF]. Due to the fact that individual DER do not provide sufficient reliable electric flexibility to be tradable in markets, aggregation is required in order to trade in organized markets. The US example of EV trading is the Delware project where an EV aggregator acts as an intermediary firm between PJM (the local SO) and flexibility service providing EVs. This aggregator sells a certain amount of capacity to the grid operator and bids this in the hourly auction for frequency regulation and for the available power capacity ($/MW-h) [START_REF] Kempton | Public policy toward GIVs for TSO services: State of the art and future trends[END_REF][START_REF] Kempton | A Test of Vehicle-to-Grid ( V2G ) for Energy Storage and Frequency Regulation in the PJM System[END_REF]. When participating in this frequency regulation, EVs receive a dispatch signal from the local TSO (PJM) and are remunerated accordingly.

If the regulation service offered by the Delaware EV aggregator has not met with the performance thresholds over a specified time period in terms of correlation (delay) and precision, PJM is able to disqualify the aggregator [START_REF] Chris | PJM manual 12: Balancing Operations[END_REF].

DER trading for system balancing and network congestion management

Balancing services are arranged a bit longer before real-time than ancillary services and in different places in Europe they are open for aggregation and demand response. In Germany many industrial loads are directly participating in the balancing mechanism, however for aggregated loads still many barriers exist to participate within the balancing markets [START_REF] Koliou | Demand response in liberalized electricity markets: Analysis of aggregated load participation in the German balancing mechanism[END_REF]. In the French system such barriers have been lowered by reduction of the minimum bidding capacities for balancing services from 50 to 10 MW. This has been done in order to motivate the entrance of smaller entities like aggregators to participate in balancing mechanisms (SEDC, 2014).

Differently, for network congestion management a French example of small load aggregation is the aggregator named Voltalis3 . Customers contracted with Voltalis receive a free box installed in their home named Bluepod, which reduces their electric heating device operation in short time intervals when Voltalis receives a signal from the TSO. The dispatch signal is mostly related to endangered electricity supply sufficiency in Brittany (a French region poorly interconnected to the other ones) and network limitations. Customers who have the box installed are automatically enrolled, but can opt-out at any time by pushing a button on the device and use their electric heater as usually. Voltalis as an aggregator is able to trade the aggregated flexibility in different markets like balancing markets and demand response mechanisms of the TSO. The customers do not receive any financial benefit when their heating device reduces their load, but observe a reduction of their normal electricity bill due to those interruptions in electricity consumption for heating.

In Sweden the DSO can incentivize load shifts by the provision of Time of Use (TOU) prices in order to deter network investments or decrease congestion by incentivizing the customer to shift the load away from peak moments [START_REF] Bartusch | Further exploring the potential of residential demand response programs in electricity distribution[END_REF][START_REF] Bartusch | Introducing a demand-based electricity distribution tariff in the residential sector: Demand response and customer perception[END_REF]. Different from the previous examples, the DSO does not trade this flexibility within a market for congestion management or deterred network investments, but this is a direct incentive arrangement between the DSO and electricity users.

DER trading in spot markets and generation capacity markets

As the first one in Europe, the French system provides a possibility for demand response trading within spot markets. This is possible since 2014 wherein demand response can be traded in the dayahead market through the NEBEF mechanism4 . In 2017 it is foreseen that demand response will also be tradable in capacity markets in France (RTE, 2013a). Furthermore the French TSO organizes an annual tender dedicated to specifically to Demand Response providers. Other possibilities of spot market and capacity trading markets for DR? Check.

Incentives for efficient operation of Distributed Energy Resources

As the distribution system transitions from a passive network of consumers to a more actively managed system of network users with diverse consumption and production behaviors, price signals will play a crucial role in shaping the interactions between the physical components of the distribution system and network users (Arriaga & Bharatkumar, 2014). There are different ways to create efficient interaction in electricity systems; tariffs are one of the incentive methods that are being described in this section.

A one-sided economic approach on settlement of incentives for DER management would be the application of real-time nodal pricing that would incorporate both grid and supply constraints at each moment in time, incentivizing upward or downward adjustments for all DER [START_REF] Sotkiewicz | Nodal Pricing for Distribution Networks : Efficient Pricing for Efficiency Enhancing DG[END_REF]. However, due to the fact that each flexibility source has its own technical requirements and abilities to provide flexibility services, non-singular approach; a combination of tariffs, contracts and direct control.

Broadly speaking, a distinction is made between price based and controllable methods for demand response [START_REF] Pfeifenberger | Demand response review. … for the Alberta Electric System Operator[END_REF], also referred to as price based and interruptible demand response [START_REF] Muratori | Role of residential demand response in modern electricity markets[END_REF] or as direct and indirect methods of load modification. Next to tariffs therefore, direct control and other contract arrangements are methods by which efficient operation of DER could be incentivized.

Price based methods for DER management

Price-based demand response is incentivized by exposing the DER user to a time-varying electricity rate, also called a dynamic rate. The theory of dynamic tariffs for demand response has already been discussed in 1989 by David and Lee for large industrial electricity users [START_REF] David | Dynamic tariffs: theory of utility-consumer interaction[END_REF]. Table 3 presents an overview of those tariff options with definitions. In this table a distinction is made between basic dynamic pricing options and those which specifically incentivize adjustments of users' normal consumption patterns (also called baseline consumption adjustments). The basic pricing options leave most freedom to the user, without requiring extra information on baseline consumption levels. Options for such pricing methods are 1) Time-of-Use pricing (TOU), 2) Real-Time Pricing (RTP) and 3) Critical Peak Pricing (CPP). The more specific incentives for baseline adjustments are 4) Peak Time Rebates (PTR), 5) Interruptible capacity programs (ICAP) and 6) Emergency demand response [START_REF] Newsham | The effect of utility time-varying pricing and load control strategies on residential summer peak electricity use: A review[END_REF]. Those options require baseline consumption information penalizing or remunerating for specific load adjustments. With RTP, the user receives a changing price per time step (for example 15 minutes) and the customer will shift electricity consumption accordingly. With critical peak pricing, only in specific hours per day a higher price is presented to the customer. Electricity customers receive an ex-ante notice of these moments in time and can therefore plan their consumption [START_REF] Koliou | Development of Demand Side Response in Liberalized Electricity Markets : Policies for Effective Market Design in Europe[END_REF]. Critical peak pricing toghether with the options for baseline adjustments are specifically incentivizing the shift of electricity consumption away from a specific moment in time. A driver for such incentives could relate to the for example high wholesale market prices or jeopardized system reliability [START_REF] Koliou | Development of Demand Side Response in Liberalized Electricity Markets : Policies for Effective Market Design in Europe[END_REF].

Basic Dynamic Pricing Options Specific incentives for baseline adjustment

Direct load control for DER management

With controllable or incentive-based demand response, a central actor like the system operator, aggregator or even retailer could make the end-user agree to automatically curtail its electricity demand under certain circumstances. This means that a central actor has direct access to the load and is able to reduce or increase this as required for the system or for portfolio management purposes. Load shedding refers to the "switching off" of entire network zones from electricity supply in order to sustain total system operation [START_REF] Newsham | The effect of utility time-varying pricing and load control strategies on residential summer peak electricity use: A review[END_REF]. With brown outs, the system operator slightly reduces voltage frequency in order reduce the needed electricity transport capacity and generation capacity but to maintain electricity supply quality within limitations [START_REF] Blume | Electric power system basics: for the nonelectrical professional[END_REF]. The named approaches are more contractual and introduce constraints for the supplied flexibility, while the price-based approaches described in Section 5.1 leave the customer to decide in real-time regarding the supply of flexibility (DOE, 2006). Consequently, direct control methods are probably more suited for short term services, or services which require a very precise location of activation like voltage control and congestion management. Table 4 provides an overview of different incentives presented in Section 5 and relates them to their suitability to DER types and markets for flexibility.

Techno-economic alignment of incentives

Depending on the type of consumer and the DER installed, a specific contract or price incentive might be appropriate to incentivize DER interactions that are efficient for the overall system. With focus on the technical activation time of DER and possible incentives, Table 4 provides an overview of possible arrangements. For grid interactions which require response between 1 to 30 minutes before real-time direct load control would be suited in order to secure reliability of supply. Appropriate DER for such short notification time periods would be most DER except for CHP units due to their longer response times. For longer notification times of 30 min to 1 hour all other pricing methods could be suited and decisions should be further depend on socio-technical issues like the urge of required response, the direction and user characteristics like price elasticity and the availability of home automation for example. All DER types would be appropriate for supplying flexibility for longer than 1 hour of activation time. For the very long term, critical peak pricing and time of use pricing is appropriate due to the possibility to settle those prices on a yearly basis, as the case in France with the tempo tariff 7 . Differently, PV units provide in some way possible flexibility but in combination with storage unit their reliability would be improved.

Discussion

This paper provided an overview of different DER and their technical abilities to provide flexibility services to the central grid. The activation of the demand side is technically bounded by the technical abilities, but furthermore should take into account other variables which are discussed here.

Socio-economic effects on appropriate incentives

The presented work provides a good starting point of deciding between appropriate DER for flexibility management in electricity systems. However, in reality, appropriate incentives are not only technical-economic dependent, but also socio-economic affected by individual factors like consumption behavior, sustainability drive, relative price differences and price elasticity of the user [START_REF] Goulden | Smart grids, smart users ? The role of the user in demand side management[END_REF]. This work did not look into those issues but pilot projects which takes into account those variables could provide more insight in those issues and could furthermore specifically direct which incentive would be most appropriate depending on consumer and DER characteristics.

Transitioning towards decentralized system operation

Beside the socio-economic context, also the regulatory environment of the electricity system at stake will affect the decisions for deciding on appropriate incentives. Flexibility trading options presented in & Saari, 2006;[START_REF] Kamat | Multi-settlement Systems for Electricity Markets : Zonal Aggregation under Network Uncertainty and[END_REF][START_REF] Pudjianto | Virtual power plant and system integration of distributed energy resources[END_REF]. For example, distribution congestion management and distribution balancing are not possible within centralized market designs and till now these externalities have been integrated in the final electricity price. Consequently, DER penetration could call for parallel a decentral arrangement of trading in markets which could enable flexibility trading for electricity flows from different voltage levels and time periods. Beside the fact that decentralized management would yield benefits for more cost-causality based incentives, it is discussed by Goulden et al. that moving away from centralized system towards decentralized system furthermore encourages a new approach on consumption, moving from "passive energy consumer" towards an "active energy citizen" [START_REF] Goulden | Smart grids, smart users ? The role of the user in demand side management[END_REF]. Therefore the use of current markets for DER management might be a transition phase possibly towards a decentralized techno-economic management approach of the electricity system. Figure 2 presents a conceptual presentation of the arrangement of such a decentralized system based on arising system challenges and opportunities with DER integration. Depending on the electricity market design and the level of sector liberalization, one or more of the actors in the sector could decide on tariffs, control and other incentives. Therefore insight in the role(s) of the DSO, electricity retailer, supplier, (independent) aggregator and other involved third parties are crucial in order to understand the process of incentive settlement. Challenges that arise have not been dealt before with, for example the ones related to load aggregation. Due to the fact that there are minimum trading values for the balancing and other free markets, DER should be bundled to simultaneously provide significant tradable amounts of flexibility in those other markets. However, challenges arise when this is done by independent actors due to competing plans of electricity suppliers and demand response aggregators [START_REF] Eurelectric | Designing fair and equitable market rules for demand response aggregation[END_REF].

7

Conclusions and Policy Recommendations

This paper presented a review of existing Distributed Energy Resources' (DER) abilities to provide flexibility services and options to incentivize such interactions. With a central management approach on the electricity system, flexibility services from DER could be traded within traditional centrally arranged markets for flexibility for securing reliability of supply. A one-sided economic approach on settlement of incentives for DER management would be the application of real-time nodal pricing that would incorporate both grid and supply constraints in order to incentivize efficient interactions for upward or downward adjustments at each moment in time for every DER. However, due to the fact that each flexibility source has its own technical requirements and abilities to provide flexibility services, the authors of this paper argue that DER require a non-singular approach; a combination of tariffs, contracts and direct control.

Next to central management of DER, also decentralized management of DER could be possible through for example local aggregators and/or virtual power plant arrangements. The potential for this type of management is arising, especially with the increasing roll-out of smart meters, sensing and control at distribution level and upcoming risks due to over-voltage and congestions with the penetration distributed generation (DG). Decentralized management of the system could reduce local needs for expenses in the network and support the reliability of supply with high level of DG. Consequently, the options presented in this paper are suited for a transition phase of centralized towards more decentralized operation, but further decentralization might require new market design. Therefore future research should take into account such options of market design and the roles of traditional and new actors within such markets. Depending on the current and expected challenges within electricity systems, policy should anticipate the required DER transactions and design markets that incentivize arrangements that will benefit the system from an economic, sustainability and reliability perspective.

Acknowledgements

Figure 1 :

 1 Figure 1: Characterization of flexibility products[START_REF] Eid | Aggregation of Demand Side flexibility in a Smart Grid : A review for European Market Design[END_REF]

Figure 2 :

 2 Figure 2: Conceptual presentation of relationships between electricity system needs and decentralization 6.3 Settlement of incentives and control: which roles for different actors?

Table 1 : The different DER and their technical characteristics

 1

Flexibility direction Flexibility characteristic (power vs energy) Availability ratio

						Predictability	Technical response time	Grid 1	Ref.
									(Lee, Li, & Hui,
									2011; Lu, Xie,
	Electrical Consumption	dishwasher) Lighting loads (W) Dispatchable, residential loads (washing machines,	Unidirectional Unidirectional	New LED systems: energy types older lightings: power types Power type 5s < t r < 5min	20% < a r < 50% during peak hours a r < 0.1 low max power ratios t r due to max off time	Good High	Second Second	DS DS	Huang, Puyleart, & Yang, 2008; Samarakoon, Ekanayake, & Jenkins, 2012) (Lu et al., 2008; Samarakoon et al., 2012)
		Electrical							(Samarakoon et
		heating/ Cooling (continuous	Unidirectional	Power type t r ≈ 15min	0.4 < a r < 1	High	Second	DS	al., 2012; Tomiyama, Daniel, & Ihara,
		loads)							1998)
	Bi-directional	Electrochemi cal Energy Storage (EES) (kW-MW) Electric Vehicle	Bidirectional Unidirectional or	Power & Energy types 4s < t r < 10h Power & Energy types	a r ≈ 100% 50% < a r < 90%	Perfect High	Second to Minute Second	DS or TS DS	(Divya & Østergaard, 2009; Yang et al., 2011) (Kempton et al., 2009; Pearre, Kempton,
		(kW)	Bidirectional	30 min < t r < 6 h					Guensler, &
									Elango, 2011)
	Generation	PV Unit Micro-CHP unit (kW)	Unidirectional (production mode) Unidirectional (production mode)	Not flexible Energy type	25%< a r <41% a r ≈ 100%	Good a few hours ahead Perfect	Minute Rather slow (0.15kW/min)	DS DS	(International Energy Agency, 2013) (Houwing et al., 2010)

Table 2 : Traditional markets for balancing from a single node perspective (based on the French case study)

 2

Table 3 : Possible dynamic pricing options for DER management (David & Lee, 1989; Faruqui & Sergici, 2009; Hakvoort & Koliou, 2014)

 3

Table 4 : Relationship of notification times, appropriate incentives and markets for DER flexibility trading

 4

  Table 2 are all presented in the framework of centralized system management, generally managed a system operator. However, as discussed by others, decentralized management approach could open up possibilities for locational pricing, local balancing and optimization at community level with virtual power plant (VPP) approach (Alanne

	Notification time before real-time	Appropriate incentives or control method for DER management	Related markets for electric flexibility trading 8	Appropriate DER
	< One minute	Direct load control (DLC)	Frequency control (primary, secondary, tertiary reserves), voltage control	EV, Continuous loads (heating/cooling, lightning), EES
	1-15 minutes	Direct load control	Network restauration, voltage control	EV, Continuous loads (heating/cooling) , EES
	15-30 min	Direct load control	Network restauration (HV/LV), Balancing market, Portfolio balancing	EV, EES, CHP units Continuous loads (heating/cooling), dispatchable loads
	1 hour	ICAP, Emergency demand response, Real time pricing, Direct load control, Peak time rebates, Critical Peak Pricing	Balancing market	EV, EES, CHP units Continuous loads (heating/cooling), dispatchable loads
		ICAP, Emergency demand		EV, EES, CHP units
	1 -48 hour	response, Real time pricing, Direct load control, Peak time rebates, Critical Peak	Spot Market (Day ahead and Intraday market)	Continuous loads (heating/cooling), dispatchable loads, PV
		Pricing		units with storage
				EV, EES, CHP units
	Year ahead	Critical peak pricing, Time of use pricing	Deferring network investments (HV/LV), generation investment peak reduction	Continuous loads (heating/cooling), dispatchable loads, PV
				units with storage

7 See http://residential.edf.com/energy-at-home/offers/electricity/tarif-bleu-56121.html 8 Composed with insight from report

(International Energy, 2008)

by

  Cherrelle Eid has been awarded an Erasmus Mundus PhD Fellowship and would like to express her gratitude towards all partner institutions within the program as well as the European Commission for their support. Paul Codani benefits from the support of the Chair "PSA Peugeot Citroen Automobile: Hybrid technologies and Economy of Electromobility", the so-called Armand Peugeot Chair led by CentraleSupélec, ESSEC Business School and sponsored by Peugot Citroen Automobile. The authors would like to thank Andrea Benets for her input regarding electricity market design in France. 8 References Aghaei, J., & Alizadeh, M.-I. (2013). Demand response in smart electricity grids equipped with renewable energy sources: A review. Renewable and Sustainable Energy Reviews, 18, 64-72. Alanne, K., & Saari, A. (2006). Distributed energy generation and sustainable development. Renewable and Sustainable Energy Reviews, 10(6), 539-558. Arriaga, I. P.-, & Bharatkumar, A. (2014). A Framework for Redesigning Distribution Network Use of System Charges Under High Penetration of Distributed Energy Resources New principles for new problems DRAFT FOR PUBLICATION . DO NOT DISTRIBUTE .

Nissan Leaf: 24kWh; Renault Zoe: 22kWh; BMW i3: 19kWh

Information on Voltalis via: www.voltalis.fr

See https://clients.rte-france.com/lang/fr/clients_distributeurs/services_clients/effacements.jsp

Note that these values relate to the French system and are very system dependent.

This paper focuses mainly on DER connected at distribution level, however demand response from large industrial units which is mostly connected at high voltage levels can be also seen as a part of DER and therefore taken into account here.