Hayani Meriem

Ali Mechkouri

George Chaaba

Philippe Moraru

Veron

M Hayani Mechkouri

A Chaaba

G Moraru

P Veron

PREDICTION OF MULTI-DIMENSIONAL MILLING BEHAVIOR

Keywords: multi-dimensional, milling process, behavior, Chatter Vibrations, RDEs, Matlab software. I

come

INTRODUCTION

Milling is a manufacturing process where material is removed by chip which is affected by vibrations. Chatter is one of the three vibration types [1], [START_REF] Tobias | Vibraciones en máquinas-herramientas[END_REF] and consists on remanufacturing of the manufactured surface. Firstly, the importance of chatter in vibrations behavior has been demonstrated experimentally by Nicolson [START_REF] Nicolson | Experiments with a Lathe Tool Dynamometer[END_REF] then Taylor warned industrials about its lower effect on machined piece quality in 1907 [START_REF] Taylor | On the art of cutting metals[END_REF]. Effectively, based on mechanistic method, Tlusty proved the direct impact of chatter on chip thickness [START_REF] Tlusty | Dynamics of Cutting Forces in End Milling[END_REF]. Afterwards, he developed a basic nonlinear system of cutting force in machining chatter [START_REF] Tlusty | Basic Non-Linearity in Machining Chatter[END_REF] which underlies the strong relationship between cutting force and chip thickness [START_REF] Martelloti | An analysis of the milling process[END_REF]. Whereas many works have been developed to avoid chatter [START_REF] Yao | On-line chatter detection and identification based on wavelet and support vector machine[END_REF] or to suppress it [START_REF] Al-Regib | Programming spindle speed variation for machine tool chatter suppression[END_REF]- [START_REF] Altintas | In-process detection and suppression of chatter in milling[END_REF] by varying tool spindle speed, for example, others have been conducted to predict it in time [START_REF] Zhongqun | Solution and Analysis of Chatter Stability for End Milling in the Timedomain[END_REF] and frequency [START_REF] Otto | Extension of Tlusty׳s law for the identification of chatter stability lobes in multi-dimensional cutting processes[END_REF] domains. Since computing stability via simulation milling process in time domain is time consuming, researchers are putting greater focus on predicting stability in frequency domain. The most popular criterion used to achieve this target is the tau-decomposition criterion which stability lobes diagram is based on [START_REF] Tlusty | Stability Lobes in Milling[END_REF], [START_REF] Altintaş | Analytical prediction of stability lobes in milling[END_REF]. Stability lobes diagram is still used to predict chatter in manufacturing by drilling [START_REF] Mehrabadi | Investigating chatter vibration in deep drilling, including process damping and the gyroscopic effect[END_REF], turning [START_REF] Altintas | Chatter Stability of Metal Cutting and Grinding[END_REF] and milling [START_REF] Tlusty | Dynamics of High-Speed Milling[END_REF], [START_REF] Ahmadi | Stability lobes in milling including process damping and utilizing Multi-Frequency and Semi-Discretization Methods[END_REF] systems thanks to its ability to give general overview idea about process stability in constant machining conditions by considering only one degree of freedom at a time. Nevertheless, general milling is a multi-dimensional cutting process. Hence the importance of substituting the Tau-decomposition by a suitable criterion. Therefore, this works aims to analyze stability in the Lyapunov sense [START_REF] Stépán | Retarded dynamical systems: stability and characteristic functions[END_REF] of nonlinear differential equations describing chatter dynamical system in multi-dimensional milling process. The second section aims to demonstrate the crossing from milling kinematics to cutting forces. So, milling tool is discretized into elemental cutting edges where differential cutting forces are modeled, computed and summed. The third consists on representing a multi-dimensional chatter milling process by a system of three retarded differential equations (RDEs). After that, we prove that the characteristic function associated is a quasi-polynomial whose degree is superior to one. Then, looking back on theorems defining stability of RDEs on the basis of the characteristic quasipolynomial root finder, some criteria are defined to determine the stability type (stable, asymptotically stable, unstable) of milling process. Hence, milling system behavior is determined in frequency domain. Afterwards, a case study of chatter milling process is borrowed from bibliography to apply the developed method of behavior prediction and Matlab software is used in computing phase. Finally, the outputs of this method are compared to milling simulation outcomes computed in time domain.

II.

MILLING KINEMATICS Machining process is assumed to be a series of workpiece/tool interactions where material is removed as chip. Particularly, in milling, tool is rotating around his axis with angular velocity Ω and moving, simultaneously, at feed rate𝒇 𝒛 to remove chip thanks to its 𝑧 teeth indexed by 𝐿 discretized into elemental cutting edges. Indeed, milling tool is discretized into 𝑛𝑘 disks of defined thickness 𝑑𝑏 𝑘 referenced by 𝑘 where each disk is containing 𝑧 elemental cutting edges indexed by (𝑘, 𝐿).

Fig. 1 : Discretization of milling tool

In the middle plane of the 𝑘 𝑡ℎ elemental disk, we define a 𝑘 𝑡ℎ med-plane perpendicular to the tool axis where each elemental cutting edge of the (𝑘, 𝐿) 𝑡ℎ edges (1 ≤ 𝐿 ≤ 𝑧) is located geometrically by an elemental angle engagement 𝜃 𝑘,𝐿 0 measured from 𝑌 ⃗ axis of a reference frame tied to the spindle(𝑂, 𝑋 , 𝑌 ⃗ , 𝑍 ) and 𝜅 𝑘 . 𝜅 𝑘 represents the axial immersion of the k th disk med-plane. Milling kinematics induce to material removal and a new surface forming in each workpiece/tool interaction called local swept surface (Error! Reference source not found.). To shift from an interaction to the next, milling tool is rotating with an angle Ω * 𝑡 where t is the time factor. Then, at the time 𝑡, the instantaneous angular position 𝜃 𝑘,𝐿 of the (𝑘, 𝐿) 𝑡ℎ elemental cutting edge is deducted as:

𝜃 𝑘,𝐿 = 𝜃 𝑘,𝐿 0 + Ω * 𝑡 (1) In the angular position 𝜃 𝑘,𝐿 (𝑡), the (𝑘, 𝐿) 𝑡ℎ elemental cutting edge, whose length is noted 𝑑𝑆 𝑘 , is applying a differential force 𝑑𝐹𝑡𝑎𝑟 ̅̅̅̅̅̅̅̅ 𝑘,𝐿 to remove the infinitesimal chip whose thickness and width are, respectively ℎ 𝑘,𝐿 and 𝑑𝑏 𝑘 . The (𝑘, 𝐿) 𝑡ℎ differential force 𝑑𝐹𝑡𝑎𝑟 ̅̅̅̅̅̅̅̅ 𝑘,𝐿 is applied on cutting edge middle and expressed as (2) based on linear model of milling force components [START_REF] Kaymakci | Unified cutting force model for turning, boring, drilling and milling operations[END_REF] in a reference system attached to the elemental tooth.

𝑑𝐹𝑡𝑎𝑟 ̅̅̅̅̅̅̅̅

Fig. 2 : Elemental milling force components

Computation of linear model milling force depends not only on elemental cutter and chip geometry but also on tangential, radial and axial ploughing and shearing coefficients [21] [22] which are representing force applied respectively per unit of area or length.

III.

CHATTER MODELLING In milling, discontinuity of edges engagement on material creates so-called cutting delay r. This phenomenon is called chatter and parameter r is defining duration between two successive teeth passage. In case of symmetric milling tool, this factor depends on teeth passing frequency N and number z.

𝑟 = 60 𝑧 * N

Fig. 3: Chip removal in chatter milling

Figure 1 shows the impact of chatter on instantaneous chip thickness. Indeed, total chip thickness is the sum of dynamic thickness h 𝑑 induced by chatter and quasi-static thickness h 𝑠𝑡 removed by rigid removal system [START_REF] Eksioglu | Discrete-time prediction of chatter stability, cutting forces, and surface location errors in flexible milling systems[END_REF] (3).

ℎ = ℎ 𝑠𝑡 + h 𝑑 (3)
In the context of a discretized milling tool into elemental cutting edges, we note ℎ 𝑘,𝐿 the instantaneous elemental chip thickness [START_REF] Eksioglu | Discrete-time prediction of chatter stability, cutting forces, and surface location errors in flexible milling systems[END_REF]express thicknesses quasi-static and dynamic elementary, respectively ℎ 𝑘,𝐿 𝑠𝑡 and ℎ 𝑘,𝐿 𝑑 as suite: Elemental dynamic chip thickness ℎ 𝑘,𝐿 𝑑 (𝑡) is due to the delay in passage of the (k, L) th edge in matter and elemental cutting edge position at time t defined by 𝜅 𝑘 and 𝜃 𝑘,𝐿 (𝑡). Resulting chatter phenomenon has an impact on the three-dimensional geometry of the instant chip expressed by the instant terms(𝛥𝑋(𝑡), 𝛥𝑌(𝑡), 𝛥𝑍(𝑡)) :

∆ ̅ (𝑡) = { 𝛥𝑋(𝑡) = 𝑋(𝑡) -𝑋(𝑡 -𝑟) 𝛥𝑌(𝑡) = 𝑌(𝑡) -𝑌(𝑡 -𝑟) 𝛥𝑍(𝑡) = 𝑍(𝑡) -𝑍(𝑡 -𝑟)
Chatter is a dynamic phenomenon represented by distance between two vibratory displacements, respectively 𝑋(𝑡), 𝑌(𝑡)and 𝑍(𝑡), spaced out in time by r delay in the direction of respective reference axes X, Y and Z of spindle tool repair (𝛥𝑋(𝑡), 𝛥𝑌(𝑡), 𝛥𝑍(𝑡)). We must remind that this section aim is to determine instantaneous global milling force 𝐹 𝑋𝑌𝑍 ̅̅̅̅̅̅ (𝑡) expression. For this purpose, differential milling forces are sorted, so that associated to non-engaged elementary cutting edges (zero) are eliminated thanks to the window function 𝑔(𝜃 𝑘,𝐿 ) defined as:

𝑔(𝜃 𝑘,𝐿 ) = { 1 𝑖𝑓 (𝑘, 𝐿) 𝑡ℎ 𝑡𝑜𝑜𝑡ℎ 𝑖𝑠 𝑒𝑛𝑔𝑎𝑔𝑒𝑑 0 𝑖𝑓 𝑛𝑜𝑡
Then, all non-zero differential milling forces are projected on a unique reference frame (milling tool reference frame, for example) thanks to 𝑇 ̿ 𝑘,𝐿 (𝑡) in order to be summed finally.

𝐹 𝑋𝑌𝑍 ̅̅̅̅̅̅ (𝑡) = ∑ ∑ 𝑑𝐹 𝑋𝑌𝑍 ̅̅̅̅̅̅̅̅ 𝑘,𝐿 (𝑡) 𝑧 𝐿=1 𝑛𝑘 𝑘=1 = ∑ ∑ 𝑇 ̿ 𝑘,𝐿 (𝑡) * 𝑔(𝜃 𝑘,𝐿 ) * 𝑑𝐹 𝑡𝑎𝑟 ̅̅̅̅̅̅̅̅̅ 𝑘,𝐿 (𝑡) 𝑧 𝐿=1 𝑛𝑘 𝑘=1 (4)
Where:

𝑇 ̿ 𝑘,𝐿 (𝑡) = { sin(𝜅 𝑘 ) *
Finally, the sum obtained by formulas (2)-( 5) is expressed as in [START_REF] Tlusty | Basic Non-Linearity in Machining Chatter[END_REF]: )

𝐹 𝑋𝑌𝑍 ̅̅̅̅̅̅ (𝑡) = {
Instantaneous milling force F XYZ ̅̅̅̅̅̅ (t) includes instantaneous cutting force a ̅(t) due to rigid body motion of the cutting tool and inducing to quasi-static chip thickness removal, besides of instantaneous chatter force α ̿(t) * ∆ ̅ (t) which depends on chatter vibration displacements and delay, implicitly. Chatter is a vibration type which is due to system flexibility. Then, instantaneous chatter milling force α ̿(t) * ∆ ̅ (t) is time response of dynamics milling system. So, on the basis of Newton's second law applied to the flexible milling system, a system of retarded differential equations is developed [START_REF] Li | Cutting force prediction and analytical solution of regenerative chatter stability for helical milling operation[END_REF]:

𝑀 ̿ [ 𝑋 ̈(𝑡) 𝑌 ̈(𝑡) 𝑍 (𝑡) ] + 𝐶 ̿ [ 𝑋 (𝑡) 𝑌 ̇(𝑡) 𝑍 (𝑡) ] + 𝐾 ̿ [ 𝑋(𝑡) 𝑌(𝑡) 𝑍(𝑡) ] = 𝛼 ̿(𝑡) * ∆ ̅ (𝑡) (7)
Where:

𝑀 ̿ = [ 𝑀 𝑥 0 0 0 𝑀 𝑦 0 0 0 𝑀 𝑧 ] ; 𝐶 ̿ = [ 𝐶 𝑥 0 0 0 𝐶 𝑦 0 0 0 𝐶 𝑧 ] ; 𝐾 ̿ = [ 𝐾 𝑥 0 0 0 𝐾 𝑦 0 0 0 𝐾 𝑧 ]
There, vibratory milling dynamics subsystem is built from transfer functions of decoupled system on a modal basis. Indeed, chatter milling dynamics is represented by motion equations of three oscillators in which each one is vibrating along one of spindle reference frame directions.

As we can see, factors of directional dynamic coefficients matrix 𝛼 ̿(𝑡) are periodic at rotation frequency as milling tool rotates and engages the same elementary cutting edges in each interaction workpiece/tool. Then, based on immersion conditions and number of elemental cutting edges engaged, time varying dynamic coefficients can be expanded by their Fourier series components. It is straightforward to approximate the matrix directional coefficients by the constant term in the Fourier series. Then the expression (7) becomes:

𝑀 ̿ [ 𝑋 ̈(𝑡) 𝑌 ̈(𝑡) 𝑍 (𝑡) ] + 𝐶 ̿ [ 𝑋 ̇(𝑡) 𝑌 ̇(𝑡) 𝑍 (𝑡) ] + 𝐾 ̿ [ 𝑋(𝑡) 𝑌(𝑡) 𝑍(𝑡) ] = 𝛼 0 ̿̿̿ 2 * ∆ ̅ (𝑡) (8)
Where:

𝛼 0 ̿̿̿ = 1 𝑇 ∫ 𝛼 ̿(𝑡)𝑑𝑡 𝑇 0 𝑇 = 1 𝑁
The constant matrix directional coefficients 𝛼 0 ̿̿̿ is time invariant but engagement dependent. So, each factor 𝛼 𝑖𝑗 0 can be expanded as:

𝛼 𝑖𝑗 0 = ∑ ∑ 1 𝑇 ∫ 𝛼 𝑖𝑗 𝑑𝑡 𝑇 0 𝑧 𝐿=1 𝑛𝑘 𝑘=1 𝑖, 𝑗 ∈ {𝑥, 𝑦, 𝑧} (9)
IV. BEHAVIOR PREDICTION Hale and Lunel [START_REF] Hale | Introduction to Functional Differential Equations[END_REF] have defined five types of stability. On determination of milling process stability, a criterion is mandatory. As we have seen, chatter milling dynamics is modeled by system retarded differential equations where coefficients of vibratory displacements and its derivatives are classified into categories of constant and space-varying parameters. On the assumption of constant machining conditions and system parameters dynamics, the only coefficients related to kinematic collision between tool and workpiece are matrix 𝛼 0 ̿̿̿ coefficients. This is why stability analysis is mandatory in each interaction workpiece/tool. Of course, the literature suggests several methods to study dynamic systems stability involving delay [START_REF] Stépán | Retarded dynamical systems: stability and characteristic functions[END_REF], [START_REF] Hsu | Application of the Tau-Decomposition Method to Dynamical Systems Subjected to Retarded Follower Forces[END_REF]- [START_REF] Ruan | On the zeros of transcendental functions with applications to stability of delay differential equations with two delays[END_REF]]. We note, moreover, approaches based on the analysis of characteristic values associated to differential equations system modelling delayed process. Hence, Pontryagin method is used to determine asymptotic stability in Lyapunov sense of chatter milling process from characteristic function zeros of differential equations system.

Stability: types and criteria

Dynamic system stability can be determined from its related characteristic function zeros based on the following theorem [START_REF] Bellman | Differential-Difference Equations[END_REF], [START_REF] Wouk | Introduction to the Theory of Differential Equations with Deviating Arguments (L. E. El'sgol'ts)[END_REF] : Theorem 1. If all roots have negative real parts, then the trivial solution of ( 8) is asymptotically stable. Theorem 2. If at least one root has a positive real part, then the trivial solution of ( 8) is unstable. Theorem 3. If there are simple purely imaginary roots and the remaining roots have negative real parts, then the trivial solution of ( 8) is stable. Theorem 4. If there is a multiple root among the purely imaginary roots, then the trivial solution of ( 8) is unstable. From these theorems, characterization of chatter process stability depends on computing and evaluating all characteristic function roots. Remember also that chatter characteristic function is an exponential polynomial because of Laplace transform of delayed functions. Therefore, stability system evaluation is time consuming. So, Quasi-Polynomial Mapping Based Rootfinder (QPmR) is proposed to decrease the amount of computing power and rendering time. In a similar way to that of Routh-Hurwitz criterion applied on ordinary differential equations systems, QPmR is a function for computation and analysis of spectrum of characteristic quasi-polynomials of both retarded and neutral time delay systems [START_REF] Vyhlídal | Quasipolynomial mapping based rootfinder for analysis of time delay systems[END_REF]- [START_REF] Vyhlidal | Mapping Based Algorithm for Large-Scale Computation of Quasi-Polynomial Zeros[END_REF].

Stability computing

A frequency approach is proposed to evaluate stability of system based on analyzing the characteristic function of differential equations system modeling chatter system. Therefore, a new delayed differential equations system is redefined from the following dimensionless parameters:

𝜏 = 𝑡 𝑟 ; 𝑟 = 1 𝜐 ; λ x = √ K x M x ⁄ 2πυ ; λ y = √ K y M y ⁄ 2πυ ; λ z = √ K z M z ⁄ 2πυ ζ x = C x 2√K x M x ; ζ y = C y 2 √ K y M y ; 𝜁 𝑧 = 𝐶 𝑧 2√𝐾 𝑧 𝑀 𝑧
With 𝜏 dimensionless time parameter,𝜆 𝑥 , 𝜆 𝑦 , 𝜆 𝑧 dimensionless frequencies and 𝜁 𝑥 , 𝜁 𝑦 , 𝜁 𝑧 damping ratios along the tree axes x, y and z, respectively, of Galilean reference frame. Thanks to Laplace transform, RDEs system (8) is expressed in frequency domain [START_REF] Altintas | In-process detection and suppression of chatter in milling[END_REF] where 𝑠 is the complex number and 𝐼𝑑 the identity matrix.

𝐼𝑑 𝑠 2 + 4𝜋 [ 𝜆 𝑥 𝜁 𝑥 0 0 0 𝜆 𝑦 𝜁 𝑦 0 0 0 𝜆 𝑧 𝜁 𝑧 ] s + 4π 2 [ 𝜆 𝑥 2 0 0 0 𝜆 𝑦 2 0 0 0 𝜆 𝑧 2 ] + 4𝜋 2 [ 𝜆 𝑥 2 𝐾 𝑥 𝛼 𝑥𝑥 0 𝜆 𝑦 2 𝐾 𝑦 𝛼 𝑥𝑦 0 𝜆 𝑧 2 𝐾 𝑧 𝛼 𝑥𝑧 0 𝜆 𝑥 2 𝐾 𝑥 𝛼 𝑦𝑥 0 𝜆 𝑦 2 𝐾 𝑦 𝛼 𝑦𝑥 0 𝜆 𝑧 2 𝐾 𝑧 𝛼 𝑦𝑧 0 𝜆 𝑥 2 𝐾 𝑥 𝛼 𝑧𝑥 0 𝜆 𝑦 2 𝐾 𝑦 𝛼 𝑧𝑦 0 𝜆 𝑧 2 𝐾 𝑧 𝛼 𝑧𝑧 0 ] (𝑒 -𝑠 -1) = 0 (10)
Characteristic function P(s) of delayed vibratory milling system (8) is defined as: |𝑖, 𝑗 ∈ {𝑥, 𝑦, 𝑧}}. Then, stability of RDEs system (8) in the Lyapunov sense can be evaluated on the basis of characteristic function roots. In order to achieve this, we use Matlab software to develop an algorithm (Error! Reference source not found.) where we use the function QPmR [START_REF] Vyhlidal | Mapping Based Algorithm for Large-Scale Computation of Quasi-Polynomial Zeros[END_REF], [START_REF] Vyhlídal | QPmR -Quasi-Polynomial Root-Finder: Algorithm Update and Examples[END_REF] to compute characteristic quasi-polynomials roots of RDEs systems. Finally, stability milling process is determined by applying the above theorems. ̿̿̿ matrix elements, we need data related to cutting process parameters and tool geometry. Also, cutting system dynamics will allow us to get a closer look at its stability. For this, we exploit the experimental conditions of Li et al. [START_REF] Li | Cutting force prediction and analytical solution of regenerative chatter stability for helical milling operation[END_REF] where a helical milling tool having a helical angle of 30°, a diameter of 10 mm and two flutes is grooving on an AL7075-T6 aluminum alloy part withconstant feed 𝑓 𝑧 = 0.05 𝑚𝑚 𝑑𝑒𝑛𝑡 ⁄ . Cutting law parameters and dynamics parameters associated to milling system are defined as: A Two-DOF structural dynamics of milling system is considered. These two degrees of freedom are the displacements 𝑋(𝑡) and 𝑌(𝑡) in the two orthogonal directions, respectively 𝑋 and 𝑌 ⃗ axes. Finally chatter vibration is described by two equations of motion placed in the following matrix form:

P(s) = | | s 2 + 4πλ x ζ x s + 4π 2 λ x 2 + 4π 2 λ x 2 α xx 0 K x (e -𝑠 -
[ 𝑋 ̈(𝑡) 𝑌 ̈(𝑡) ] + 4𝜋 [ 𝜁 𝑥 𝜔 𝑥 0 0 𝜁 𝑦 𝜔 𝑦 ] [ 𝑋 (𝑡) 𝑌 ̇(𝑡) ] + 4𝜋 2 [ 𝜔 𝑥 2 0 0 𝜔 𝑦 2 ] [ 𝑋(𝑡) 𝑌(𝑡) ] = 2𝜋 2 [ 𝜔 𝑥 2 𝛼 𝑥𝑥 0 𝐾 𝑥 𝜔 𝑥 2 𝛼 𝑥𝑦 0 𝐾 𝑦 𝜔 𝑥 2 𝛼 𝑦𝑥 0 𝐾 𝑥 𝜔 𝑥 2 𝛼 𝑦𝑦 0 𝐾 𝑦 ] ∆ ̅ (𝑡)
Under these conditions, we define three configurations for the torque (rotation frequency, cutting depth). For each configuration, we examine roots of the associated characteristic function using a computer program written on Matlab script.

VI. RESULTS AND DISCUSSION

For each of the three configurations mentioned in Error! Reference source not found., we inject tool parameters, cutting law coefficients and process parameters associated with in the input of the script. After that, we evaluate the Lyapunov stability of the RDEs system [START_REF] Yao | On-line chatter detection and identification based on wavelet and support vector machine[END_REF]. Based on roots analysis of quasi-polynomial characteristic function, the QPmR function returns an asymptotic stability for configuration (1) and milling process instability for configurations (2) and (3). In order to evaluate QPmR function outputs, we use dde23 Matlab function to solve the system of DDEs with constant delays representing the two degree-of-freedom structural dynamics of milling system. Then, we obtain the signals of the regenerative milling system displacements with respect to the three configurations (1), ( 2), (3) on the axes 𝑋 and 𝑌 ⃗ .

Fig. 5: Vibratory displacements of the configurations (1), (2), (3)

Error! Reference source not found. proves the nature of milling system stability evaluated through the frequency approach. So, output signals of the configuration (1) are periodic at tooth passing and bounded. At the output of the configuration (2), we remark that the envelopes of the signals are linear. Nevertheless, at the simulation of the configuration (3), we obtain signals with exponential envelopes. In conclusion, the numerical results (extracts from the simulation of the regenerative milling) are qualitatively consistent with the analytical prediction (evaluation of the characteristic roots).

VII.

CONCLUSION Stability lobes find its limits in predicting stability of multi-dimensional chatter cutting process. Then, a new approach is proposed to determine chatter milling process stability in sense of Lyapunov by looking back to theorems of RDEs system stability in frequency domain. For a case study of milling process, three configurations of machining conditions are integrated in the inputs of the QPmR Matlab function to determine milling stability for each configuration. Impact of chatter on stability milling system is visualized in time domain by resolving the RDEs system thanks to the dde23 Matlab function.

VIII. REFERENCES

  ℎ 𝑘,𝐿 𝑠𝑡 = 𝑠𝑖𝑛 𝜅 𝑘 * 𝑠𝑖𝑛 𝜃 𝑘,𝐿 * 𝑓 𝑧 ℎ 𝑘,𝐿 = ℎ 𝑘,𝐿 𝑠𝑡 + 𝑠𝑖𝑛 𝜅 𝑘 * 𝑠𝑖𝑛 𝜃 𝑘,𝐿 * [𝛥𝑥(𝑡) + 𝛥𝑦(𝑡) * 𝑐𝑜𝑡 𝜃 𝑘,𝐿 -𝛥𝑧(𝑡) * 𝑐𝑜𝑡 𝜅 𝑘 ] ⏟

Fig. 4 :

 4 Fig. 4: Milling process behavior computation V. APPLICATION In the aim of calculating 𝛼 0̿̿̿ matrix elements, we need data related to cutting process parameters and tool geometry. Also, cutting system dynamics will allow us to get a closer look at its stability. For this, we exploit the experimental conditions of Li et al.[START_REF] Li | Cutting force prediction and analytical solution of regenerative chatter stability for helical milling operation[END_REF] where a helical milling tool having a helical angle of 30°, a diameter of 10 mm and two flutes is grooving on an AL7075-T6 aluminum alloy part withconstant feed 𝑓 𝑧 = 0.05 𝑚𝑚 𝑑𝑒𝑛𝑡 ⁄ . Cutting law parameters and dynamics parameters associated to milling system are defined as:

  𝐾𝑐 𝑟 , 𝐾𝑐 𝑡 , 𝐾𝑐 𝑎 Radial, tangential and axial cutting coefficients 𝐾𝑒 𝑟 , 𝐾𝑒 𝑡 , 𝐾𝑒 𝑎 Radial, tangentiel and axial de sharing coefficients

	Ω	Angular velocity
	N	rotation frequency
	𝑧	Teeth number
	𝑓 𝑧	Feed rate
	ℎ	Uncut chip thickness
	𝜃 𝑘,𝐿	Instantaneous engagement angle
	𝜅 𝑘	Axial immersion angle
	𝑑𝑏 𝑘	Elemental uncut chip width
	𝑑𝑆 𝑘	Elemental cutting edge length

  𝑘,𝐿 = { 𝑑𝐹𝑡 ̅̅̅̅̅ 𝑘,𝐿 = 𝐾𝑐 𝑡 * ℎ 𝑘,𝐿 * 𝑑𝑏 𝑘 + 𝐾𝑒 𝑡 * 𝑑𝑆 𝑘 𝑑𝐹𝑎 ̅̅̅̅̅ 𝑘,𝐿 = 𝐾𝑐 𝑎 * ℎ 𝑘,𝐿 * 𝑑𝑏 𝑘 + 𝐾𝑒 𝑎 * 𝑑𝑆 𝑘 𝑑𝐹𝑟 ̅̅̅̅̅ 𝑘,𝐿 = 𝐾𝑐 𝑟 * ℎ 𝑘,𝐿 * 𝑑𝑏 𝑘 + 𝐾𝑒 𝑟 * 𝑑𝑆 𝑘

	(2)
	}

  𝐾𝑒 𝑟 * 𝑠𝑖𝑛 𝜃 𝑘,𝐿 * 𝑠𝑖𝑛 𝜅 𝑘 -𝐾𝑒 𝑡 * 𝑐𝑜𝑠 𝜃 𝑘,𝐿 -𝐾𝑒 𝑎 * 𝑠𝑖𝑛 𝜃 𝑘,𝐿 * 𝑐𝑜𝑠 𝜅 𝑘 ) * 𝑑𝑆 𝑘 𝐾𝑒 𝑟 * 𝑐𝑜𝑠 𝜃 𝑘,𝐿 * 𝑠𝑖𝑛 𝜅 𝑘 + 𝐾𝑒 𝑡 * 𝑠𝑖𝑛 𝜃 𝑘,𝐿 -𝐾𝑒 𝑎 * 𝑐𝑜𝑠 𝜃 𝑘,𝐿 * 𝑐𝑜𝑠 𝜅 𝑘 ) * 𝑑𝑆 𝑘

	 𝛼 𝑋𝑍 = ∑ 𝑛𝑘 𝑘=1	∑ 𝑔(𝜃 𝑘,𝐿 ) * 𝑑𝑏 𝑘 * 𝑧 𝐿=1
	(𝐾𝑐 𝑟 * 𝑠𝑖𝑛 𝜃 𝑘,𝐿  𝛼 𝑌𝑋 = ∑ ∑ 𝑔(𝜃 𝑘,𝐿 ) * 𝑑𝑏 𝑘 * (-𝐾𝑐 𝑟 * 2 * 𝑠𝑖𝑛 2𝜅 𝑘 2 + 𝐾𝑐 𝑡 𝑠𝑖𝑛 2𝜃 𝑘,𝐿 * 𝑐𝑜𝑠 𝜅 𝑘 + 𝐾𝑐 𝑎 * 𝑠𝑖𝑛 𝜃 𝑘,𝐿 2 𝑠𝑖𝑛 2𝜃 𝑘,𝐿 * 𝑠𝑖𝑛 𝜅 𝑘 2 + 𝐾𝑐 𝑡 * 𝑠𝑖𝑛 𝜃 𝑘,𝐿 2 * 𝑐𝑜𝑠 𝜅 𝑘 ) 𝑧 𝐿=1 𝑛𝑘 𝑘=1 2 𝑠𝑖𝑛 𝜅 𝑘 -𝐾𝑐 𝑎 * 𝑠𝑖𝑛 2𝜃 𝑘,𝐿 2 * 𝑠𝑖𝑛 2𝜅 𝑘 ) 2  𝛼 𝑌𝑌 = ∑ ∑ 𝑔(𝜃 𝑘,𝐿 ) * 𝑑𝑏 𝑘 * (-𝐾𝑐 𝑟 * 𝑐𝑜𝑠 𝜃 𝑘,𝐿 2 * 𝑠𝑖𝑛 𝜅 𝑘 2 + 𝐾𝑐 𝑡 * 𝑠𝑖𝑛 2𝜃 𝑘,𝐿 2 * * 𝑧 𝐿=1 𝑛𝑘 𝑘=1 2 𝑠𝑖𝑛 𝜅 𝑘 -𝐾𝑐 𝑎 * 𝑐𝑜𝑠 𝜃 𝑘,𝐿 𝑠𝑖𝑛 2𝜅 𝑘 2 * 2
								𝑎 𝑋	𝛼 𝑋𝑋 𝛼 𝑋𝑌 𝛼 𝑋𝑍	𝛥𝑋(𝑡)	(6)
								𝑎 𝑌	} + [	𝛼 𝑌𝑋 𝛼 𝑌𝑌 𝛼 𝑌𝑍	] {	𝛥𝑌(𝑡)	}
								𝑎 𝑍	𝛼 𝑍𝑋 𝛼 𝑍𝑌 𝛼 𝑍𝑍	𝛥𝑍(𝑡)
	In order to simplify, we set:
					𝑎 ̅(𝑡) = {	𝑎 𝑋 𝑎 𝑌	} ; 𝛼 ̿(𝑡) = [	𝛼 𝑋𝑋 𝛼 𝑋𝑌 𝛼 𝑋𝑍 𝛼 𝑌𝑋 𝛼 𝑌𝑌 𝛼 𝑌𝑍	] ; ∆ ̅ (𝑡) = [	𝛥𝑋(𝑡) 𝛥𝑌(𝑡)	]
							𝑎 𝑍	𝛼 𝑍𝑋 𝛼 𝑍𝑌 𝛼 𝑍𝑍	𝛥𝑍(𝑡)
	Where :
	 𝑎 𝑋 = ∑ 𝑛𝑘 𝑘=1	∑ 𝑔(𝜃 𝑘,𝐿 ) * 𝑧 𝐿=1
	[	𝑓 𝑧 * (-𝐾𝑐 𝑟 * 𝑠𝑖𝑛 𝜃 𝑘,𝐿 +(-] 2 * 𝑠𝑖𝑛 𝜅 𝑘 2 -𝐾𝑐 𝑡 𝑠𝑖𝑛 2𝜃 𝑘,𝐿 2 * 𝑠𝑖𝑛 𝜅 𝑘 -𝐾𝑐 𝑎 * 𝑠𝑖𝑛 𝜃 𝑘,𝐿 2 * 𝑠𝑖𝑛 𝜅 𝑘 ) * 𝑑𝑏 𝑘
	 𝑎 𝑌 = ∑ 𝑛𝑘 𝑘=1	∑ 𝑔(𝜃 𝑘,𝐿 ) * 𝑧 𝐿=1
	[	𝑓 𝑧 * (-𝐾𝑐 𝑟 * +(-] 𝑠𝑖𝑛 2𝜃 𝑘,𝐿 2 * 𝑠𝑖𝑛 𝜅 𝑘 2 + 𝐾𝑐 𝑡 * 𝑠𝑖𝑛 𝜃 𝑘,𝐿 2 * 𝑠𝑖𝑛 𝜅 𝑘 -𝐾𝑐 𝑎 * 𝑠𝑖𝑛 2𝜃 𝑘,𝐿 2 * 𝑠𝑖𝑛 2𝜅 𝑘 2 ) * 𝑑𝑏 𝑘
	 𝑎 𝑍 = ∑ 𝑛𝑘 𝑘=1  𝛼 𝑋𝑋 = ∑ 𝑛𝑘 𝑘=1 ∑ 𝑔(𝜃 𝑘,𝐿 ) * [ 𝑧 𝐿=1 ∑ 𝑔(𝜃 𝑘,𝐿 ) * 𝑑𝑏 𝑘 * 𝑓 𝑧 * (-𝐾𝑐 𝑟 * +(-𝐾𝑒 𝑟 * 𝑐𝑜𝑠 𝜅 𝑘 -𝐾𝑒 𝑎 * 𝑠𝑖𝑛 𝜅 𝑘 ) * 𝑑𝑆 𝑘 𝑠𝑖𝑛 2𝜅 𝑘 2 * 𝑠𝑖𝑛 𝜃 𝑘,𝐿 -𝐾𝑐 𝑎 * 𝑠𝑖𝑛 𝜅 𝑘 2 * 𝑠𝑖𝑛 𝜃 𝑘,𝐿 ) * 𝑑𝑏 𝑘 𝑧 𝐿=1	]
	(-𝐾𝑐 𝑟 * 𝑠𝑖𝑛 𝜃 𝑘,𝐿  𝛼 𝑋𝑌 = ∑ ∑ 𝑔(𝜃 𝑘,𝐿 ) * 𝑑𝑏 𝑘 * 2 * 𝑠𝑖𝑛 𝜅 𝑘 2 -𝐾𝑐 𝑡 𝑧 𝐿=1 𝑛𝑘 𝑘=1	𝑠𝑖𝑛 2𝜃 𝑘,𝐿 2	* 𝑠𝑖𝑛 𝜅 𝑘 -𝐾𝑐 𝑎 * 𝑠𝑖𝑛 𝜃 𝑘,𝐿	2 * 𝑠𝑖𝑛 𝜅 𝑘 )
	(-𝐾𝑐 𝑟 *	𝑠𝑖𝑛 2𝜃 𝑘,𝐿 2	* 𝑠𝑖𝑛 𝜅 𝑘	2 -𝐾𝑐 𝑡 * 𝑐𝑜𝑠 𝜃 𝑘,𝐿	2 * 𝑠𝑖𝑛 𝜅 𝑘 -𝐾𝑐 𝑎 *	𝑠𝑖𝑛 2𝜃 𝑘,𝐿 2	* 𝑠𝑖𝑛 𝜅 𝑘 )

  Delayed function is converted to exponential function through Laplace transform. Then, characteristic function 𝑃(𝑠) is an exponential polynomial where 𝑠 6 𝑒 3 is the principal term.Coefficients of the characteristic quasi-polynomial 𝑃(𝑠) are formulas which call parameters related to structural dynamics of milling system besides the coefficients matrix 𝛼 0 ̿̿̿ of vibratory displacements. So, characteristic

	function	coefficients	determination	are		based	on	calculating	𝛼 0 ̿̿̿
	components {𝛼 𝑖𝑗	0 = ∑ 𝑛𝑘 𝑘=1	∑ 𝑧 𝐿=1	1 𝑇 ∫ 𝛼 𝑖𝑗 𝑑𝑡 𝑇 0
		4π 2 λ y K y 2 α yx	0	(e -s -1)	1)	4π 2 λ x K x 2 α xy s 2 + 4πλ y ζ y s + 4π 2 λ y 2 + 0 (e -𝑠 -1) 4π 2 λ y 2 α 𝑦𝑦 K y	0	(e -s -1)	4π 2 λ x K x 2 α xz 4π 2 λ y K y 2 α yz	0 0	(e -s -1) (e -s -1)	|
		4π 2 λ z K z 2 α zx	0	(e -s -1)			4π 2 λ z K z 2 α 𝑧𝑦	0	(e -s -1)	s 2 + 4πλ z ζ z s + 4π 2 λ z	2 +	4π 2 λ z K z 2 α zz	0	(e -s -1)	|

Table 1 : Milling process configurations

 1

	Configurations	Rotation frequency [𝑡𝑟 𝑚𝑖𝑛 ⁄	]	Cutting depth [𝑚𝑚] (𝑑𝑏 𝑘 * 𝑛𝑘)
	(1)	1800		𝑑𝑏 𝑘 = 0.15 𝑚𝑚 ; 𝑛𝑘 = 20
	(2)	5000		𝑑𝑏 𝑘 = 0.05 𝑚𝑚 ; 𝑛𝑘 = 33
	(3)	10000		𝑑𝑏 𝑘 = 0.1 𝑚𝑚 ; 𝑛𝑘 = 25

  [1] R. N. Arnold, « Cutting Tools Research: Report of Subcommittee on Carbide Tools: The Mechanism of Tool Vibration in the Cutting of Steel », Proc. Inst. Mech. Eng., vol. 154, n o 1, p. 261-284, juin 1946.

ℎ 𝑘,𝐿 𝑑 (𝑡)


Zebo Zhou

Jin Wu

Jinling Wang
email: jinling.wang@unsw.edu.au

Hassen Fourati
email: hassen.fourati@gipsa-lab.grenoble-inp.fr

Optimal, Recursive and Sub-Optimal Linear Solutions to Attitude Determination from Vector Observations for GNSS/Accelerometer/Magnetometer Orientation Measurement

Keywords: attitude determination, GNSS receiver, Wahba's problem, vector observations, autonomous navigation

The integration of the Accelerometer and Magnetometer (AM) provides continuous, stable and accurate attitude information for land-vehicle navigation without magnetic distortion and external acceleration. However, magnetic disturbance and linear acceleration strongly degrade the overall system performance. As an important complement, the Global Navigation Satellite System (GNSS) produces the heading estimates, thus it can potentially benefit the AM system. Such a GNSS/AM system for attitude estimation is mathematically converted to a multi-observation vector pairs matching problem in this paper. The optimal and sub-optimal attitude determination and their time-varying recursive variants are all comprehensively investigated and discussed. The developed methods are named as the Optimal Linear Estimator of Quaternion (OLEQ), Suboptimal-OLEQ (SOLEQ) and Recursive-OLEQ (ROLEQ) for different application scenarios. The theory is established based on our previous contributions, and the multi-vector matrix multiplications are decomposed with the eigenvalue factorization. Some analytical results are proven and given, which provides the reader with a brand new viewpoint of the attitude determination and its evolution. With the derivations of the two-vector case, the n-vector case is then naturally formed. Simulations are carried out showing the advantages of the accuracy, robustness and time consumption of the proposed OLEQs, compared with representative methods. The algorithms are then implemented using the C++ programming language on the designed hardware with a GNSS module, three-axis accelerometer and three-axis magnetometer, giving an effective validation of them in real-world applications. The designed schemes have proven their fast speed and good accuracy in these verification scenarios.

Introduction

Attitude determination (or estimation) from vector observation pairs is a significant technology in aerospace engineering and related geodetic applications [START_REF] Zhou | Equality constrained robust measurement fusion for adaptive kalman-filter-based heterogeneous multi-sensor navigation[END_REF][START_REF] Zhou | Integrated Navigation System for a Low-Cost Quadrotor Aerial Vehicle in the Presence of Rotor Influences[END_REF][START_REF] Chang | M-estimator for the 3D symmetric Helmert coordinate transformation[END_REF]. For instance, inertial navigation, having an important role in modern military applications, has a high demand for attitude determination accuracy for the initial alignment [START_REF] Wu | A New Technique for INS/GNSS Attitude and Parameter Estimation Using Online Optimization[END_REF][START_REF] Chang | Pseudo Open-Loop Unscented Quaternion Estimator for Attitude Estimation[END_REF][START_REF] Han | Quantization and colored noises error modeling for inertial sensors for GPS/INS integration[END_REF][START_REF] Chang | Initial Alignment by Attitude Estimation for Strapdown Inertial Navigation Systems[END_REF]. A typical attitude-measuring system integrating a three-axis Accelerometer with a three-axis Magnetometer (AM) is extensively applied for real-time, continuous, stable and accurate attitude estimation for various navigation campaigns [START_REF] Marins | An extended Kalman filter for quaternion-based orientation estimation using MARG sensors[END_REF][START_REF] Tian | An adaptive-gain complementary filter for real-time human motion tracking with MARG sensors in free-living environments[END_REF]. For most applications, AM sensors are very efficient for low-cost attitude determination [START_REF] Wu | Attitude Determination Using a Single Sensor Observation : Analytic Quaternion Solutions and Property Discussion[END_REF]. However, the magnetometer is easily disturbed by unknown and unexpected magnetic fields in electromagnetic signal-contaminated environments. On the other hand, for large-scale region applications, the reference magnetic vector is no longer a constant vector and needs to be corrected in a timely manner with other additional heading information. Otherwise, the overall system performance will be heavily degraded. Moreover, accelerometers inevitably suffer from biases, thus leading to inaccurate attitude estimation. Therefore, auxiliary sensors are necessary to overcome such a problem.

Alternatively, the Global Navigation Satellite System (GNSS) provides precise position and velocity information [START_REF] Du | Inertial aided cycle slip detection and identification for integrated PPP GPS and INS[END_REF]. It has been successfully used for land, marine and aircraft attitude determination applications [START_REF] Aleshechkin | Algorithm of GNSS-based attitude determination[END_REF]. Traditional methods integrate GNSS with inertial sensors and simultaneously estimate the orientation with synchronized position, velocity and attitude loops [START_REF] Grewal | Kalman Filtering: Theory and Practice Using MATLAB[END_REF][START_REF] Niu | An Accurate Land-Vehicle MEMS IMU/GPS Navigation System Using 3D Auxiliary Velocity Updates[END_REF][START_REF] Niu | A universal approach for processing any MEMS inertial sensor configuration for land-vehicle navigation[END_REF]. However, this leads to the risk of degrading the accuracy of the attitude solution associated with position and velocity estimation loops. Thus, Gebre-Egziabher and Elkaim [START_REF] Gebre-Egziabher | MAV attitude determination by vector matching[END_REF] proposed an independent attitude estimation loop by means of vector matching. Compared with GNSS antenna arrays, which compute highly precise solutions of baselines by using carrier-phase measurement, a single GNSS antenna is preferred in many low-cost and low-power consumption platforms. It mainly uses the simultaneous velocity vector information generated by GNSS Doppler [START_REF] Zhang | Short note: On the relativistic Doppler effect for precise velocity determination using GPS[END_REF][START_REF] Zhou | Optimal Doppler-aided smoothing strategy for GNSS navigation[END_REF]. Indirectly, this produces important complementary orientation information, thus potentially benefiting the AM sensor system, especially for land vehicles. In addition, integrating high-rate sensors also contributes to seismogeodetic systems [START_REF] Geng | A new seismogeodetic approach applied to GPS and accelerometer observations of the 2012 Brawley seismic swarm: Implications for earthquake early warning[END_REF][START_REF] Geng | Recovering coseismic point ground tilts from collocated high-rate GPS and accelerometers[END_REF].

Efficient attitude estimation strategy is very crucial for a GNSS/AM integrated multi-sensor system. In essence, it can be mathematically converted into a multi-observation vector pairs matching problem. Representative methods are mainly about solutions to the famous Wahba's problem [START_REF] Wahba | A Least Squares Estimate of Satellite Attitude[END_REF], which aims to obtain the optimal attitude determination results using weighted least-squares. Among these algorithms, Shuster's QUaternionESTimator(QUEST [START_REF] Shuster | Three-axis attitude determination from vector observations[END_REF]), Markley's Singular Value Decomposition (SVD [START_REF] Markley | Attitude Determination using Vector Observations and the Singular Value Decomposition[END_REF]) and Mortari's algorithm (ESOQ [START_REF] Mortari | A Closed-Form Solution to the Wahba Problem[END_REF]) are the most frequently-used ones, which are mostly inspired by Davenport's q-method [START_REF] Davenport | A Vector Approach to the Algebra of Rotations with Applications[END_REF][START_REF] Chang | Error analysis of Davenport's q method[END_REF]. Our newly-developed Fast Linear Attitude Estimator (FLAE [START_REF] Wu | Fast Linear Quaternion Attitude Estimator Using Vector Observations[END_REF]) obtains the fastest Wahba's solution, to our existing knowledge. Some other interesting approaches are proposed, as well, investigating the other aspects of this problem, e.g., Yang's analytical method, the Riemannian-manifold algorithm and Forbes' Linear Matrix Inequality (LMI) solution [START_REF] Yang | An analytic solution to Wahbas problem[END_REF][START_REF] Yang | Attitude determination using Newton's method on Riemannian manifold[END_REF][START_REF] Forbes | Linear-Matrix-Inequality-Based Solution to Wahba's Problem[END_REF].

There are still some weight-less algorithms for multi-sensor attitude determination. They are usually used in applications like vision attitude determination, where the a priori information of weights can hardly be accurately determined. For example, using the nonlinear special orthogonal group on SO3 [START_REF] Mahony | Nonlinear complementary filters on the special orthogonal group[END_REF], it is able to obtain the attitude quaternion from arbitrary pairs of vectors. Via optimization approaches like the Gradient-Decent Algorithm (GDA [START_REF] Madgwick | Estimation of IMU and MARG orientation using a gradient descent algorithm[END_REF]) and the Gauss-Newton Algorithm (GNA [START_REF] Tian | An adaptive-gain complementary filter for real-time human motion tracking with MARG sensors in free-living environments[END_REF]), we may also compute stable solutions. Apart from these, a famous analytical method was proposed by Arun et al. where Singular Value Decomposition (SVD) is employed to calculate the attitude matrix [START_REF] Arun | Least-Squares Fitting of Two 3-D Point Sets[END_REF]. Similarly, it is then introduced in computing both the attitude and translation vector in machine vision systems [START_REF] Ye | 6-DOF Pose Estimation of a Robotic Navigation Aid by Tracking Visual and Geometric Features[END_REF].

For Wahba's problem, it has been shown that most existing algorithms are based on Davenport's q-method. For a long time, the attitude-solving process has been fixed in this framework, which aims to find the largest eigenvalue of the Davenport matrix K. Is it possible to seek another quite different attitude determination approach? The answer is positive, and in this paper, three novel quaternion attitude determination algorithms from pairs of vector measurements are proposed in the sense of optimal, time-recursive and sub-optimal formulations. The main contributions are listed below:

1.

The main theory is based on our previous research contributions and is extended to linear arbitrary pairs vector measurements for GNSS/AM attitude application.

2.

Three estimators, i.e., the Optimal Linear Estimator of Quaternion (OLEQ), Recursive OLEQ (ROLEQ) and Sub-optimal OLEQ (SOLEQ), are derived. We also proposed acceleration techniques to make the algorithms faster.

3.

Simulations and real experiments are carried out, which verify the accuracy, flexibility, robustness and time consumption of various algorithms for GNSS/AM attitude determination. Detailed comparisons with representative methods demonstrate the superiority of the proposed OLEQs.

This paper is structured as follows: Section 2 briefly introduces the GNSS-AM sensor system along with its functional and stochastic models formulated by vector pair matching. Section 3 contains the problem formulation and starts with the quaternion estimation from a single sensor observation. In Section 4, the two-vector attitude determination theory along with the n-vector theory are given accordingly. Section 5 involves the numerical examples and the real field test where comparisons between the proposed SOLEQ and other representative methods are given. Finally, Section 6 consists of concluding remarks and future work.

Fundamentals of the GNSS/AM System

For attitude determination, we require GNSS receivers and AM sensor arrays for low-cost and accurate solutions. First, considering the motion behaviours of land vehicles, a vector pair for GNSS velocity can be established as:

v b G = C b r v r G + ε b G (1)
where v denotes velocity vector; the subscript G denotes the observation source 'GNSS'; the symbols r and b represent the navigation frame (r-frame, North-East-Down (NED)) and body frame (b-frame, forward-right-down), respectively; the vector is transformed from the r-frame to the b-frame; ε b G is the Gaussian white noise with variance of R G . For the land-vehicle application, it should be pointed out that the velocities in the r-frame and the b-frame are respectively:

v r = C r e v e , v b G =     (v e x ) 2 + v e y 2 + (v e z ) 2 0 0     (2)
where e denotes an Earth-Centred-Earth-Fixed (ECEF) coordinate frame (e-frame), i.e., WGS-84; C r e denotes the transformation from the e-frame to the r-frame, which can be computed according to the GNSS position and Earth ellipsoid metrics in advance.

AM sensors consist of a three-axis MEMS accelerometer and three-axis magnetometer. The accelerometer gives the specific force measurement of a rigid body, and the magnetometer provides the users with the sensed local Earth geomagnetic field according to the IGRFmodel [START_REF] Abdelrahman | Simultaneous spacecraft attitude and orbit estimation using magnetic field vector measurements[END_REF].

z A = C b r (µ r -G r ) + b A + ε A , ε A ∼ (0, R A ) z M = C b r M r + ε M , ε M ∼ (0, R M ) (3)
where the subscripts A and M denote the accelerometer and magnetometer sources, respectively; z denotes the observed vector; b denotes the accelerometer bias; G r = [0, 0, -g] T , where g is the gravity, is a function of the geo-location; the normalized Earth's magnetic field reference vector M r = [cosα, 0, -sinα] T , where α is the local dip angle; µ is the linear acceleration vector, which is usually treated as an external disturbance; ε is the Gaussian white noise with the variance of R. For simplicity, two points need to be clarified that: (1) the bias term has been obtained and compensated for the accelerometer readings. [START_REF] Zhou | Integrated Navigation System for a Low-Cost Quadrotor Aerial Vehicle in the Presence of Rotor Influences[END_REF] The linear acceleration µ is estimated by using GNSS velocity information with differentiation between two adjacent epochs. Then, the vectors and matrices for the system model in the b-frame and the r-frame can be merged into the following multi-vector pair equations:

D b i = C b r D r i + ε i , i = 1, 2, 3 (4)
where:

    ε 1 ε 2 ε 3    =    ε G ε A ε M    (5)
Neglecting the cross-correlations between sensors, the stochastic model of the system is:

   ε G ε A ε M    ∼ (0, R), R =    R G R A R M    (6)

Problem Formulation

The conventional Wahba's problem aims to find the optimal attitude matrix from vector observation pairs in the sense of least-squares, such that:

L(C) = n ∑ i=1 a i D b i -CD r i 2 , n = 2, 3, • • • (7)
in which C denotes the optimal Direction Cosine Matrix (DCM);

D b i = D b x,i , D b y,i , D b z,i ⊤ and D r i = D r x,i , D r y,i , D r z,i ⊤
are the i-th pair of normalized vector observations from the body frame b and the reference frame r, respectively; a i is the weight of the i-th sensor output, which is given by the standard deviations of the input vectors from σ 1 to σ n :

a i = 1 σ 2 i n ∑ k=1 1 σ 2 k ( 8
)
provided that the variance information such as shown in Equation ( 6) is predetermined [START_REF] Shuster | Three-axis attitude determination from vector observations[END_REF]. This ensures the biggest variance scalar with respect to the smallest weight. Wahba's problem is solved via many representative methods [START_REF] Markley | Quaternion attitude estimation using vector observations[END_REF]. Many of these algorithms solve the problem via eigenvalue decompositions analytically or numerically [START_REF] Yang | An analytic solution to Wahbas problem[END_REF][START_REF] Yang | Attitude determination using Newton's method on Riemannian manifold[END_REF]. When there is only one pair of vector observations, Wahba's solutions fail to obtain the optimal quaternion since there will be ambiguous quaternions corresponding to the maximum eigenvalue of one [START_REF] Markley | How to estimate attitude from vector observations[END_REF][START_REF] Markley | Humble problems[END_REF]. In our previous contribution [START_REF] Wu | Fast Complementary Filter for Attitude Estimation Using Low-Cost MARG Sensors[END_REF], the continuous stable quaternion solution to an accelerometer-based attitude determination system is derived.

Quaternion from a Single Sensor Observation

Considering an attitude determination model from a pair of vector observations:

D b = CD r (9)
this section deals with the attitude determination from a single pair of sensor observations. Note that the DCM is decomposed with quaternions q = (q 0 , q 1 , q 2 , q 3 ) ⊤ [START_REF] Wu | Fast Complementary Filter for Attitude Estimation Using Low-Cost MARG Sensors[END_REF] via: C = (P 1 q, P 2 q, P 3 q) [START_REF] Wu | Attitude Determination Using a Single Sensor Observation : Analytic Quaternion Solutions and Property Discussion[END_REF] in which:

P 1 =   
q 0 q 1 -q 2 -q 3 -q 3 q 2 q 1 -q 0 q 2 q 3 q 0 q 1    , P 2 =    q 3 q 2 q 1 q 0 q 0 -q 1 q 2 -q 3 -q 1 -q 0 q 3 q 2    , P 3 =    -q 2 q 3 -q 0 q q 1 q 0 q 3 q q 0 -q 1 -q 2 q    [START_REF] Du | Inertial aided cycle slip detection and identification for integrated PPP GPS and INS[END_REF] In this section, the theory is extended to an arbitrary sensor with a similar approach as in [START_REF] Wu | Fast Complementary Filter for Attitude Estimation Using Low-Cost MARG Sensors[END_REF]. Inserting Equation [START_REF] Wu | Attitude Determination Using a Single Sensor Observation : Analytic Quaternion Solutions and Property Discussion[END_REF] into Equation ( 9) yields:

D b = (D r x P 1 + D r y P 2 + D r z P 3 )q (12)
It has been proven in [START_REF] Wu | Fast Complementary Filter for Attitude Estimation Using Low-Cost MARG Sensors[END_REF] that:

P ⊤ 1 = P † 1 , P ⊤ 2 = P † 2 , P ⊤ 3 = P † 3 ( 13
)
where † stands for the Moore-Penrose pseudo-inverse. In fact, another property can be deduced from the following theorem:

Theorem 1. P 1 P ⊤ 2 + P 2 P ⊤ 1 = P 1 P ⊤ 3 + P 3 P ⊤ 1 = P 2 P ⊤ 3 + P 3 P ⊤ 2 = 0 3×3 .
Proof. This can be proven via brute-force calculation:

P 1 P ⊤ 2 =    0 -q 1 0 + q 2 1 + q 2 2 -q 2 3 2q 0 q 1 + 2q 2 q 3 q 1 0 -q 2 1 -q 2 2 + q 2 3 0 2q 0 q 2 -2q 1 q 3 -2q 0 q 1 -2q 2 q 3 -2q 0 q 2 + 2q 1 q 3 0    P 2 P ⊤ 1 =    0 q 1 0 -q 2 1 -q 2 2 + q 2 3 -2q 0 q 1 -2q 2 q 3 -q 1 0 + q 2 1 + q 2 2 -q 2 3 0 -2q 0 q 2 + 2q 1 q 3 2q 0 q 1 + 2q 2 q 3 2q 0 q 2 -2q 1 q 3 0    (14a)
P 1 P ⊤ 3 =    0 2q 0 q 1 -2q 2 q 3 q 1 0 -q 2 1 + q 2 2 -q 2 3 -2q 0 q 1 + 2q 2 q 3 0 -2q 1 q 2 -2q 0 q 3 -q 1 0 + q 2 1 -q 2 2 + q 2 3 2q 1 q 2 + 2q 0 q 3 0    P 3 P ⊤ 1 =    0 -2q 0 q 1 + 2q 2 q 3 -q 1 0 + q 2 1 -q 2 2 + q 2 3 2q 0 q 1 -2q 2 q 3 0 2q 1 q 2 + 2q 0 q 3 q 1 0 -q 2 1 + q 2 2 -q 2 3 -2q 1 q 2 -2q 0 q 3 0    (14b)
P 2 P ⊤ 3 =    0 2q 0 q 2 + 2q 1 q 3 -2q 1 q 2 + 2q 0 q 3 -2q 0 q 2 -2q 1 q 3 0 q 1 0 + q 2 1 -q 2 2 -q 2 3 2q 1 q 2 -2q 0 q 3 -q 1 0 -q 2 1 + q 2 2 + q 2 3 0    P 3 P ⊤ 2 =    0 -2q 0 q 2 -2q 1 q 3 2q 1 q 2 -2q 0 q 3 2q 0 q 2 + 2q 1 q 3 0 -q 1 0 -q 2 1 + q 2 2 + q 2 3 -2q 1 q 2 + 2q 0 q 3 q 1 0 + q 2 1 -q 2 2 -q 2 3 0    (14c) Lemma 1. Let K(q) = (D r x P 1 + D r y P 2 + D r z P 3 ), K ⊤ (q) = K † (q).
Proof. We have:

K(q)K ⊤ (q) = (D r x P 1 + D r y P 2 + D r z P 3 )(D r x P 1 + D r y P 2 + D r z P 3 ) ⊤ = (D r x ) 2 P 1 P ⊤ 1 + D r y 2 P 2 P ⊤ 2 + (D r z ) 2 P 3 P ⊤ 3 +D r x D r y (P 1 P ⊤ 2 + P 2 P ⊤ 1 ) + D r x D r z (P 1 P ⊤ 3 + P 3 P ⊤ 1 ) + D r y D r z (P 2 P ⊤ 3 + P 3 P ⊤ 2 ) = (D r x ) 2 + D r y 2 + (D r z ) 2 I 3×3 = I 3×3 (15)
which gives Lemma 1.

Hence, Equation ( 12) can be transformed into:

D b = K(q)q ⇒ K † (q)D b = q ⇒ K ⊤ (q)D b = q (16)
Since K(q) is in the form of a quaternion, K ⊤ (q)D b can be expanded by:

K ⊤ (q)D b = D r x P ⊤ 1 D b + D r y P ⊤ 2 D b + D r z P ⊤ 3 D b (17)
Furthermore, it can be obtained that:

P ⊤ 1 D b =      D b x q 0 + D b z q 2 -D b y q 3 D b x q 1 + D b y q 2 + D b z q 3 D b z q 0 + D b y q 1 -D b x q 2 -D b y q 0 + D b z q 1 -D b x q 3      =      D b x 0 D b z -D b y 0 D b x D b y D b z D b z D b y -D b x 0 -D b y D b z 0 -D b x      q = M 1 q (18a) P ⊤ 2 D b =      D b y q 0 -D b z q 1 + D b x q 3 -D b z q 0 -D b y q 1 + D b x q 2 D b x q 1 + D b y q 2 + D b z q 3 D b x q 0 + D b z q 2 -D b y q 3      =      D b y -D b z 0 D b x -D b z -D b y D b x 0 0 D b x D b y D b z D b x 0 D b z -D b y      q = M 2 q (18b) P ⊤ 3 D b =      D b z q 0 + D b y q 1 -D b x q 2 D b y q 0 -D b z q 1 + D b x q 3 -D b x q 0 -D b z q 2 + D b y q 3 D b x q 1 + D b y q 2 + D b z q 3      =      D b z D b y -D b x 0 D b y -D b z 0 D b x -D b x 0 -D b z D b y 0 D b x D b y D b z      q = M 3 q (18c)
Then, we have:

K ⊤ (q)D b = D r x M 1 q + D r y M 2 q + D r z M 3 q = Wq ( 19
)
where W is given by:

W = D r x M 1 + D r y M 2 + D r z M 3 (20)
Therefore, the attitude determination problem is shifted to:

Wq = q (21)
Theorem 2. For any normalized vector observation pair D r , D b , W 2 = I, where I is the four-dimensional identity matrix.

Proof. The characteristic polynomial of W is given by:

p(λ) = (λ -1) 2 (λ + 1) 2 (22)
Then, with the Cayley-Hamilton theorem, we substitute W for its eigenvalues and obtain:

W 2 -I 2 = 0 (23)
and as W 2 -I is real symmetric, it should be 0, which completes the proof.

According to our previous contribution [START_REF] Wu | Fast Complementary Filter for Attitude Estimation Using Low-Cost MARG Sensors[END_REF], Equation ( 21) can be treated as an iterative dynamical system:

q(n) = Wq(n -1) (24)
where q(n) denotes the quaternion for the n-th iteration. Furthermore, as has been proven, the discrete-time system can be converted to a corresponding continuous system. If W 2 = I, the stable solution to the continuous-time system is calculated as:

q = W + I 2 q rand ( 25
)
where q rand denotes a randomly-chosen unit quaternion. This provides us with a new approach to obtaining the measurement quaternion from a single strapdown sensor.

Optimal Linear Estimator of the Quaternion

With Equation ( 9), it is possible to list all the single rotation equations together as follows:

               D b 1 = CD r 1 D b 2 = CD r 2 . . . D b n = CD r n ⇒                √ a 1 D b 1 = √ a 1 CD r 1 √ a 2 D b 2 = √ a 2 CD r 2 . . . √ a n D b n = √ a n CD r n ( 26
)
which can in turn be converted to:

             √ a 1 q = √ a 1 W 1 q √ a 2 q = √ a 2 W 2 q . . . √ a n q = √ a n W n q (27)
This system is rewritten as:

       √ a 1 I √ a 2 I . . . √ a n I        q =        √ a 1 W 1 √ a 2 W 2 . . . √ a n W n        q (28)
Noticing that:

( √ a 1 I, √ a 2 I, • • • , √ a n I)        √ a 1 I √ a 2 I . . . √ a n I        = n ∑ i=1 a i I = I (29)
hence the pre-multiplication leads to:

q = n ∑ i=1 a i W i q (30)
It should be then noted that Equation ( 28) is free of noises. In this case, the maximum eigenvalue of n ∑ i=1 a i W i in engineering practice is very close to the noise-free theoretical result of one. Corresponding to Equation [START_REF] Mortari | A Closed-Form Solution to the Wahba Problem[END_REF], based on Brouwer's fixed-point theorem, it is possible to recursively obtain the normalized optimal quaternion by rotating a randomly given initial quaternion over and over again to infinity. This is something similar to the power method of the numerical eigenvalue factorization, but can be accelerated with geometric series-like iterations as follows:

R = n ∑ i=1 a i W i R 2 = R • R . . . R 2 j = R 2 j-1 • R 2 j-1 (31)
where R denotes the rotation operator over the Hamilton space H. In fact, the above iterations can hardly be achieved when the maximum eigenvalue approaches one. The reason is that at this time, the power R 2 approaches I, as well. A more robust way is shown to solve this problem. Considering both sides of Equation ( 28), we find out that the right side is in fact the mixture of solutions to single vector observation pairs. As mentioned in Equation ( 25), a stable and continuous solution to each single equation can be done by pre-multiplying 1 2

(W i + I). Substituting 1 2 (W i + I) for W i , i.e.,        √ a 1 I √ a 2 I . . . √ a n I        q =            1 2 √ a 1 (W 1 + I) 1 2 √ a 2 (W 2 + I) . . . 1 2 √ a n (W n + I)            q ( 32
)
the quaternion evolutes are from the mixture of the single optimal solution. Here, the rotation operator is revised as:

R = n ∑ i=1 1 2 a i (W i + I) = 1 2 I + n ∑ i=1 a i W i ( 33
)
This equals the least-square of the set of pre-computed single rotated quaternions, which is definitely faster and more robust than rotation from a randomly given initial quaternion.

The covariance of the calculated optimal quaternion is equal to that of QUEST under the condition of least-square optimality:

Σ Q,OLEQ = 1 4 n ∑ i=1 1 σ 2 i I - n ∑ i=1 D b D b ⊤ Σ q,OLEQ = [q OLEQ ] 0 0 1×3 0 3×1 Σ Q,OLEQ [q OLEQ ] ⊤ (34)
where Q = (q 1 , q 2 , q 3 ) ⊤ is the vector part of the quaternion.

[q] defines the following matrix:

[q] =      q 3 -q 2 q 1 q 0 q 2 q 3 -q 0 q 1 -q 1 q 0 q 3 q 2 -q 0 -q 1 -q

2 q 3      (35)

Variant 1: Recursive-OLEQ

We have seen from the above formulations that for each epoch, the vector observations are batch processed thoroughly. When used in aerospace electronic systems, the measured vector observation pairs in neighbouring time epochs are usually continuous since they are always smoothed by the sum filters and Low-Pass Filters (LPF). Therefore, with this consideration, the attitude quaternion can be propagated from the last estimated one using the rotation operator described before. In this way, the quaternions are recursively computed with much fewer computations, and the accuracy is maintained. Furthermore, for highly-reliable attitude determination systems, high-precision gyroscopes are usually employed. This provides us with a second-stage acceleration scheme, inspired by the conventional recursive algorithms like filter QUEST, REQUEST, etc. [START_REF] Shuster | Filter QUEST or REQUEST[END_REF], by which we may first rotate the estimated quaternion in the last time epoch with the zero-order angular transition matrix by: [START_REF] Markley | Quaternion attitude estimation using vector observations[END_REF] where:

q k = Φ k,k-1 q k-1
Φ k,k-1 = I + ∆t 2 [Ω×] (37)
in which ∆t is the sampling time and [Ω×] composed by the angular rate from the gyroscope, which is detailed in many classical literature works [START_REF] Valenti | A Linear Kalman Filter for MARG Orientation Estimation Using the Algebraic Quaternion Algorithm[END_REF]. After this, even a single rotation by R would then be very accurate. The one-step covariance matrix of the obtained quaternion is calculated by:

Σ q k ,ROLEQ = RΦ k,k-1 Σ q k-1 ,ROLEQ Φ ⊤ k,k-1 R ⊤ (38)

Variant 2: SOLEQ

Two-Vector Case

When there are two pairs of vector observations, regardless of the weights of the respective sensors, it is natural that the attitude quaternion may be computed via:

q = 1 4 (W 1 + I) (W 2 + I) q rand ( 39
)
where W 1 and W 2 are for the first and second sensors. Equation ( 25) is reasonable since we have:

1 2 (W + I) 2 = 1 4 W 2 + 2W + I = 1 4 (2W + 2I) = 1 2 (W + I) (40)
However, for the two-vector case, one can write:

1 4 (W 1 + I) (W 2 + I) 2 = 1 4 (W 1 W 2 + W 1 + W 2 + I) 2 = 1 4 (W 1 + I) (W 2 + I) (41)
This reflects that the accurate quaternion can be recursively computed with:

     q k = 1 4 (W 1 + I) (W 2 + I) q k-1 q k = q k q k ( 42
)
which exits when the Euler distance q kq k-1 is less than one predetermined threshold η. Note that the above initialization procedure is equivalent to the following process:

         q j = 1 4 (W 1 + I) (W 2 + I) j q rand q j = q j q j ( 43
)
where q rand is a randomly chosen quaternion; j is chosen to make sure q jq j-1 < η. Let us define the following transformation operators:

A = 1 4 (W 1 + I) (W 2 + I) (44a)
B = A ⊤ = 1 4 (W 2 + I) (W 1 + I) (44b)
Theorem 3. For the two-vector attitude determination case, the steady-state evolution in Equation ( 43) is not affected by the transformation operators' order, such that:

     q j = A j q rand q j = q j q j ⇔      q j = B j q rand q j = q j q j , j → +∞ (45)
Proof. The integrated transformation can be computed by:

AB = 1 16 (W 1 + I) (W 2 + I) 2 (W 1 + I) = 1 8 (W 1 + I) (W 2 + I) (W 1 + I) = 1 2 A (W 1 + I) = 1 2 (W 1 + I) B (46)
and accordingly

BA = 1 2 B (W 2 + I) = 1 2 (W 2 + I) A (47)
These generate:

ABA = 1 2 (W 1 + I) BA = 1 4 (W 1 + I) (W 2 + I) A = A 2 (48)
and:

ABAB = A 2 B = A (AB) = 1 2 A 2 (W 1 + I) (49)
Finally, we have:

(AB) j = 1 2 A j (W 1 + I) = 1 2 (W 1 + I) B j (50)
Then, U is computed by:

U = HH ⊤ =      0 0 0 0 0 0 0 0 0 0 u 1 u 2 0 0 u 3 u 4      =      0 0 0 0 0 0 0 0 0 0 h 2 1 + h 2 2 h 1 h + h 2 h 4 0 0 h 1 h 3 + h 2 h 4 h + h 2 4      (62)
where u 2 = u 3 . Using this, we have:

(AB) j = V 1 UV ⊤ 1 V 1 UV ⊤ 1 • • • V 1 UV ⊤ 1 = V 1 U j V ⊤ 1 ( 63
)
U can be decomposed with eigenvalue decomposition as well, such that:

U = V U S U V T U ( 64
)
which can be analytically given by:

S U = diag 0, 0, - 1 2 √ ∆ + u 1 + u 4 2 , + 1 2 √ ∆ + u 1 + u 4 2 V U =      1 0 0 0 0 1 0 0 0 0 -u 1 + u 4 + √ ∆ 2u 2 u 1 -u 4 + √ ∆ 2u 2 0 0 1 1      (65)
with:

∆ = u 2 1 -2u 1 u 4 + 4u 2 2 + u 2 4 (66)
Letting:

λ U,1 = - 1 2 √ ∆ + u 1 + u 4 2 λ U,2 = + 1 2 √ ∆ + u 1 + u 4 2 (67)
U j is finally computed by:

U j = V U S j V T U = V U diag(0, 0, λ j U,1 , λ j U,2 )V T U (68)
Therefore:

(AB) j = V 1 U j V ⊤ 1 = V 1 V U diag(0, 0, λ j U,1 , λ j U,2 )V ⊤ U V ⊤ 1 ( 69
)
The required computation of V i is computed as follows:

1 2 (W i + I) = Ṽi S Ṽ-1 i ( 70
)
where Ṽi (x, y) stands for the element of Ṽi in the x-th row and the y-th column, the details of which are given by Equation (71):

Ṽi (1, 1) = + 1 V i (D b x,i -D r x,i )(D b z,i -D r z,i ) Ṽi (1, 2) = - 1 V i (D b x,i -D r x,i )(D b y,i -D r y,i ) Ṽi (1, 3) = + 1 V i (D b x,i + D r x,i )(D b z,i + D r z,i ) Ṽi (1, 4) = - 1 V i (D b x,i + D r x,i )(D b y,i + D r y,i ) (71a) Ṽi (2, 1) = + 1 V i (D b x,i -D r x,i )(D b y,i + D r y,i ) Ṽi (2, 2) = + 1 V i (D b x,i -D r x,i )(D b z,i + D r z,i ) Ṽi (2, 3) = + 1 V i (D b x,i + D r x,i )(D b y,i -D r y,i ) Ṽi (2, 4) = + 1 V i (D b x,i + D r x,i )(D b z,i -D r z,i ) (71b) Ṽi (3, 1) = 1 Ṽi (3, 2) = 0 Ṽi (3, 3) = 1 Ṽi (3, 4) = 0 Ṽi (4, 1) = 0 Ṽi (4, 2) = 1 Ṽi (4, 3) = 0 Ṽi (4, 4) = 1 (71c)
in which the factor is computed by:

V i = D b y,i 2 + D b z,i 2 -D r y,i 2 -D r z,i 2 (72)
Related information can also be acquired from [START_REF] Wu | Fast Linear Quaternion Attitude Estimator Using Vector Observations[END_REF]. It should be noted that:

Ṽi Ṽ⊤ i = Ṽ⊤ i Ṽi = I (73)
Thus, the Gram-Schmidt orthogonalization is applied to Ṽi , enabling [42].

V i V i ⊤ = V i ⊤ V i = I
A typical commitment to achieve this is to implement the QR decomposition [START_REF] Yang | Spacecraft attitude determination and control: Quaternion based method[END_REF], such that:

V i R = Ṽi (74)
where R denotes an invertible upper triangular matrix. If q rand = (1, 0, 0, 0), the suboptimal quaternion is equal to the normalized first column of (AB) j .

n-Vector Case

Corresponding to the above notations and derivations, the n-vector case's transformation operators are defined by:

A = V 1 SV ⊤ 1 • • • V i SV ⊤ i • • • V n SV ⊤ n B = V n SV ⊤ n • • • V i SV ⊤ i • • • V 1 SV ⊤ 1 ( 75
)
Defining:

H = SV ⊤ 1 n-1 ∏ i=2 V i SV ⊤ i V n S = 1 2 n-2 SV ⊤ 1 n-1 ∏ i=2 (W i + I) V n S ( 76
)
we have:

U = SV ⊤ 1 • • • V i SV ⊤ i • • • V n S 2 V ⊤ n • • • V i SV ⊤ i • • • V 1 S = HH ⊤ (77)
Then:

AB = V 1 UV ⊤ 1 ⇒ (AB) j = V 1 U j V ⊤ 1 = V 1 V U diag(0, 0, λ j U,1 , λ j U,2 )V ⊤ U V ⊤ 1 (78)
Accordingly, the normalized first column of (AB) j constitutes the attitude quaternion.

The Effect of Power Order

In this sub-section, we will show that the selection of j is in fact not influential with respect to the final result at all. Letting:

g 1 = -u 1 + u 4 + √ ∆ 2u 2 g 2 = u 1 -u 4 + √ ∆ 2u 2 (79)
we have:

(AB) j (x, y) = V 1 V U diag(0, 0, λ j U,1 , λ j U,2 )V ⊤ U V ⊤ 1 = [V 1 (y, 4) + V 1 (y, 3)g 1 ] [V 1 (x, 4) + V 1 (x, 3)g 1 ] λ j U,1 + [V 1 (y, 4) + V 1 (y, 3)g 2 ] [V 1 (x, 4) + V 1 (x, 3)g 2 ] λ j U,2 (80)
Here, we also have:

λ U,2 > λ U,1 > 0 (81) since: u 1 + u 4 = h 2 1 + h 2 2 + h 2 3 + h 2 4 > 0 ( 82
)
and:

∆ -(u 1 + u 4 ) 2 = -4u 1 u 4 + 4u 2 2 = 4 (h 1 h 3 + h 2 h 4 ) 2 -h 2 1 + h 2 2 h 2 3 + h 2 4 = 4 h 2 1 h 2 3 + h 2 2 h 2 4 + 2h 1 h 2 h 3 h 4 -h 2 1 h 2 3 -h 2 2 h 2 3 -h 2 1 h 2 4 -h 2 2 h 2 4 = -4(h 2 h 3 -h 1 h 4 ) 2 < 0 (83)
Therefore, with increasing iteration numbers, the item multiplied by λ j U,1 gradually vanishes in the results. The limiting result of AB j turns out to be:

lim j→+∞ (AB) j (x, y) = [V 1 (y, 4) + V 1 (y, 3)g 2 ] [V 1 (x, 4) + V 1 (x, 3)g 2 ] λ j U,2 (84)
Additionally, the quaternion solution is not about which column of AB j , and the result is the normalization of the following vector:

⌢ q =      V 1 (1, 4) + V 1 (1, 3)g 2 V 1 (2, 4) + V 1 (2, 3)g 2 V 1 (3, 4) + V 1 (3, 3)g 2 V 1 (4, 4) + V 1 (4, 3)g 2      (85)

Discussion of OLEQs

The three derived OLEQs can be used on different occasions. The OLEQ incorporates the weights so that the determination results are optimal in the sense of least-squares. With the aid of the gyroscope, the ROLEQ can achieve faster and smoother estimates. The meaning of the proposed SOLEQ is that it has a very simple linear expression that may generate short and concise analytic results for certain sensor combinations. Furthermore, when required in applications where the weights can hardly be accurately determined, e.g., vision attitude determination, the SOLEQ could be an alternative choice, as well. The attitude determination results of the three OLEQs are evaluated in the following experimental section.

Simulations and Experiments

Simulation: Common Cases

In this sub-section, the sensor observations are simulated with random reference vectors and true DCM in:

D b = C true D r + ε ( 86
)
where ε is the noise item, which are not supposed to have any correlation with each other and be subject to Gaussian distribution. The reference vectors and the standard deviations of noise items are selected according to the classical test samples by Markley (see Table 1), where the reference DCM C true is: We first observe the attitude RMSEs. From the results of OLEQ, QUEST and FLAE, it is noticeable that they have a similar attitude determination accuracy. Combining the same statistics in Table 5, the proposed OLEQ is well verified for its optimality. From the presented results, we see that SOLEQ has larger attitude errors and loss function values compared to the other optimal methods. The proposed SOLEQ is sub-optimal as it actually approximates the attitude estimator where the weights are ignored. Therefore, this simulation scenario has validated the correctness and optimality of the proposed OLEQ and SOLEQ.

C true =    0.

Simulation: Extreme Cases

D r 1 =     1 0 0     , D r 2 =     -0.99712 0.07584 0     , D r 3 =     -0.99712 -0.07584 0     σ 1 = 1 arcsec, σ 2 = σ 3 = 1 deg (88)
Conventional Wahba's solutions face a dilemma when exposed to some critical pairs of vector observations. For instance, Markley and Mortari give an example where the sensors are configured by Equation (88) [START_REF] Markley | How to estimate attitude from vector observations[END_REF]. In such a scenario, the root of the characteristic polynomial of the Davenport matrix cannot be easily obtained by Newton iterations. The internal reason is given by Cheng [START_REF] Cheng | Improvement to the Implementation of the QUEST Algorithm[END_REF] showing that it results in numerical loss according to the specific CPU word storage length. A flexible transformation of the characteristic polynomial is proposed then to significantly boost the convergence. As such configurations indeed happen in engineering practice, there is a necessity to evaluate the proposed schemes by comparisons with representative solvers. With similar simulation techniques as previously mentioned, the vectors are simulated with given reference vectors and standard deviations by rotation of C true . Here, the QUEST algorithm is revised to Cheng's form. First, we mainly compare the two iterative methods QUEST and OLEQ because in our existing paper [START_REF] Wu | Fast Linear Quaternion Attitude Estimator Using Vector Observations[END_REF], QUEST and FLAE have been proven to have a very similar behaviour when faced with this extreme case. Here, the iteration stops when the Euclidean norm of the neighbouring attitude quaternion difference is less than 1 × 10 -8 . For QUEST, the maximum iteration number is set to 50. The obtained results are depicted in Figure 1. We can see that the supervised QUEST can obtain accurate quaternion solutions within several iterations. Actually, before Cheng's improvement, QUEST may exceed the maximum iteration number from time to time. The proposed OLEQ, however, shows better performance dealing with this extreme case. Furthermore, the final mean loss function values of the two algorithms are computed as 4.9890 × 10 -11 and 2.8391 × 10 -10 , which reveals that the proposed OLEQ can not only obtain faster solutions, but leads to smaller loss function values, compared with supervised QUEST. As is known to everyone, QUEST is the most representative Wahba's solution using Davenport's q-method. Many other algorithms like ESOQand FOAMactually have the same performance as QUEST. Therefore, in this way, OLEQ is proven to be faster and much more robust than the whole class of algorithms based on Davenport's q-method. This also shows that the presented novel attitude evolution method owns brand new abilities.

Experiment: Accelerometer-Magnetometer Case

In this sub-section, we conduct an experiment where the accelerometer-magnetometer combination is adopted. Such a sensor combination is extensively applied nowadays in low-cost attitude estimation schemes. The accelerometer is pre-calibrated using six-direction bias cancelling, while the magnetometer is calibrated using the method proposed by Y . Wu et al. [START_REF] Wu | On Calibration of Three-Axis Magnetometer[END_REF].

The hardware is constituted by a battery, a high-end Attitude and Heading Reference System (AHRS) with a high precision internal accelerometer, magnetometer and gyroscope, a transmitter for remote data transmission and a micro-controller for the implementation of the algorithm using the C++ language on the FreeRTOS. With the designed hardware platform shown in Figure 2, we collect a dataset with 10,000 samples.

The main purpose of this sub-section is to validate the performances of the proposed OLEQ, SOLEQ and ROLEQ since the AHRS has angular rate readouts. The compared results with the reference angles from representative methods are obtained (see Figures 345). Note that here, the weights between the accelerometer and magnetometer for Wahba's solution are chosen as 0.63 and 0.37 according to their respective noise characteristics. Yet, the local magnetometer's reference vector is calculated as M r = (0.60311, 0, -0.79766) ⊤ in Wuxi, Jiangsu Province, China. In Figure 3, the reference angles, QUEST solutions and SOLEQ solutions have been presented. The results indicate that the proposed sub-optimal estimator can estimate the attitude angles with similar macroscopic accuracy. Detailed attitude errors are plotted in Figures 4 and5. It is noticed that in general, these algorithms have the same errors with respect to the reference. Figure 5 shows that the overall attitude accuracy of ROLEQ is slightly smaller than the others. This is because it is first processed using the angular rate data, which essentially plays the role of a smoothing procedure. results using QUEST and SOLEQ.

Experiment: GNSS-Aided Attitude Determination for Land Vehicles

The GNSS receiver is widely employed in the attitude determination tasks for land and unmanned aerial vehicles. In this experiment, we use a designed rover (see Figure 6) to validate the feasibility of the proposed algorithm for GNSS attitude determination.

The rover is equipped with the aforementioned navigation computer and employs an external ublox M8N GNSS module with a serial-com connection to the board at the sampling frequency of 5 Hz. This rover is controlled by a handheld 2.4-GHz transmitter, and the onboard Pixhawk autopilot generates PID controlling commands to the servos and motors according to internal measurements. In this experiment, the rover is ran on a playground of UESTC, and we pick up one period of data in which the GNSS velocity is valid. In the data history, the magnetometer was distorted by outer unknown electromagnetic and ferromagnetic fields. Furthermore, during the execution process, raw measurements from the gyroscope and accelerometer are also logged with a speed of 250 Hz. The raw data are shown in Figure 7.

According to the sensor noise characteristics, the weights of the accelerometer-magnetometer combination are 0.9 and 0.1, respectively; while for the accelerometer-magnetometer-GNSS combination, the weights are chosen as 0.474, 0.05 and 0.474, respectively. The reference vector of the magnetometer is determined by the initial GNSS positioning results with the IGRF model [START_REF] Abdelrahman | Simultaneous spacecraft attitude and orbit estimation using magnetic field vector measurements[END_REF]. By making use of algorithms including QUEST, FLAE, OLEQ and SOLEQ, the computation results are summarized in Figures 8 and9.

We specifically add the GNSS velocity norm under the Euler angle results to illustrate the influence of the velocity scalar on the attitude determination results. In principle, when the vehicle is not moving with relative discriminative velocities, the GNSS receiver cannot give accurate speed estimates. Therefore, it is shown in the initial stage of the attitude determination results where the GNSS takes part, in that the determination accuracy of the yaw angle is seemingly very poor. As the velocity increases, the accuracy is improved very quickly accordingly. The accelerometer-magnetometer combination is largely distorted by magnetic disturbances. The integrated results of roll, pitch and, moreover, the yaw angles are influenced, generating very obvious differences with reference angles. With the aid of GNSS velocity, the corresponding attitude determination accuracy is not damaged because Wahba's solution balances the sensor inputs by the weights. It is observed that the OLEQ is validated to have almost the same accuracy for normal sensors in the aforementioned section, and in this section, such behaviour holds, as well. The SOLEQ, however, does not employ the weights and therefore produces relatively worse estimates, but for GNSS case, it is still better than the accelerometer-magnetometer estimates. The RMSEs of these scenarios are shown in Table 6. The results demonstrate to us the validity of the proposed algorithms, especially OLEQ.

Computation Time

From another point of view, the time consumption of various algorithms needs to be investigated. The time consumption is assessed on the embedded platform with the C++ programming language to ensure fairness. A rough evaluation is done with two pairs of vector observations in a few samples, which shows direct time consumption results (see Figure 10). As shown in the figure, for two pairs of vector observations, the three proposed algorithms' computation times are between QUEST and FLAE. However, from the expressions of the algorithms presented before, the number of vector observations is influential with respect to the final time consumption. Hence, with the simulation samples, each algorithm is again tested 20,000-times with different vector observation numbers. The time consumption is averaged and plotted in Figure 11. The results indicate that the algorithms are all linear, owning time complexities of O(n). QUEST, OLEQ and ROLEQ join at the vector observation number of 20. For common tasks, such a number covers most amounts of sensors. Although FLAE has the least time consumption, it can not overcome the drawbacks of extreme cases as well as OLEQ. That is to say, the proposed algorithms can replace the original algorithms for faster and more robust attitude determination.

Conclusions

This paper revisits the attitude determination from vector observations for a GNSS/AM case study. Novel linear algorithms are designed to obtain accurate attitude estimates in the sense of least-squares. Handling it in this manner, the computed quaternion is identical or suboptimal with respect to the conventional Wahba's solutions including QUEST and FLAE. We have found out that:

1.

Numerical simulations exhibit that the proposed OLEQs have a similar accuracy with representative solvers.

2.

It is also demonstrated that faced with extreme cases, OLEQs show much more robustness with less computation iterations.

3.

The computation speeds of OLEQs are tested, revealing that they belong to computationally-efficient algorithms. 4.

Moreover, a real vehicular experiment of a GNSS/AM system is designed and conducted showing the effectiveness of the proposed OLEQs in real-world embedded applications.

The MATLAB source codes of the proposed algorithms are uploaded as open-source files at https://github.com/zarathustr/OLEQs. The presented approach provides the reader with a brand new viewpoint of attitude evolution and hopefully will benefit related multi-sensor attitude determination applications.

Figure 1 .

 1 Figure 1. Iteration numbers of QUEST and OLEQ when faced with an extreme case.

Figure 2 .

 2 Figure 2. Designed hardware for the implementation of the algorithm.

Figure 3 .

 3 Figure 3. Experiment

Figure 4 .

 4 Figure 4. Experiment results of OLEQ, SOLEQ and representative algorithms using sampled data and different algorithms.

Figure 5 .

 5 Figure 5. Experiment results of Recursive OLEQ (ROLEQ) and representative algorithms using sampled data and different algorithms.

Figure 6 .

 6 Figure 6. The designed multi-functional rover for the validation of the proposed algorithms.

Figure 7 .

 7 Figure 7. Raw sensor measurements from the memory logging.

Figure 8 .

 8 Figure 8. Attitude determination from the accelerometer, magnetometer and velocity output of the GNSS receiver, by means of QUEST and FLAE.

Figure 9 .

 9 Figure 9. Attitude determination from the accelerometer, magnetometer and velocity output of the GNSS receiver, by means of OLEQ and SOLEQ.

Figure 10 .

 10 Figure 10. Time consumption of various algorithms.

Figure 11 .

 11 Figure 11. Time consumption of algorithms with respect to the number of vector observations.

Table 3 .

 3 Pitch RMSE (deg). 10 -05 5.7335 × 10 -05 4.0108 × 10 -05 4.0108 × -05 2 5.2860 × 10 -05 5.6736 × 10 -05 5.2860 × 10 -05 5.2860 × -05 3 4.0104 × 10 -01 5.6744 × 10 -01 4.0104 × 10 -01 4.0104 × -01 4 5.3887 × 10 -01 5.7656 × 10 -01 5.3887 × 10 -01 5.3887 × -01 5 3.9149 × 10 -01 3.9149 × 10 -01 4.4335 × 10 -01 1.2561 × +01 6 4.0121 × 10 -05 5.7809 × 10 -05 4.0121 × 10 -05 4.0121 × -05 7 5.3398 × 10 -05 5.7657 × 10 -05 5.3398 × 10 -05 5.3398 × -05 8 3.6755 × 10 -01 5.7326 × 10 -01 3.6755 × 10 -01 3.6755 × -01 9 4.5938 × 10 -01 5.7880 × 10 -01 4.5938 × 10 -01 4.5938 × -01 10 5.7186 × 10 -05 5.7186 × 10 -05 5.7184 × 10 -05 5.7186 × -05 11 5.7845 × 10 -05 5.7845 × 10 -05 5.7846 × 10 -05 5.7844 × -05 12 4.9161 × 10 -01 5.7554 × 10 -01 7.9376 × 10 +00 3.9359 × +00

	Case	OLEQ	SOLEQ	QUEST	FLAE
	1	4.0108 ×

Table 4 .

 4 Yaw RMSE (deg). 10 -05 6.1084 × 10 -05 4.3587 × 10 -05 4.3587 × -05 2 4.8694 × 10 -05 6.1015 × 10 -05 4.8694 × 10 -05 4.8694 × -05 3 4.4127 × 10 -01 6.0868 × 10 -01 4.4127 × 10 -01 4.4127 × -01 4 4.8593 × 10 -01 6.1289 × 10 -01 4.8593 × 10 -01 4.8593 × -01 5 2.5186 × 10 -01 2.5186 × 10 -01 2.8536 × 10 -01 1.7459 × +01 6 3.6421 × 10 -05 6.1651 × 10 -05 3.6421 × 10 -05 3.6421 × -05 7 4.8748 × 10 -05 6.0826 × 10 -05 4.8748 × 10 -05 4.8748 × -05 8 3.9812 × 10 -01 6.1557 × 10 -01 3.9812 × 10 -01 3.9812 × -01 9 4.9366 × 10 -01 6.1163 × 10 -01 4.9366 × 10 -01 4.9366 × -01 10 6.1834 × 10 -05 6.1834 × 10 -05 6.1836 × 10 -05 6.1844 × -05 11 6.2069 × 10 -05 6.2069 × 10 -05 6.2069 × 10 -05 6.2077 × -05 12 3.1726 × 10 -01 6.1275 × 10 -01 5.1968 × 10 +00 2.8814 × +00

	Case	OLEQ	SOLEQ	QUEST	FLAE
	1	4.3587 ×

Table 5 .

 5 Loss function values. 10 -13 1.0130 × 10 -12 5.0651 × 10 -13 5.0651 × -13 2 2.4901 × 10 -13 4.9802 × 10 -13 2.4901 × 10 -13 2.4901 × -13 3 4.9338 × 10 -05 9.8666 × 10 -05 4.9338 × 10 -05 4.9338 × -05 4 2.5369 × 10 -05 5.0736 × 10 -05 2.5369 × 10 -05 2.5369 × -05 5 5.0582 × 10 -13 5.0582 × 10 -13 6.5095 × 10 -13 8.5878 × -10 6 5.0422 × 10 -13 9.4333 × 10 -10 5.0422 × 10 -13 5.0422 × -13 7 2.4728 × 10 -13 8.8239 × 10 -09 2.4728 × 10 -13 2.4728 × -13 8 4.8216 × 10 -05 1.1593 × 10 -04 4.8216 × 10 -05 4.8216 × -05 9 2.5327 × 10 -05 5.0651 × 10 -05 2.5327 × 10 -05 2.5327 × -05 10 1.4827 × 10 -12 1.7575 × 10 -12 4.8431 × 10 -10 1.7195 × -12 11 4.8573 × 10 -13 4.8573 × 10 -13 2.4106 × 10 -10 9.2333 × -13 12 5.0105 × 10 -13 5.0105 × 10 -05 1.3143 × 10 -10 3.4336 × -11

	Case	OLEQ	SOLEQ	QUEST	FLAE
	1	5.0651 ×

Table 6 .

 6 RMSEs of the experiment (in deg).

	Roll Case	OLEQ SOLEQ QUEST FLAE
	Acc-MagExperiment	17.2	16.8	17.2	17.2
	Acc-Mag-GNSS Experiment	10.1	9.7	10.2	10.2
	Pitch Case	OLEQ SOLEQ QUEST FLAE
	Acc-Mag Experiment	15.9	14.6	15.9	15.9
	Acc-Mag-GNSS Experiment	14.4	13.2	14.4	14.4
	Yaw Case	OLEQ SOLEQ QUEST FLAE
	Acc-Mag Experiment	154.3	155.9	154.3	154.3
	Acc-Mag-GNSS Experiment	16.7	44.8	16.7	16.7

Remote Sens. 2018, 10, 377; doi:10.3390/rs10030377 www.mdpi.com/journal/remotesensing

Acknowledgments: This research is supported by the National Natural Science Foundation of China under Grant Numbers of 41604025, 41704029 and is also supported by the State Key Laboratory of Geodesy and Earth's Dynamics (Institute of Geodesy and Geophysics, CAS) Grant No. SKLGED2018-3-2-E. F. L. Markley provided constructive suggestions bout the internal theory of this paper. Yunheng Inc., Chengdu gave us the instrumentation services during the experimental validation. We genuinely thank them for their support.

Author Contributions: Zebo Zhou, Jin Wu, Jinling Wang and Hassen Fourati all conceived the idea of the paper, Zebo Zhou and Jin Wu conducted all the equation derivations and performed the experiments. Zebo Zhou and Jin Wu wrote the paper. Jinling Wang and Hassen Fourati improved the paper writing.

This proves that the mixed steady-state transformation (AB) j can be achieved by independent transformations from A or B, which completes the proof. Following this theorem, the problem is to compute the power A j . In fact, A is formed with 1 2 (W 1 + I) and 1 2 (W 2 + I). Their respective eigenvalue decompositions can be given by:

where V and D are constituted by eigenvectors and eigenvalues respectively as (W i + I) is real symmetric [START_REF] Horn | Matrix Analysis[END_REF]. Since W 1 and W 2 are in the same form, the eigenvalue matrices are equal to each other, i.e.,

Then, A can be rewritten as:

Identically, we have:

Combining Equations ( 53) and (54), it is obtained that:

Equation ( 55) is simplified as:

Here, we define:

Actually, it can be decomposed by:

An interesting fact is that the eigenvalue matrix S can be analytically calculated and is given by:

where diag(•) represents the diagonal matrix. This further yields H to be a matrix with the form of:

Case Reference Vectors

Noise Standard Deviations Using the simulated samples, the mean attitude Root Mean Squared Errors (RMSEs) in Euler angles are evaluated with our proposed algorithms OLEQ, SOLEQ and representative algorithms including QUEST and FLAE [START_REF] Wu | Fast Linear Quaternion Attitude Estimator Using Vector Observations[END_REF], which are shown in Tables 234. Table 5 contains the computed average Wahba's loss function values for different cases and algorithms. These algorithms are executed with MATLAB r2016 software (r2016, MathWorks, Natick, USA) on a PC 10,000 times with each data sample.