Zoltan Dombovari
email: dombovari@mm.bme.hu

Jokin Munoa
email: munoa@ideko.es

Rachel Kuske
email: rachel@math.gatech.edu

Gabor Stepan

THE EFFECT OF SLOWLY VARYING DYNAMICS ON MILLING STABILITY

Keywords: milling, regeneration, chatter, varying dynamics

In order to predict the quality and the stability properties of milling processes, the relevant dynamics reduced to the cutting edges needs to be known. However, this dynamics varies through the workspace along the tool path during a given machining operation. This is the case for large heavy duty milling operations, where the main source of the relevant dynamics is related to the otherwise slowly varying machine structure rather than to the fairly steady milling tool dynamics. The effect of slowly varying dynamic parameters is presented on milling stability when the cutting process takes place in a region of the working space where the steady-state cutting would change from stable to unstable. After the separation of the slow and fast time scales, the governing non-autonomous delay differential equation is frozen in slow-time in order to determine the timeperiodic stationary cutting solution of the milling operation for different ram extensions.

Introduction

The aim of this work is to point out the effect of slowly changing parameters on milling dynamics. It is well known that machining processes like drilling, turning or milling are subjected to regenerative effect when the past relative motion of the workpiece-tool system influences the present behavior of the operation. By modelling the geometric arrangement of the cutting edges, the corresponding regenerative delays can be identified [START_REF] Munoa | Chatter suppression techniques in metal cutting[END_REF]. Combining with an empirical cutting force charasteristics and with the dynamic model of the machine tool structure, the governing equation can be derived [START_REF] Altintas | Generalized dynamic model of metal cutting operations[END_REF][START_REF] Dombovari | General milling stability model for cylindrical tools[END_REF]. This results in delay differential equation (DDE) [START_REF] Hale | Theory of functional differential equations[END_REF] of autonomous (time-independent) or non-autonomous (time-dependent) kinds for different machining operations. In case of milling, the governing equations are time-periodic due to the nonregular cutter-workpiece-engagement (CWE) and the rotation of the tool [START_REF] Ferry | Virtual five-axis flank milling of jet engine impellers-part I: Mechanics of five-axis flank milling[END_REF]. In this time-periodic case, the instability of corresponding time-periodic stationary solution refers to unstable milling operation that leads to chatter [START_REF] Otto | Extension of Tlusty's law for the identification of chatter stability lobes in multi-dimensional cutting processes[END_REF]. By using the Floquet theory on the linearized variational system [START_REF] Farkas | Periodic Motions[END_REF], stability charts can be constructed usually in the parameter space of spindle speed n and depth of cut a. Between the stable and unstable domains, the stability boundaries correspond to either (secondary) Hopf or period doubling losses of stability.

The above mentioned methodology is capable to predict chatter for constant parameters; however, in reality, one or more parameters may be slowly varying during the machining operations. For example, large machines are well known to have varying dynamic behavior, thus, slowly moving cutter through their workspace is subjected to slowly varying dynamic properties. In five axis milling, even rigid compact machines operate in slowly changing environment during complex 3D tool motions, while the varying geometry along the tool path also affects the CWE in time.

In mathematical terms, the slowly changing variable introduces a permanent non-cyclic time dependency in the originally time-periodic milling model. This means that the DDE cannot be handled using Floquet theory, or at least, not in a straightforward manner. The slowly changing ordinary differential equation (ODE) models have already revealed the effect of slowly changing parameters on the corresponding stability loss and bifurcation [START_REF] Haberman | Slowly-varying jump and transition phenomena associated with algebraic bifurcation problems[END_REF][START_REF] Erneux | Imperfect bifurcation with a slowly-varying control parameter[END_REF]. In case of Hopf bifurcations, by using slow time scale, a shift of the stability boundary can be identified by considering the accumulative effect of the variational dynamics around the slowly changing stationary solution [START_REF] Baer | The Slow Passage Through a Hopf Bifurcation: Delay, Memory Effects, and Resonance[END_REF].

This work intends to apply the slow time scale methodology in the time-periodic DDE model of milling operations to show how those slowly varying parameters affect the classical chatter predictions. The paper considers a simple 1Dof model of the milling operation as a demonstrative example. After having a mathematical assumption for the slow-time deviation from the frozen-parameter case, a first order delay differential equation is derived for the amplitude and phase deviations. Stability of the slowly changing dynamics is determined by considering the exponential growth of the amplitude deviation from the frozentime solutions. The single degree of freedom model describes 'naively' the cantilever-like structural arrangement of a heavy-duty milling machine. The model represents the slowly changing dynamics during milling processes performed in ram-axis-direction. It is shown that ram-directional motion can shift stability boundaries only because of the varying dynamics.

1DoF milling model

where 𝜔 n , 𝜉, 𝜏 and 𝐾 c,𝑡 stand for the natural angular frequency, damping ratio, regenerative delay and tangential cutting coefficient, respectively. The regenerative delay is originated in the tooth passing frequency Ω 𝑍 = 𝑍Ω as 𝜏 = 2π/Ω 𝑍 , where Ω (rad/s) = 2π𝑛 rpm / (60 s/min) is the angular velocity of the tool. The milling cutter is assumed to be equally spaced with cutting edges, that is, the principle period equals with the regenerative delay as 𝑇 p = 𝜏 (see Figure 1a). The stability diagram is usually depicted in the parameter space spanned by the spindle speed 𝑛 and axial depth of cut a as demonstrated in Figure 1b. As shown in [START_REF] Iglesias | Analytical expressions for chatter analysis in milling operations with one dominant mode[END_REF] and [START_REF] Zatarain | Stability analysis of milling with irregular pitch tools by the implicit subspace iteration method[END_REF], the time periodicity of (1) appears in

𝐵(𝑡) ≡ 𝐵(𝑡 + 𝑇 p ) = 𝐔 ⊺ ∑ 𝑔(𝜑 𝑖 (𝑡)) sin 𝜅 𝐓(𝜑 𝑖 (𝑡))𝛋 c ⊗ 𝐧(𝜑 𝑖 (𝑡)) 𝑍 𝑖=1 𝐔, (2)
and in the stationary excitation

𝐺(𝑡) ≡ 𝐺(𝑡 + 𝑇 p ) = -𝑎 𝐔 ⊺ ∑ 𝑔(𝜑 𝑖 (𝑡)) sin 𝜅 𝐓(𝜑 𝑖 (𝑡)) (𝐊 e + 𝐊 c ⊗ 𝐧(𝜑 𝑖 (𝑡)) [ 𝑓 𝑍 0 0 ]) 𝑍 𝑖=1 . (3)
The mass normalized mode shape vector, lead angle and the feed per tooth are denoted by 𝐔, 𝜅 and 𝑓 𝑍 . The edge normal, the edge coefficients and the cutting coefficients are stored in 𝐧, 𝐊 e and 𝐊 c = 𝐾 c,𝑡 𝜿 c vectors defined in local edge (t, r, a) system, while the transformation between (t, r, a) and (x, y, z) is realized by 𝐓(𝜑) [START_REF] Dombovari | The effect of serration on mechanics and stability of milling cutters[END_REF]. The edge angular position and the effect of CWE are traced by 𝜑 𝑖 (𝑡) and 𝑔(𝜑) [START_REF] Dombovari | The effect of serration on mechanics and stability of milling cutters[END_REF].

Linear stability of the time periodic stationary solution 𝑄(𝑡) = 𝑄(𝑡 + 𝑇 p ) of ( 1) can be easily determined by various methods, like multi-frequency solution [START_REF] Merdol | Multi Frequency Solution of Chatter Stability for Low Immersion Milling[END_REF][START_REF] Budak | Analytical Prediction of Chatter Stability in Milling Part I: General Formulation[END_REF][START_REF] Zatarain | Estudio comparativo de los modelos matemáticos de chatter en fresado: monofrecuencia, multifrecuencia y simulación en el tiempo[END_REF] and semi-discretization [START_REF] Insperger | Semi-Discretization for Time-Delay Systems: Stability and Engineering Applications[END_REF]. The presented example with 90 deg lead angle 𝜅 and zero helix angle 𝜂 (see Table 1) was calculated by first order semidiscretization and using triangularization algorithm.

Table 1: Parameters of full immersion milling process used in this paper [START_REF] Dombovari | Prediction of multiple dominant chatter frequencies in milling processes[END_REF].

𝑍 𝜅 (deg) 𝜂 (deg) feed direction mode direction 4 90 0 [1 0 0] [0 1 0] 𝐾 c,𝑡 (MPa) 𝐾 c,𝑟 (MPa) 𝐾 c,𝑎 (MPa) 𝜔 n (Hz) 𝜉 (%) 𝑘 (N/

Slowly changing milling model

A real milling process is usually subjected to time-dependent slowly changing parameters. In this manner, we can modify the DDE (1) by considering some slowly changing parameters with respect to the so-called "slow time" 𝑠 ≔ 𝜀𝑡, where the general rate of change is denoted by ε . Thus,

𝑦(𝑡) + 2 𝜉(𝑠)𝜔 n (𝑠)𝑦̇(𝑡) + 𝜔 n 2 (𝑠)𝑦(𝑡) = 𝑎(𝑠)𝐾 c,𝑡 𝐵(𝑡, 𝑠) (𝑦(𝑡 -𝜏(𝑠)) -𝑦(𝑡)) + 𝐺(𝑡, 𝑠). (4)
Surely, not all parameters are changing at the same rate during a given process. For example, in the general form of (4) one can model a cutting process defined for a ramp-like workpiece [START_REF] Shi | Theory of finite-amplitude machine-tool instability[END_REF] with slowly changing 𝑎 ≔ 𝑎(𝑠) and 𝐺(𝑡) ≔ 𝐺(𝑡, 𝑠) (see path A in Figure 1b and description in Subsection 4.1). In the meantime, the changing ram extension can be followed by slowly changing 𝜔 n ≔ 𝜔 n (𝑠), 𝜉 ≔ 𝜉(𝑠), 𝐵(𝑡) ≔ 𝐵(𝑡, 𝑠)

and 𝐺(𝑡) ≔ 𝐺(𝑡, 𝑠) due to 𝐔 ≔ 𝐔(𝑠) (see path B in Figure 1b and description in Section 4.2). Note that 𝑦 ≡ 𝑞; it has been introduced only to distinguish the solution of (1) and that of (4).

In any case, the governing DDE (4) form of the slowly changing milling operation is no longer exactly timeperiodic. It behaves as a general non-autonomous (time-dependent) system with two different time scales described by the (real) 'fast time' t and the (introduced) 'slow time' s.

Stationary Solution

It is straightforward to assume that (4) must have a slowly changing but in this case not time-periodic stationary solution 𝑌(𝑡) ≠ 𝑌(𝑡 + 𝑇 p ). Perturbation 𝑥 is introduced in the solution as 𝑦(𝑡): = 𝑌(𝑡, 𝑠) + 𝑥(𝑡).

(

) 5
Abusing the notation, the real time and slow time dependencies are dropped for a while. Substituting the assumption ( 5) into (4) one can get the following form

𝑌 ̈+ 𝑥̈+ 2𝜉𝜔 n (𝑌 ̇+ 𝑥) + 𝜔 n 2 (𝑌 + 𝑥) = 𝑎 𝐾 c,𝑡 𝐵(𝑌 𝜏,𝜀𝜏 + 𝑥 𝜏,𝜀𝜏 -𝑌 -𝑥) + 𝐺, (6)
where

𝑌 𝜏,𝜀𝜏 : = 𝑌(𝑡 -𝜏, 𝑠 -𝜀𝜏) ≈ 𝑌(𝑡 -𝜏, 𝑠) -𝜀𝜏 𝜕 𝜕𝑠 𝑌(𝑡 -𝜏, 𝑠) =: 𝑌 𝜏,0 -𝜀𝜏 𝑌 𝑠 𝜏,0 .
This results in

𝑌 ̈+ 𝑥̈+ 2𝜉𝜔 n (𝑌 ̇+ 𝑥) + 𝜔 n 2 (𝑌 + 𝑥) = 𝑎 𝐾 c,𝑡 𝐵(𝑌 𝜏,0 -𝜀𝜏 𝑌 𝑠 𝜏,0 + 𝑥 𝜏,𝜀𝜏 -𝑌 -𝑥) + 𝐺. (7)
Approximating the stationary solution 𝑌(𝑡, 𝑠) of ( 4) with the stationary solution of (1) for fixed (frozen) slow time s, that is, by 𝑌(𝑡, 𝑠) ≈ 𝑄(𝑡; 𝑠) = 𝑄(𝑡 + 𝑇 p ; 𝑠), ( 7) is simplified to

𝑥̈+ 2𝜉𝜔 n 𝑥̇+ 𝜔 n 2 𝑥 = 𝑎 𝐾 c,𝑡 𝐵(𝑥 𝜏,𝜀𝜏 -𝑥) -𝑎 𝐾 c,𝑡 𝜀𝜏 𝑄 𝑠 𝜏,0 . (8)
Considering the time periodicity in 𝑄 𝑠 𝜏,0 , (8) has the same form as [START_REF] Hale | Theory of functional differential equations[END_REF]. By following the same procedure, a frozen-time stationary solution can be determined as 𝑋(𝑡; 𝑠) herewith introduced by 𝑥(𝑡): = 𝜀 𝑋(𝑡; 𝑠) + 𝑢(𝑡).

In this successive substitution, actually, an expansion w.r.t. 𝜀 of the stationary solution of (4) can be derived, resulting in the following form

𝑢̈+ 2𝜉𝜔 n 𝑢̇+ 𝜔 n 2 𝑢 = 𝑎 𝐾 c,𝑡 𝐵(𝑢 𝜏,𝜀𝜏 -𝑢) -𝑎 𝐾 c,𝑡 𝜀 2 𝜏 𝑋 𝑠 𝜏,0 , (10)
with 𝑌(𝑡, 𝑠) ≈ 𝑄(𝑡; 𝑠) + 𝜀 𝑋(𝑡; 𝑠).

Considering 𝜀 to be sufficiently small, the last term in [START_REF] Baer | The Slow Passage Through a Hopf Bifurcation: Delay, Memory Effects, and Resonance[END_REF] can be dropped and the variational equation described by the perturbation u around the stationary solution 𝑌 has the form

𝑢̈+ 2𝜉𝜔 n 𝑢̇+ 𝜔 n 2 𝑢 = 𝑎 𝐾 c,𝑡 𝐵(𝑢 𝜏,𝜀𝜏 -𝑢). ( 12
)

Asymptotic behavior of the stationary solution

The asymptotic behavior of this non-autonomous slow/fast system can be derived by the WKB method [START_REF] Baer | The Slow Passage Through a Hopf Bifurcation: Delay, Memory Effects, and Resonance[END_REF] considering the original time periodicity slightly depending on the slow time s. According to Floquet theory [START_REF] Farkas | Periodic Motions[END_REF], the general solution of a linear time periodic system is given by an exponential and a time-periodic term [START_REF] Dombovari | Prediction of multiple dominant chatter frequencies in milling processes[END_REF]. Similarly, using the WKB method, the following general solution can be assumed for the slowly changing equation ( 12

In ( 13), the slow time s is also carried by the asymptotic behavior referring to the stability (rapid exponential growth σ(𝑠)) of the slowly varying solution u on the fast time scale t. Substitution of ( 13) into [START_REF] Zatarain | Stability analysis of milling with irregular pitch tools by the implicit subspace iteration method[END_REF]

However, keeping ε sufficiently small and consequently having small change on 𝐴 w.r.t. slow time s, one can assume 𝐴 𝑠 = 𝐴 𝑠𝑡 = 𝐴 𝑠𝑠 = 0. Thus, Fourier expansion is applied on the now exactly time-periodic 𝐴 and B as

𝐴(𝑡): = ∑ 𝐴 𝑙 ∞ 𝑙=-∞ e i 𝑙 Ω 𝑍 𝑡 and 𝐵(𝑡): = ∑ 𝐵 𝑙 ∞ 𝑙=-∞ e i 𝑙 Ω 𝑍 𝑡 . (16)
The multi-frequency approach [START_REF] Budak | Analytical Prediction of Chatter Stability in Milling Part I: General Formulation[END_REF] or Hill type of infinite expansion of the slowly changing milling dynamics can be given after substituting ( 16) into ( 14) as

((𝜀 𝜎 𝑠𝑠 + 𝜎 𝑠 2 + 2𝜉𝜔 n 𝜎 𝑠 + 𝜔 n 2 )𝐈 + [2 i 𝑙 Ω 𝑍 𝜎 𝑠 -𝑙 2 Ω 𝑍 2 + 2 𝜉 i 𝑙 Ω 𝑍 ] 𝑙=-∞ ∞ ) [ ⋮ 𝐴 -1 𝐴 0 𝐴 1 ⋮ ] = 𝑎 𝐾 c,𝑡 (e -τ 𝜎 𝑠 -1) [ ⋱ ⋮ ⋮ ⋮ ⋰ ⋯ 𝐵 0 𝐵 -1 𝐵 -2 ⋯ ⋯ 𝐵 1 𝐵 0 𝐵 -1 ⋯ ⋯ 𝐵 2 𝐵 1 𝐵 0 ⋯ ⋰ ⋮ ⋮ ⋮ ⋱ ] [ ⋮ 𝐴 -1 𝐴 0 𝐴 1 ⋮ ] . (17)
For the sake of simplicity, zeroth order (average) consideration can be derived by picking only the averages of 𝐴 and 𝐵 as in [START_REF] Altintas | Analytical Prediction of Stability Lobes in Milling[END_REF], resulting in

(𝜀 𝜎 𝑠𝑠 + 𝜎 𝑠 2 + 2𝜉𝜔 n 𝜎 𝑠 + 𝜔 n 2 -𝑎 𝐾 c,𝑡 (e -τ 𝜎 𝑠 -1)𝐵 0 )𝐴 0 = 0. (18)
In [START_REF] Shi | Theory of finite-amplitude machine-tool instability[END_REF], 𝐴 0 has non-trivial solution if the following is satisfied

𝜀 𝜎 𝑠𝑠 + 𝜎 𝑠 2 + 2𝜉𝜔 n 𝜎 𝑠 + 𝜔 n 2 -𝑎 𝐾 c,𝑡 (e -τ 𝜎 𝑠 -1)𝐵 0 = 0. ( 19
)
Note that, in [START_REF] Dombovari | Prediction of multiple dominant chatter frequencies in milling processes[END_REF], at least one parameter is slow-time dependent. Also, [START_REF] Dombovari | Prediction of multiple dominant chatter frequencies in milling processes[END_REF] resembles to the characteristic equation of constant parameter case [START_REF] Munoa | Chatter suppression techniques in metal cutting[END_REF]. Moreover, in point of 𝜎 𝑠 , [START_REF] Dombovari | Prediction of multiple dominant chatter frequencies in milling processes[END_REF] has the form of a singular perturbed algebraic equation Γ(𝜎 𝑠 ), which is released by dynamic term 𝜎 𝑠𝑠 with 𝜀, that is

𝜀 𝜎 𝑠𝑠 + Γ(𝜎 𝑠 ) = 0. ( 20
)
Obviously, if the rate of parameter change approaches zero, ( 19) is restricted to the constant parameter characteristic equation for 𝜎 𝑠 originated from (13).

Stability criteria

The general solution for the slow-time system in [START_REF] Dombovari | The effect of serration on mechanics and stability of milling cutters[END_REF] suggests that, in case of negative real part 𝜎, the perturbation introduced in (9) dies out, while positive real part σ induces rapid explosion w.r.t. s. That means, the first order nonlinear ODE representation of the characteristic equation in [START_REF] Dombovari | Prediction of multiple dominant chatter frequencies in milling processes[END_REF] has to be integrated by using a proper initial condition when s = 0. This can be done by using numerical algorithms like Runge-Kutta method. For the sake of initial condition, we assume at s = 0 the system is frozen for a while in slow time s, then

𝜎(𝑠) ≈ 𝜆𝑠 ⟹ 𝜎 𝑠 (𝑠) ≈ 𝜆 ⟹ 𝜎 𝑠 (0) ≈ 𝜆, ( 21
)
where 𝜆 is the frozen time characteristic exponent. The real part of the cumulated value of 𝜎 𝑠 over slow time has the form

re (∫ 𝜎 𝑠 (𝜁)d𝜁 𝑠 0 ) (22)
that refers to the stability properties of ( 17) as shown in [START_REF] Baer | The Slow Passage Through a Hopf Bifurcation: Delay, Memory Effects, and Resonance[END_REF]. If the function ( 22) crosses zero, the stability property of the slowly changing stationary solution 𝑌(𝑡, 𝑠) (11) will flip.

Keeping in mind the 'singular perturbation'-like structure of (20), formula [START_REF] Gatti | Applied Structural and Mechanical Vibrations: Theory, Methods and Measuring Instrumentation[END_REF] has interesting theoretical consequences. On the one hand, sufficiently small 𝜀 induces shifting of the onset of unstable motion independently of 𝜀, as shown by only integrating the solution of Γ(𝜎 𝑠 ) = 0 from [START_REF] Altintas | Analytical Prediction of Stability Lobes in Milling[END_REF]. On the other hand, [START_REF] Gatti | Applied Structural and Mechanical Vibrations: Theory, Methods and Measuring Instrumentation[END_REF] suggests that the accumulated stability is overtaken by the accumulated instability, which results in the shift on the onset of the unstable motion, which very much depends on the initial values of the slowly changing variables [START_REF] Baer | The Slow Passage Through a Hopf Bifurcation: Delay, Memory Effects, and Resonance[END_REF].

Case Studies

In this section we provide two distinct examples, in which the behavior of the dynamic bifurcation analysis might have a relevance. It is important to emphasize that these are completely artificial examples. Real case study might result in different significance of the explained effect.

Slowly varying depth of cut

In the literature there are various measurement examples, when the test required pre-manufactured workpiece with a gentle slope [START_REF] Shi | Theory of finite-amplitude machine-tool instability[END_REF][START_REF] Dombovari | Prediction of multiple dominant chatter frequencies in milling processes[END_REF][START_REF] Stepan | Identification of cutting force characteristics based on chatter experiments[END_REF]. These measurements typically are aimed to present stability limits or so-called nonlinear hysteresis phenomenon, the direct consequence of subcritical Hopf bifurcation of the stationary milling process. The slowly changing parameter in (4) the axial depth of cut is defined as

𝑎(𝑠) = 𝑎 min + 𝜀 𝑡, (23)
with 𝑎 min = 1 mm and 𝑎 max = 10 mm. In Figure 2ab) one can follow the accumulation of the critical eigenvalue that shows according to the criterion [START_REF] Gatti | Applied Structural and Mechanical Vibrations: Theory, Methods and Measuring Instrumentation[END_REF]. The system loses its stability after the constant parameter limit a c at the dynamic one a din . The loss of stability can be realized as the time domain solution in Figure 2c) "escapes" the stationary solution calculated by using simply the frozen time stationary solution 𝑄(𝑡; 𝑠). One can realize the linear dependency of the axial depth of cut a on the stationary solution in Figure 2c). It can be also realized in Figure 2c) the actual onset point where the solution escapes from the stationary solution is a bit ahead of the predicted position. This suggests deeper dependency of the rate of change 𝜀 on the dynamics which needs further, more detailed study of the problem.

Slowly varying ram overhang

Obviously, the presented relations in [START_REF] Dombovari | On the bistable zone of milling processes[END_REF] are far to be true in a real machine. Instead, the real behavior should be characterized by measuring frequency response functions (FRFs) in many overhang positions. In this manner the real dynamics and even the damping can be determined. This artificial example is calculated with

𝐿(𝑠) = 𝐿 min + 𝜀 𝑡, (25)
with 𝐿 min = 100 mm and 𝐿 max = 122.

In Figure 3ab) the accumulated effect of the slowly changing parameter L is shown on the critical eigenvalue. Similarly to Figure 2ab), the dynamic bifurcation limit L din and the constant parameter stability limit L c are not the same. A shift appears confirmed in Figure 3c) by time domain simulation. Here, the stationary solution changes drastically (blue and red envelope curves in Figure 3c) by varying the ram overhang L. The amplitude of the stationary solution shrinks because by changing the L the system goes away from the stability pockets (see path B in Figure 1b), where "resonance" causing high but finite gain on the amplitude [START_REF] Bachrathy | Surface properties of the machined workpiece for helical mills[END_REF]. Although the stationary solution shrinks, the system becomes more flexible causing the onset of the unstable motion after the constant parameter limit L c at L din and the amplitude goes to infinite (without modelling the threshold fly-over effect [START_REF] Dombovari | On the bistable zone of milling processes[END_REF]).

Conclusion

There are several industrial problems, including ones related to cutting technologies, which may involve dynamic processes on (very) different time scales. In the present work, we study milling operations that clearly have fast time-periodic dynamics. In the meantime, there exists a slow rate of change of some system parameters originated in the slowly varying structural dynamics as the tool moves in the working space of the milling machine.

In this paper, the so-called dynamic bifurcation phenomenon has been introduced for the analysis of milling stability. A new generalized governing equation was derived, with which the stability of the slowly changing milling dynamics can be predicted. Two simplified case studies presented slowly varying behavior. The nontrivial, sometimes counter-intuitive theoretical predictions based on the analysis of the new governing equations were confirmed by time domain simulation results, although some parameter domains still need further and deeper study.

The results may have industrial relevance when the milling cutter moving in the workspace has varying reduced dynamics. The results are somewhat counter intuitive. On the one hand, the accumulated stability, in theory, does not or weakly depends on the rate of change, which suggests that shifting of the stability appears even for extremely small rate of change. On the other hand, the shift (e.g., a din  a c ) carries the initial parameter value, since this has direct effect on the accumulated stability behavior.

It is important to emphasize that the above presented results apply only in cases when the cutting operation starts from stable stationary cutting; at this point, the results do not explain the transition backward from chatter to stable stationary cutting, which requires different modelling techniques that are applicable also for quasi-periodic oscillations.

Acknowledgement

The research was partially funded by the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013) ERC Advanced grant agreement No. 340889 and the scholarship of the Hungarian Academy of Sciences (BO/00589/13/6).

Figure 1 :

 1 Figure 1: a) sketch of the actual mechanical model of milling; b) dimensionless stability lobe diagram with dimensionless spindle speed υ and depth of cut α [11]. The paths A and B represent examples related to slowly increasing depth of cut a (see subsection 4.1) and ram overhang L (see subsection 4.2), respectively.Simple one DoF model of the regenerative milling process is considered here with the following second order time periodic DDE

  ): 𝑢(𝑡): = 𝑢(𝑡, 𝑠) = e σ(𝑠) 𝜀 𝐴(𝑡, 𝑠), where 𝐴(𝑡, 𝑠) = 𝐴(𝑡 + 𝑇 p , 𝑠).

Figure 2 .

 2 Figure 2. The real part of the critical eigenvalue 𝜎(𝑠) a) and its derivative 𝜎 𝑠 (𝑠) b) solving (17) w.r.t. s slowly changing axial depth of cut a(s). In panel c) time domain simulation y(t) and the frozen time stationary solution 𝑄(𝑡; 𝑠) are shown. Here 𝜀 = 1 × 10 -4 (m/s).

Figure 3 .

 3 Figure 3. The real part of the critical eigenvalue 𝜎(𝑠) a) and its derivative 𝜎 𝑠 (𝑠) b) solving (17) w.r.t. s under slowly changing ram extension (cantilever overhang) L(s). In panel c) time domain simulation y(t) and the frozen time stationary solution 𝑄(𝑡; 𝑠) is shown. 𝜀 = 1 × 10 -4 (m/s).Considering the ram as simply a cantilever beam the following dependencies can be determined w.r.t. the ram overhang 𝐿 by using simple Euler beam theory[START_REF] Gatti | Applied Structural and Mechanical Vibrations: Theory, Methods and Measuring Instrumentation[END_REF] 𝜔 n = 𝐴 𝜔 n 1 𝐿 2 , 𝑘 = 𝐴 𝑘 1 𝐿 , 𝐔 = 𝐀 𝐔 1 √𝐿 , and consequently 𝐵 0 ∼ 1 𝐿 . (24)

  𝐴 ̇= 𝐴 𝑡 + 𝜀 𝐴 𝑠 and 𝐴 ̈= 𝐴 𝑡𝑡 + 2𝜀 𝐴 𝑡𝑠 + 𝜀 2 𝐴 𝑠𝑠 .

					leads to a
	partial differential equation as
	𝜀 𝜎 𝑠𝑠 𝐴 + 𝜎 𝑠 2 𝐴 + 2𝜎 𝑠 𝐴 ̇+ 𝐴 ̈+ 2𝜉𝜔 n (𝜎𝐴 + 𝐴 ̇) + 𝜔 n 2 𝐴 = 𝑎 𝐾 c,𝑡 𝐵 (e -τ 𝜎 𝑠 -1) 𝐴	(14)
	by using the assumption e	σ(𝑠-𝜀 𝜏) 𝜀	≈ e	σ(𝑠) 𝜀 e -τ 𝜎 𝑠 (𝑠) . In accordance with (13), one can obtain


Kichun Jo
email: kichun.jo@gmail.com

Sungjin Cho
email: sungjincho215@gmail.com

Chansoo Kim
email: chansoo7857@gmail.com

Paulo Resende
email: paulo.resende@valeo.com

Benazouz Bradai
email: benazouz.bradai@valeo.com

Fawzi Nashashibi
email: fawzi.nashashibi@inria.fr

Myoungho Sunwoo
email: msunwoo@hanyang.ac.kr

Cloud Update of Tiled Evidential Occupancy Grid Maps for the Multi-Vehicle Mapping

Keywords: Dempster-Shafer theory, occupancy grid map, cloud service, multi-vehicle mapping, LiDAR

Nowadays, many intelligent vehicles are equipped with various sensors to recognize their surrounding environment and to measure the motion or position of the vehicle. In addition, the number of intelligent vehicles equipped with a mobile Internet modem is increasing. Based on the sensors and Internet connection, the intelligent vehicles are able to share the sensor information with other vehicles via a cloud service. The sensor information sharing via the cloud service promises to improve the safe and efficient operation of the multiple intelligent vehicles. This paper presents a cloud update framework of occupancy grid maps for multiple intelligent vehicles in a large-scale environment. An evidential theory is applied to create the occupancy grid maps to address sensor disturbance such as measurement noise, occlusion and dynamic objects. Multiple vehicles equipped with LiDARs, motion sensors, and a low-cost GPS receiver create the evidential occupancy grid map (EOGM) for their passing trajectory based on GraphSLAM. A geodetic quad-tree tile system is applied to manage the EOGM, which provides a common tiling format to cover the large-scale environment. The created EOGM tiles are uploaded to EOGM cloud and merged with old EOGM tiles in the cloud using Dempster combination of evidential theory. Experiments were performed to evaluate the multiple EOGM mapping and the cloud update framework for large-scale road environment.

Introduction

Intelligent vehicles, such as self-driving cars and industrial mobile robots, require surrounding environment information to perform efficient and safe driving. A behavior decision system in the intelligent vehicle uses the environment information to determine the optimal behavior in consideration of passenger comfort and energy efficiency. The environment information is also used to find the safest vehicle motion that avoids collisions with obstacles [START_REF] Wolcott | Fast LIDAR localization using multiresolution Gaussian mixture maps[END_REF][START_REF] Guo | A Low-Cost Solution for Automatic Lane-Level Map Generation Using Conventional In-Car Sensors[END_REF][START_REF] Joshi | Generation of Accurate Lane-Level Maps from Coarse Prior Maps and Lidar[END_REF][START_REF] Wolcott | Robust LIDAR localization using multiresolution Gaussian mixture maps for autonomous driving[END_REF]. An occupancy grid map is one way to model the driving environment by dividing space into grid (or voxel) cells and informing the occupancy state of each cell. The occupied states indicate whether the physical space of the cell is occupied or free of static objects.

Probabilistic approaches have been applied to infer the occupancy state of grid cells with consideration of sensor noise characteristics [START_REF] Thrun | A real-time algorithm for mobile robot mapping with applications to multi-robot and 3D mapping[END_REF][START_REF] Grisetti | Improved Techniques for Grid Mapping With Rao-Blackwellized Particle Filters[END_REF][START_REF] Schütz | Laser-Based Hierarchical Grid Mapping for Detection and Tracking of Moving Objects[END_REF].These papers provide the probabilistic grid mapping using a Bayesian approach using a LiDAR beam model. Each cell state of the probabilistic occupancy grid map (POGM) models the posterior probability of occupancy. Although the probabilistic approaches of POGM are able to manage the sensor noise well, there are limits to conditions of sensor blind spots and dynamic objects. For the sensor blind spots, the occupancy state is unknown because of no sensor measurement for the detection area. For the dynamic objects, the occupancy state is conflict because the sensor measurement changes over time. In the probabilistic approach, the unknown and conflict state are only represented by a probability of 0.5, which is not clear to represent the state of the cell.

To overcome the limits of the POGM, evidential (Dempster-Shafer) theory has been applied to model the occupancy grid map [START_REF] Smets | Analyzing the combination of conflicting belief functions[END_REF][START_REF] Wang | Multi-sensors data fusion algorithm of improved Dempster-Shafer evidence reasoning[END_REF][START_REF] Peri | Fundamentals of the Dempster-Shafer Theory[END_REF]. The evidential framework can explicitly represent each cell as the occupied, free, unknown and conflict state. There are many previous studies to introduce the Dempster-Shafer approach for modeling of evidential occupancy grid map (EOGM) [START_REF] Karaman | Comparison of different grid based techniques for real-time map building[END_REF][START_REF] Cao | Laser-scanner grid map building based on Dempster-Shafer evidence theory[END_REF]. Karaman, et al. compared the performance of various types of approaches(including Dempster-Shafer approach) for environment modeling, and Cao, et al. described the 2D grid map building process using Demspter-Shafer evidence theory. In this literature, simple laser beam model and Demspter-Shafer approach based information fusion for map building were provided. To create the EOGMs for the surrounding environment of the vehicle, SLAM (Simultaneous Localization And Mapping) is applied basically. Celmens, et al. applied the evidential theory to SLAM application. This paper evaluated the evidential theory based SLAM with collision-free path planning. Evidential FastSLAM based on Rao-Brackwellized particle filter was used to build the EOGM and evaluated by comparing with POGM in [START_REF] Reineking | Evidential FastSLAM for Grid Mapping[END_REF]. Credibilist SLAM simultaneously estimates the vehicle pose and builds the EOGM using LiDAR [START_REF] Trehard | Credibilist simultaneous Localization And Mapping with a LIDAR[END_REF][START_REF] Trehard | On line mapping and global positioning for autonomous driving in urban environment based on evidential SLAM[END_REF]. This paper applies GraphSLAM which is a state-of-the-art method to obtain the pose of a vehicle based on the graph optimization technique [START_REF] Grisetti | A Tutorial on Graph-Based SLAM[END_REF].

The intelligent vehicles equipped with perception sensors (such as LiDAR and Stereo camera), GPS receiver, and motion sensors are able to create the EOGM for their passing trajectory. If these vehicles have a connection to the Internet, the created EOGMs can be shared with other vehicles through a cloud map service. The cloud sharing of the multi-vehicle EOGMs provides many benefits to operate intelligent vehicles [START_REF] Guo | A Low-Cost Solution for Automatic Lane-Level Map Generation Using Conventional In-Car Sensors[END_REF][START_REF] Redzic | A Location Cloud for Highly Automated Driving[END_REF][START_REF] Meiling | The construction method of GIS for autonomous vehicles[END_REF][START_REF] Seif | Autonomous Driving in the iCity-HD Maps as a Key Challenge of the Automotive Industry[END_REF]. The intelligent vehicle can plan the global behaviors and motions by reflecting the distant area's EOGM downloaded from the cloud. The accuracy and reliability of EOGM in the cloud will be improved by merging the multiple EOGMs from several intelligent vehicles. Also, vehicles without perception sensors can use the EOGM generated by other vehicles if they only have an Internet connection.

For the cloud sharing of multi-vehicle EOGMs, there are two technical issues to be addressed. First, the EOGM must have a standard data format to be shared with other vehicles and to deal with the large-scale area. In the small area such as an indoor space, the EOGM is able to be implemented with single rectangular matrix data structure. However, the single rectangular matrix is not efficient format to cover the large-scale area for multiple vehicles. Second, evidential update method of EOGM is required to merge the new EOGMs updated from multiple on-road vehicles and old EOGM in the cloud. The evidential update also must take into account the aging factor of the old EOGM.

The main contribution of this paper is to propose a cloud update framework of evidential occupancy grid maps (EOGMs) for multiple intelligent vehicles in the large-scale real road environment. Each intelligent vehicle creates the EOGM for their moving trajectory based on GraphSLAM using LiDARs, motion sensors, and low-cost GPS receiver. The created EOGMs are uploaded and merged to EOGM in the map cloud. A standard tiling method of geodetic quad-tree tile system is applied to manage the EOGMs. The tiling method provides a common interface and data format to manage the large-scale EOGMs. Dempster combination rule of the evidential theory provides a theoretical basis to merge the new EOGM tiles to old EOGM tiles with consideration of the aging factor of old tiles. Experiments were performed to evaluate the cloud update framework of multi-vehicle EOGMs in the large-scale area. This paper is organized as follows. Section 2 presents overall cloud update framework of multiple vehicle EOGMs. Section 3 describes the EOGM mapping based on the GraphSLAM. Section 4 describes the cloud update of EOGM based on a geodetic quad-tree tile system. Section 5 provides experimental results of the proposed cloud update framework. The final section provides conclusion and future works.

Cloud Update Framework for Multi-Vehicle

Figure 1 shows the overall framework for the EOGM cloud update of multi-vehicle mapping. The static objects for the driving environment are mapped into EOGM. Since the EOGM is based on the evidential theory, it is robust to the disturbance of dynamic objects such as moving vehicles, pedestrians, and motorcycles. The large-scale EOGM is split into local tiles based on a geodetic quad-tree tile system. The geodetic quad-tree tile system provides a standard tiling method for the large-scale environment. The cloud update framework consists of two steps. The first step is multi-vehicle mapping. The multiple vehicles create EOGM tiles for their driving trajectory. Each vehicle must be equipped with perception sensor (such as LiDAR), motion sensor, and low-cost GPS receiver for the EOGM tiles generation. GraphSLAM is used to create the optimal EOGM tile of the passing trajectory when the vehicle is out of the tile boundary. For the initialization of GraphSLAM, the position from low-cost GPS receiver and EOGM tiles downloaded from the cloud is used. Section 3 describes the details of the GraphSLAM-based EOGM mapping. The second step is cloud update of the EOGM tiles from multi-vehicle. The multiple vehicles upload their created EOGM tiles when the generation is complete and the network to cloud is available. The uploaded EOGM tiles are merged into the old EOGM tiles in the cloud. Since the EOGM tiles follow the standard tiling method of the geodetic quad-tree tile system, spatial matching between the uploaded EOGM tiles and the old EOGM tiles in the cloud is straightforward. The evidential framework of EOGM provides a theoretical basis to merge the new EOGM tiles to old EOGM tiles with taking into account the aging factor of the old tiles. Section 4 describes the details of the cloud update of EOGM tiles. Based on the iteration of the two-step framework for the EOGM cloud update, the EOGMs in the cloud keeps up-to-date with the change of the static environment.

Mapping of Evidential Occupancy Grid Map

Evidential Occupancy Grid Map (EOGM)

An occupancy grid map (OGM) models the driving space using discrete grids or voxels. This paper uses two-dimensional (2D) grid map to model the driving space due to the constraint of network bandwidth and storage space, but the same process of 2D mapping can be extended to three-dimensional (3D) voxel map. The unit of the discrete space is a cell. The cell contains information about whether the corresponding space is occupied by static obstacles. Probabilistic approaches can be applied to infer the cell occupancy. The probability of a cell is zero when the cell is a free state, whereas the probability is one when the cell is an occupied state. However, the probability cannot deal with the situation of unknown (no sensor information for the cell due to blind spot) and conflict (different sensor data for the same cell due to moving objects). In the probabilistic approach, the cell probability of both situations is set to 0.5, which is ambiguous probability to represent the unknown and conflict.

To handle the unknown and conflict state, an evidential approach based on Dempster-Shafer (DS) is applied to build an evidential occupancy grid map (EOGM). The cell of OGM has to consider the two states of free (F) and occupied (O) as a frame of discernment Ω = {F, O}. In the EOGM based on the DS theory, the states of the cell can be extended to the power set 2 Ω = {F, O, Ω, ∅}, which is the set of all subsets of the Ω = {F, O}. This means that the EOGM can consider two more cell states of Ω and ∅ than the probabilistic occupancy grid map (POGM). For elements of the power set, a mass function m can be used to quantify the evidence of the hypothesis. The mass of m(F) and m(O) represent the evidence of the cell is free and occupied, respectively. The mass of m(Ω) that considers the union of F and O represents the cell is unknown, and m(∅) represents the cell is conflicted by a different source of information. Based on the definition of the mass function in the evidential framework, the sum of mass functions for the power set must be one.

EOGM Mapping

The four masses of each cell in the EOGM can be determined based on the sensing information of time-of-flight (TOF) range sensors, such as LiDAR, radar, and ultrasonic sensors. Before the mapping process, the sensing information reflected by ground is eliminated by pre-processing. This paper uses the LiDAR as a representative of the TOF sensor. The EOGM mapping process consists of two steps: local EOGM mapping using LiDAR sensor model and global EOGM mapping based on the accumulation of the local EOGMs.

Generation of Local EOGM Using LiDAR

The local EOGM represents an instant evidential grid map generated by a single scan of the LiDAR. Due to the scanning mechanism of the LiDAR, the local EOGM can be modeled as the polar grid map having the origin at the center of the LiDAR, as shown in Figure 2 [START_REF] Trehard | Credibilist simultaneous Localization And Mapping with a LIDAR[END_REF]. In the angular sector (a) of Figure 2, the red cells on the polar grid map indicate the occupied states, which can be inferred by the echo of LiDAR beam against potential obstacles. The array of green cells before the occupied cell in distance axis indicates the state of free, and the array of blue cells after the occupied cell represents the unknown state. For the multiple echoes in one angular sector (Figure 2b), all the echoes of LiDAR set the cell as occupied, but the only array of cells located before the occupied cell that has minimum distance is set by free, and the remaining cells are set to the unknown state. The state of cells for the sensor model can be used to determine the initial value of mass functions for the certain cell {r, θ}, denoted by [m L{r,θ} (F), m L{r,θ} (O), m L{r,θ} (Ω), m L{r,θ} (∅)], using

m L{r,θ} (O) = λ, if occluded 0, otherwise m L{r,θ} (F) = λ, if f ree 0, otherwise (1)
m L{r,θ} (Ω) = 1, if unknown 1 -λ, otherwise

Global EOGM Generation

Since the local EOGM is temporal information about one scan of LiDAR and contains noise from moving objects, the local EOGM cannot be used as an environment model to represent a drivable space and static obstacle. The local EOGM should be incrementally integrated to the global EOGM to stably provide the drivable space and static obstacle with removing the moving object disturbance, as shown in Figure 3. The global EOGM is represented by the grid plane of Cartesian coordinates.

In the beginning, all cells in the global EOGM are initialized by the vacuous mass function

m G{i,j} (F) = 0, m G{i,j} (O) = 0 (2) m G{i,j} (Ω) = 1, m G{i,j} (∅) = 0
that represents the no prior information. After the initialization, the local EOGMs of successive LiDAR scans are incrementally merged into the global EOGM. For the merging of a local EOGM to global EOGM, the local polar EOGM is converted to the local Cartesian EOGM which has the same grid configuration with global EOGM. The converted local EOGM in the Cartesian coordinates is then translated and rotated to align to the corresponding global EOGM by using the global pose (position and heading) from the GraphSLAM with the LiDAR extrinsic calibration parameters. The LiDAR extrinsic parameters are calibrated before the mapping process.

The mass function of each cell in the converted local EOGM can be represented by m L{i,j} that has the same index (i, j) with the global EOGM. Dempster combination rule of evidence theory is used to merge the local EOGM m L{i,j},t at time t with the global EOGM m G{i,j},t-1 at time t -1 :

m G{i,j},t = m G{i,j},t-1 ⊕ m L{i,j},t . (3)
The Dempster combination rule consists of the conjunctive combination rule (4) and the Dempster normalization (5) [START_REF] Moras | Moving Objects Detection by Conflict Analysis in Evidential Grids[END_REF].

∀A ⊆ Ω, m 1∩2 (A) = ∑ B∩C=A|B,C⊆Ω m 1 (B) • m 2 (C) (4) m 1⊕2 (A) = m 1∩2 (A) 1 -m 1∩2 (∅) , ∀A ⊆ Ω, A = ∅ (5) m 1⊕2 (∅) = 0
The m 1∩2 (A) means the result of conjunctive combination of m 1 and m 2 for state A, and m 1⊕2 (A) is the result from Demspter combination.

Wasting memory

GraphSLAM

To generate the global EOGM by incrementally integrating the local EOGMs, the pose (position and heading) of a vehicle must be known. The state-of-the-art method to obtain the pose of a vehicle is GraphSLAM based on the graph optimization technique [START_REF] Grisetti | A Tutorial on Graph-Based SLAM[END_REF]. The SLAM problem can be represented using a graph (Figure 4) consisted of nodes that represent random variables of stochastic process and edges which represent constraints between the nodes.

The nodes of x 1:t = [x 1 , • • • , x t ] represents a sequence of the poses for the discrete time steps 1 to t. The state of pose is composed of latitude, longitude, and heading. Low-cost GPS data is used to initialize and bound the state of poses. The map node m represents the static environment that is available to be recognized by perception sensors. The global EOGM is applied as the map node. The map node is initialized by downloaded EOGM from the cloud. The nodes of u The edges are composed of a state transition model (blue arc) and a measurement model (green arc), as shown in the Figure 4. The state transition model g(u t , x t-1 ) predicts the present pose x t from the previous pose x t-1 with the input u t . A constant velocity model is used for the state transition model with the motion inputs of the vehicle speed and yaw rate. We assumes the error u,t between the pose x t and the predicted pose g(u t , x t-1 ) follows the Gaussian distribution with covariance P u . Edge constraint of the state transition model is represented by a quadratic form of the transition error u,t and covariance P u , as described in [START_REF] Schütz | Laser-Based Hierarchical Grid Mapping for Detection and Tracking of Moving Objects[END_REF].

1:t = [u 1 , • • • , u t ]
x t = g(u t , x t-1 ) + u,t (6)
(x tg(u t , x t-1 )) T P -1 u (x tg(u t , x t-1 ))

The measurement model h(x t , m) estimates the local EOGM measurement z t based on the pose x t with the global EOGM m, as described in [START_REF] Smets | Analyzing the combination of conflicting belief functions[END_REF]. To construct the edge constraint of the measurement model, the error z,t between the measurement z t and measurement model h(x t , m) must be evaluated. However, since the z t and h(x t , m) are EOGM data type, the error z,t cannot be directly calculated with an analytical function. Instead of evaluating the error z,t , the similarity evaluation between two EOGMs based on the evidential scoring function is applied to construct the measurement model constraint. The evidential scoring function estimates the degree of matching of two EOGMs based on the sum of the Occupied masses in all cells of merged EOGM, as described in [START_REF] Wang | Multi-sensors data fusion algorithm of improved Dempster-Shafer evidence reasoning[END_REF]. The merged EOGM is obtained by applying the conjunctive combination rule (4) for the two EOGMs, m 1 = z t and m 2 = h(x t , m). Based on the scoring function, the best matching pose x t which has the maximum score is evaluated, as described in [START_REF] Peri | Fundamentals of the Dempster-Shafer Theory[END_REF]. The best matching pose x t is used to construct the edge constraint with the pose x t , as described in [START_REF] Karaman | Comparison of different grid based techniques for real-time map building[END_REF]. The error between the pose x t and the best matching pose x t is assumed to follows the Gaussian distribution with covariance P z .

z t = h(x t , m) + z,t (8)
scoring(m 1 = z t , m 2 = h(x t , m)) = ∑ ∀cells m 1∩2 (O) (9)
x t = argmax x scoring(m 1 , m 2 ) ( 10
)
(x t -x t ) T P -1 z (x t -x t ) (11)
A cost equation ( 12) is derived by sum of the log-likelihood constraints [START_REF] Schütz | Laser-Based Hierarchical Grid Mapping for Detection and Tracking of Moving Objects[END_REF] and [START_REF] Karaman | Comparison of different grid based techniques for real-time map building[END_REF]. When the cost function is minimized, the entire poses

x 1:t = [x 1 , • • • , x t ] are optimized.
The cost optimization problem is solved by g2o library, which is open source c++ framework to optimize the nonlinear least squares problems [START_REF] Kümmerle | 2 o: A general framework for graph optimization[END_REF].

J = ∑ t (x t -g(u t , x t-1 )) T P -1 u (x t -g(u t , x t-1 )) + ∑ t (x t -x t ) T P -1 z (x t -x t ) (12)

Cloud Update of EOGMs

We can create the global EOGM for small areas (less than one kilometer) using single rectangular matrix. However, the single rectangular matrix has limits to represent the global EOGM for the large-scale environment. First, the rectangular matrix of global EOGM causes wasting of memory space, as shown in Figure 3. Second, coordinate conversion between the geodetic coordinate and EOGM Cartesian coordinate for large-scale area causes a conversion approximation error. The global EOGM is based on a 2D Cartesian coordinate (east and north axis with meter unit), whereas the position for large areas on the Earth is normally represented by the geodetic coordinates (latitude and longitude). For the mapping of global EOGM in large areas, the mapping position in the geodetic coordinate such as WGS84 (World Geodetic System 1984) must transform into a Cartesian coordinate of global EOGM in order to be integrating the local EOGM. The coordinate conversion from geodetic coordinate to Cartesian coordinate is performed based on the tangential plane approximation for one reference point, as shown in Figure 5a. However, when the distance between the reference point and vehicle position is larger than 10 km, the approximation error of the coordinate conversion causes a distortion of global EOGM mapping [START_REF] Jo | Fast GPS-DR Sensor Fusion Framework: Removing the Geodetic Coordinate Conversion Process[END_REF].

Reference point East

North

Reference point East

North

Approximation error

Approximation error (a) (b) To overcome the problems, a previous study presented a method to manage the grid map by tiles that split the space into many small areas [START_REF] Alsayed | PML-SLAM: A solution for localization in large-scale urban environments[END_REF]. However, since this tiling method can have a different configuration depending on the implementation of grid map, it is difficult to share or reuse the grid map by other vehicles. For the standardization of data structure of the global EOGM, this paper applies a geodetic quad-tree tile system to the global EOGM management. The geodetic quad-tree tile system is widely used for quick access and management of map in the digital map industry, such as Bing, HERE and Google map.

Geodetic Quad-Tree Tile System

The quad-tree tile system provides a common tile structure that can be used all around the world. The world can be split into hierarchical levels of tiles, which an upper-level tile is divided into four lower level tiles, as shown in Figure 6. The entire world is level 0, and the tiles of level 1 are divided by four regions of the level 0 tile, as shown in Figure 6a. The top two tiles of the level 1 are empty to keep the angular range of geodetic coordinate system (latitude: ±90 • , longitude: ±180 • ) and to be the same angular size for each tile as [START_REF] Reineking | Evidential FastSLAM for Grid Mapping[END_REF].

Degree o f width and height per tile = 360 The quad-key is a unique identification number for each tile. This number is assigned sequentially from the bottom left tile as shown in Figure 6a. The tiles of level two are also split into four tiles from one tile of level one tile. Figure 6b shows the level two tiles divided from the level one tile that contains a star mark (position of Eiffel tower). The quad-key of the level two tiles starts with the quad-key of the parent tile, and then the identification numbers are added, as shown in Figure 6b. By repeating this process, we can reach the level of tiles of the reasonable size for the global EOGM. Figure 6d shows tiles of level 16, which has latitude and longitude sizes of about 0.0055 • (width: 403.1 m, height: 610.88 m) with a quad-key 1220002130322231.

Geodetic Quad-Tree Tiling of Global EOGM

By assigning global EOGMs to the tiles of a certain level, we can manage the group of global EOGMs that can cover all around the world. To assign a global EOGM to a tile, the geodetic coordinate of the tile should be converted into Cartesian coordinate. For the coordinate conversion, the lower left corner point of the tile is set to be an origin of Cartesian coordinate that has east and north axis. Using the coordinate conversion based on the origin point, the tile size and certain position in the geodetic coordinate can be converted into the size and position (east and north) on the Cartesian coordinate. Figure 6e There are many advantages when we use the quad-tree tile as the management system of the global EOGMs. First, the quad-tree tile representation can reduce the plane approximation error of coordinate conversion because the conversion is performed based on each origin point of the tiles, as shown in Figure 5b. Second, it is very effective in terms of memory utilization because only the global EOGM in the area of interest with the certain tile size need to be created and loaded into the memory. Figure 7 shows an example of the geodetic quad-tree tile-based global EOGM, which has same local mapping configuration with the Figure 3. Compared to Figure 3 with one large rectangular global EOGM, the quad-tree tile based management is more efficient for the memory utilization because the global EOGM for areas of no interest do not need to be updated. Third, since all quad-tree tiles have a unique quad-key, the searching and data management of the global EOGM for cloud update is straightforward. Finally, the most important thing is the quad-tree tile system provides the common interface and data format for the global EOGM, which is compatible for all around the world. Based on the common interface and format, multiple vehicles around the world can share the global EOGM for the environment perception.

New tile creation

Cloud Update of EOGM Tiles

The created EOGM tiles from multiple intelligent vehicles are uploaded when the GraphSLAM and quad-tree tiling are complete and the network to cloud service is available. The uploaded EOGM tiles are matched to the exiting EOGM tiles in the cloud using the quad-key matching. Then, the uploaded EOGM tiles are merged into the matched EOGM tiles in the cloud. The merged EOGM tiles in the cloud can be downloaded to other intelligent vehicles for the environment modeling and SLAM initialization. There are many benefits when the global EOGM tiles generated by individual vehicles are shared with multiple vehicles via cloud service. First, the intelligent vehicle can obtain the environmental information about the areas of future visits by downloading the global EOGMs of the area from cloud. Second, the GraphSLAM performance of the individual vehicles for the new visited area is improved by initializing the SLAM using the downloaded global EOGMs. Finally, it is possible to obtain the accurate and reliable global EOGMs by integrating the several EOGM tiles from multiple vehicles.

To merge the uploaded EOGM tiles from multiple intelligent vehicles to the existing EOGM tiles in the cloud, Dempster combination rule of evidential theory with aging effect is applied. The aging of the old EOGM tiles in the cloud must be considered for the merging because there is a probability of changing the static environment as time passes. The aging effect can be implemented by decay factor α that can be described by

α = exp( t old -t new τ ) (14)
where t new is the time of grid map uploaded from the vehicles, t old is the time of grid map stored in the cloud and τ is the time constant that determines the decaying rate. By applying the decay factor α to the EOGM, we can obtain the aging EOGM as

α m G{i,j},t new (A) = α • m G{i,j},t old (A), A ⊂ Ω (15)
α m G{i,j},t new (Ω) = 1α + α • m G{i,j},t old (Ω)

Then, the global EOGMs in the cloud are updated using the Dempster combination rule:

m G{i,j},t new = α m G{i,j},t new ⊕ m L{i,j},t new . ( 16)

Experiments

Experimental Environments

A test vehicle is equipped with six automotive LiDARs (Valeo ScaLa) and a single GPS receiver (NovAtel receiver without any corrections). The six four-layers scanning LiDARs can cover the 360 • Field-of-view (FoV) and, the maximum detection range of LiDAR is 150 m, the horizontal resolution is 0.25 • , the vertical resolution is 0.8 • , and the distance resolution is 0.1 m. The Cartesian local EOGM covers a distance of 100 m around the vehicle with 0.2 m grid size. The confidence coefficient λ which is represented in (1) is set to 0.7 with considering of the LiDAR noise characteristics. The test vehicle's CAN (controller area network) provides the vehicle motion data from the onboard sensors (wheel speed sensors, yaw rate, steering angle sensors, and accelerometer). The GPS data used to initialize and bound the pose nodes of GraphSLAM. An embedded computer (Intel i7-4700EQ CPU 2.4 G) on the test vehicle runs GraphSLAM to create the EOGM for the recording trajectory. The generated EOGM are converted to geodetic quad-tree tile format and uploaded to the map cloud server (Amazon Web Services) via a 4G mobile network.

Tiled EOGM Generation for Large-Scale Area

An experiment was performed to evaluate the large-scale coverage of the quad-tree tiled EOGMs. The test trajectory was about 20 km (Figure 8). The red tile shows a single large tile that covers the entire trajectory, and the blue tiles represent the tiles of the quad-tree level 17. The length of one side of the level 17 tiles is about 250 m. The reason why the blue tiles exceed the red tile at the boundary is that quad-tree tiles have their own global position and size. Figure 9 shows the memory usage of various level tiles for the mapping. The smaller size (higher level) tiles improve the memory efficiency because they effectively place tiles only on trajectories without wasting space. In addition, this memory efficiency is achieved by little computational cost because only calculation of quad-key is required for tile based EOGM mapping.

Since the EOGM tiles have the same data structure as images, it is available to store the global EOGMs into image files. The amount of storage required for the entire level 17 tiles is about 654 megabytes using the bitmap format, but it can be reduced to an average of approximately 1.28 megabytes using the lossless data compression of PNG (portable network graphics). The one file of the compressed EOGM tile is about 150 kilobytes, and this size is acceptable to share with other vehicles via mobile Internet.

In addition to the memory efficiency, the coordinate conversion error from the plane approximation can be reduced by the adaptive reference point of each quad-tree tile. Since the approximation error of the coordinate conversion is difficult to evaluate using quantities value due to the dependence of the evaluation location, we evaluate the approximation performance by analyzing the distance between the reference point (origin) and conversion point. When there is one reference point (start point) for the entire trajectory, the maximum distance from the reference to the vehicle position is about 15 km. However, for the level 17 quad-tree tiles, the maximum distance is reduced to about 350 m. This means that the the quad-tree tile system reduces the approximation error through the mechanism as shown in Figure 5.

Cloud Update of EOGMs

Another experiment was performed to evaluate the EOGM GraphSLAM and the cloud update of EOGM tiles for the multiple mappings. The test site was an urban area where has moving objects and parked cars, as show in Figure 10. Although we had only one test vehicle, the cloud update experiments for the multi-vehicle was available by multiple driving, mapping, and uploading. We performed six tests at different times to create and update the multiple EOGMs, as shown in the mapping trajectories of Figure 10.

The test site can be evenly divided into several tiles of the quad-tree level 17. The tile boundary and tile ID for level 17 tiling of the test site is shown in Figure 10. Each created EOGM tiles for each mapping are uploaded to the cloud and merged into one EOGM tiles, as shown in Figure 11. For the EOGM cloud update, new EOGM tile must be matched with existing tiles in cloud. Since each tile has a unique quad-key, it is straightforward to search and manage the tiles matching using a hash table. Figure 12 shows the EOGM update over time for the red box area in the Figure 10. The EOGM for the area were updated three times: EOGM mapping 1 (09:12), mapping 2 (09:34), and mapping 5 (15:46). Figure 12a shows the aerial image. Figure 12b shows the EOGM after first cloud update of mapping 1. Although there are many moving vehicles and pedestrians, the only static objects such as parked cars and non-moving cars are mapped. Figure 12c shows the EOGM of mapping 2, and it merged to (b). The merged result of Figure 12b,c is shown in Figure 12d. The static evidence (red in the white box) of a changed objects are blurred due to EOGM merging with taken into account the aging. The time constant τ which determines the decaying rate is 24 h. This value means that the environment which does not change in 24 h can be considered as static environment. After about six hours, EOGM was generated by mapping 5 as shown in Figure 12e. Some parked cars have left and new parked car has appeared. The merged result of Figure 12d,e is shown in Figure 12f. Since six hours have passed, old maps are reflected weakly and new updated maps are more reflected.

Conclusions

This paper presents a cloud update framework of multi-vehicle EOGMs based on the GraphSLAM, geodetic quad-tree tile system, and the evidential cloud update of EOGM tiles.

(1) The evidential theory for the EOGM provides a basis to manage all types of cell states such as free, occupied, unknown, and conflict, which cannot be managed by probabilistic theory. In addition, the Dempster combination rule allows merging the uploaded EOGM tiles from the multiple vehicles to the old EOGM tiles in the cloud.

(2) The geodetic quad-tree tile system provides a common tile interface to manage the EOGMs for the entire world. The sharing of the global EOGMs of multiple intelligent vehicles is based on the geodetic quad-tree tile format. In addition, the management based on the quad-tree tile system provides many benefits, including the improvement of memory utilization, reduction of a coordinate conversion error, and acceleration of searching speed by tree based searching. Since the geodetic tiling method only requires the calculation of quad-key, the computing cost can be ignored.

(3) The proposed algorithm was evaluated through two experiments. The first experiment evaluated the ability of the EOGM tiling to manage large-scale area (more than 20 km). The second experiment verified the cloud update process of multiple EOGM tiles in consideration of aging effect. After uploading the EOGMs from individual vehicles, a post processing and quality evaluation of the uploaded EOGMs are required in cloud. Since the computing time of post processing and evaluation are not fast enough, the real-time sharing is not available currently. For our future work, we have a plan to share the EOGMs in real-time by reducing the computational time of post processing and evaluation. If the post processing technique is mature in the future, real-time download will be possible.

The presented study limits the environment model into 2D plane coordinates because of the computation and memory constraints. Future studies will extend the 2D EOGM to 3D EOVM (evidential occupancy voxel map) with a GPU parallel processing in order to provide a more detail environment model for intelligent vehicles.

Figure 1 .

 1 Figure 1. Two-step framework for EOGM cloud update of multi-vehicle mapping.

mFigure 2 .

 2 Figure 2. Local EOGM generated from the LiDAR (F: green, O: red, Ω: blue). (a) Local EOGM in polar coordinate, (b) Unit polar grid.

Figure 3 .

 3 Figure 3. Mapping of Global EOGM based on the merging of local EOGM.

  represent a sequence of the vehicle motion (such as the yaw rate and speed) for the corresponding times of the states. The nodes of z 1:t = [z 1 , • • • , z t ] represent the local EOGM measurement generated by the local EOGM mapping of LiDAR point cloud in the Section 3.2.1.

Figure 4 .

 4 Figure 4. SLAM represented by a directed graph of a dynamic Bayesian network (DBN) based on the assumptions of the Markov process.

Figure 5 .

 5 Figure 5. Coordinate conversion between geodetic coordinate and Cartesian coordinate for mapping of global EOGM. Plane approximation (a) without the quad-tree tile and (b) with the quad-tree tile.

Figure 6 .

 6 Figure 6. Geodetic quad-tree tile system for global EOGM management. tiles of level 1-16, and tile for global EOGM in Cartesian coordinate.

  shows an example of coordinate conversion for a level 16 tile contained Eiffel tower. The lower left corner of the tile is set to origin of the Cartesian coordinate. The size (360 • /2 16 ) of the tile in geodetic coordinate is converted into size (403.1 m, 610.88 m) of Cartesian coordinate. The position of Eiffel tower (48.8582 • , 2.2947 • ) in the geodetic coordinate is converted to the position (297.24 m, 222.28 m) in the Cartesian coordinate of level 16 tile. In the same way, we can create and access the global EOGMs for the specific level tiles in the worlds.

Figure 7 .

 7 Figure 7. Mapping of Global EOGM based on the quadtree tile.

Figure 8 .

 8 Figure 8. The test trajectory (20 km) to evaluate the large-scale coverage of the geodetic quad-tree EOGM tiles.

Figure 9 .

 9 Figure 9. Memory usage for the various level tiles to evaluate the memory efficiency of the quad-tree tile system.

Figure 10 .

 10 Figure 10. Trajectories and quad-tree tiles for the six test driving at different times in an urban area.

1Figure 11 .

 11 Figure 11. Cloud update of the multiple EOGM tiles for the different areas at different times. EOGM tiles of six test mapping are uploaded and merged into EOGM tiles in the cloud.

tiles Level 2 tiles Level 3 tiles

	2	3	12	13	122	123
	0	1	10	11	120	121
					North [m]
				610.88
				Conversion to
				Cartesian
				coordinate
				222.28
							East [m]
					297.24	403.10

• /2 tile level (13) Level 1 (48.8582, 2.2947) (48.8582, 2.2947) (48.8582, 2.2947)

(297.24, 222.28)

	Lat: 270			Lat: 90	Split into	Lat: 90
					four tiles
	Lat: -90 Lon: -180	Lon: 180	Split into four tiles	Lat: -90 Lon: 0	Lon: 180	Lat: 0 Lon: 0	Lon: 90
					(48.8582,
					2.2947)

Level 16 tiles (48.8582, 2.2947) 1220 0021 3032 2232 1220 0021 3032 2233 1220 0021 3032 2230 1220 0021 3032 2231

	Lat: 48.8672		Lat: 48.8617
	Lat: 48.8562	Zoom in	Lat: 48.8562
	Lon: 2.2852 Lon: 2.2961	Lon: 2.2906 Lon: 2.2961

Level 16 tile (48.8582, 2.2947)

Funding Acquisition, M.S. Funding: This work was financially supported by the BK21 plus program (22A20130000045) under the Ministry of Education, Republic of Korea, the Industrial Strategy Technology Development Program (No. 10039673, 10060068, 10079961), the International Collaborative Research and Development Program (N0001992) under the Ministry of Trade, Industry and Energy (MOTIE Korea), and National Research Foundation of Korea (NRF) grant funded by the Korean government (MEST) (No. 2011-0017495).


Tengfei Chang
email: tengfei.chang@inria.fr

Thomas Watteyne
email: thomas.watteyne@inria.fr

Xavi Vilajosana
email: xvilajosana@uoc.edu

Competition: OpenWSN, a Development Environment for 6TiSCH

6TiSCH is a standardization effort at the IETF that combines industrial performance with seamless Internet connectivity. It brings together a number of existing standards, and adds scheduling capabilities. One of these standards is IEEE802.15.4 Time Slotted Channel Hopping (TSCH), a medium access control layer technique that provides high reliability and ultra low power consumption. Another is RPL, a routing protocol that turns the topology into a multi-hop mesh network. OpenWSN is an open-source implementation of 6TiSCH. It provides a set of developing tools and supports multiple hardware platforms. OpenWSN is the basis of our solution for the EWSN 2019 dependability competition.

1 6TiSCH, the Industrial IoT Protocol Stack 6TiSCH [START_REF] Dujovne | 6TiSCH: Deterministic IP-enabled Industrial Internet (of Things)[END_REF] defines a protocol stack for Industrial IoT, depicted in Fig. 1. This protocol stack is rooted in the IEEE802.15.4 physical layer. Time Slotted channel hopping (TSCH), as Medium Access Control (MAC) layer protocol, provides ultra-high reliability and ultra-low energy consumption. The 6TiSCH protocol stack uses the 6top protocol [START_REF] Wang | 6TiSCH Operation Sublayer (6top) Protocol (6P)[END_REF], one or more 6TiSCH scheduling functions, and the 6LoWPAN adaption layer to transport large IPv6 packets into short IEEE802.15.4 frames [START_REF] Thubert | Compression Format for IPv6 Datagrams over IEEE[END_REF]. RPL [START_REF] Winter | RPL: IPv6 Routing Protocol for Low-Power and Lossy Networks[END_REF] is a routing protocol for low-power lossy networks. It allows communication both from devices inside the mesh network to its gateway ("upstream") and from the gateway to devices ("downstream"). 6TiSCH uses CoAP [START_REF] Shelby | The Constrained Application Protocol (CoAP)[END_REF] at its application layer, including for carrying the secure joining messages.

According to the competition categories, there are two traffic patterns, multipoint-to-point (MP2P) and point-tomultipoint (P2MP). RPL supports both traffic pattern. Moreover, the competition evaluation metrics are end-to-end reliability, energy consumption and latency. 6TiSCH is a natural fit for these types of traffic. Adapting the scheduling function further allows us to trade-off lower latency for higher energy consumption.

The following two sections describes TSCH and RPL, the main two protocols that influence network performance.

Time Slotted Channel Hopping (TSCH)

TSCH achieves high reliability through channel hopping, and ultra low-power consumption through time synchronization.

In a TSCH network, time is split into time slots. A pair of communicating neighbor nodes are synchronized when their time slot boundaries are aligned. Link-layer frames are exchanged within slots. Since nodes are synchronized, the receiver knows when in the slot the transmitter starts sending a frame. This time offset is TxOffset in Fig. 2. This allows a receiving node to only turns on its radio when it is about to receive a frame (withing a guard time of TxOffset to account for clock drift). This mechanism significantly the radio-on time, hence the average energy consumption of the device. All nodes share a common sense of time, as each time slot is identified by the Absolute Slot Number (ASN), an ever-increasing number. This shared notion of time allows neighbor nodes to turn on their radio only when the are scheduled to exchange a frame.

In TSCH, the 83 MHz frequency band at 2.4 GHz is split into 16 equally-spaced channels. All communication is orchestrated by a schedule. That schedule can be represented as a slot × channel matrix, as shown in Fig. 2. When neighbor nodes are about to exchange a frame, they use Eq. 1 to compute the frequency they will tune their radio to. chO f f set is the logic channel offset used by a slot. numChannels indicates the number of channels used by the network. Fig. 2 shows a 5 node network using 3 channels to communicate with one another. Channel hopping is the mechanism by which successive frames exchanged between two neighbors are done so on a different frequency. Channel hopping is effective at combating external interference and multiple-path fading [START_REF] Watteyne | Reliability Through Frequency Diversity: Why Channel Hopping Makes Sense[END_REF].

channel = (ASN + chO f f set)%numChannels (1)
3 Routing Protocol for Low-Power and Lossy Networks (RPL)

The RPL routing protocol enables both upstream and downstream communication. It organizes the network into a Directed Acyclic Graph (DAG) by assigning a rank to each node in the network. The rank of a node indicates the cost to send a packet to the root node. A node's rank increases as it is multiple hops away from the root, and because the it has to its routing parent is weak. Upstream communication a achieved when each node forward a message to its parent. Fig. 2 shows a DODAG with 5 nodes.

A node can have multiple neighbors in its parent set. Fig. 2 shows that node D has two neighbors which are closer than itself to the root: nodes B and C. In case of node B is powered off, node D can still forward its packets to node C, so no data is lost. Having multiple nodes in the parent set increases the reliability and stability of the network.

OpenWSN

OpenWSN [START_REF] Watteyne | OpenWSN: a standards-based lowpower wireless development environment[END_REF] is the reference open-source implementation of the 6TiSCH protocol stack. It mainly contains two parts: the firmware running on the embedded devices, and the software running on the PC. The firmware can be run on 11 different hardware platforms, including the TelosB board used in the competition. The software -called OpenVisualizer -contains a set of supporting scripts, including visualizing the status of the node in the network, or configuring one of the nodes as the root of the network.

The following sections introduces two features: HDLCbased serial debugging, and the build system.

HDLC-based Serial Debugging

In OpenWSN firmware provides a component called openserial. The component is responsible for forwarding the status of nodes, the packet to the root or Internet, and receiving commands from OpenVisualizer. We use High-Level Data Link Control (HDLC) as the framing mechanism on top of the raw serial communication. The status of the node includes the node's identifiers, the node's rank, the current ASN, whether the node is synchronized, as well as details about the node's schedule, neighbor and transmit queue. The OpenVisualizer gathers this status information printed by the nodes, and displays it conveniently on a dynamic web page.

A Complete Build System

OpenWSN uses SCons to build the 6TiSCH protocol stack for different target hardware platforms. Scons is a build system, similar to "Make", but written in Python. With the gcc toolchain for msp430, the source code can be compiled for the TelosB. OpenWSN also provides several bootloader scripts, allowing a user to program her devices after compiling the code and generating the image.

Conclusion

We propose to use OpenWSN for the EWSN 2019 dependability competition. This extended abstract introduces the 6TiSCH protocol, as well as its two main protocols: TSCH and RPL. We describe OpenWSN, and detail its HDLC-based serial debugging capability, and the build system used.

Figure 1 .

 1 Figure 1. In the 6TiSCH protocol stack, TSCH and RPL are the main two protocols that influence network performance.

Figure 2 .

 2 Figure 2. The TSCH schedule indicates to each now what to do in each time slot. The cell shaded blue is use by node D to send a frame to node C, and get a link-layer acknowledgement.