X Collilieux

Z Altamimi

P Rebischung

L Métivier

Coordinate kinematic models in the International Terrestrial Reference Frame releases

.

Introduction

Determining accurate coordinates at the Earth's surface is not straightforward since the reference system axes are not accessible in practice. As a consequence, a set of physical points whose coordinates are known is commonly used to materialize the terrestrial reference system and gives an implicit access to the system axes. This set of coordinates supplemented by statistical indicators is called Terrestrial Reference Frame. Those coordinates include a part of conventional information for defining the position of the origin of the frame, the orientation of its axes with respect to the crust and its scale (unit of length). One wishes that this frame be in co-rotation with the Earth for positioning objects at the Earth's surface.

As the Earth is not rigid, points at the Earth's surface exhibit relative motions which indubitably lead to coordinate variations in the Terrestrial Reference Frame whatever its definition is. Nowadays, relative positions between points are determined by space geodesy which provides the best precision for baselines at inter-continental distances. And all space geodetic measurements involve objects whose kinematic and dynamic descriptions are simplified in a celestial reference frame. Thus, the rotation of the Earth and coordinates of the tracking stations are solved for simultaneously in space geodetic software packages. This simultaneous estimation leads to a rank deficiency of the equation system to be solved since any change of the Earth's rotation vector can be canceled by an opposite rotation of the coordinate sets. Thus, constraints are added to define the reference frame of the station coordinates which make the estimated rotation vector of the Earth's motion in space fully dependent on the frame definition. [START_REF] Dermanis | Global reference frames: connecting observation to theory and geodesy to geophysics Dermanis, A[END_REF] revisited this issue of finding the most suitable reference frame in an elegant way. Assuming that the coordinates of station i, written x i (t), are known in an arbitrary frame, it is possible to define coordinates in a new frame x i '(t) using a similarity transformation x i '(t) = s(t) R(t) x i (t) + T(t). The fundamental problem of reference system consists in finding at each epoch t the scale factor s(t), the rotation matrix R(t) and the translation T(t) to be applied according to an optimal criterion. This criterion has been historically chosen in order to simplify the equation of the Earth's rotation in space [START_REF] Munk | The rotation of the Earth. A geophysical Discussion[END_REF][START_REF] Gross | Earth Rotation Variations -Long Period[END_REF] without introducing large coordinate changes in the terrestrial frames [START_REF] Dermanis | Global reference frames: connecting observation to theory and geodesy to geophysics Dermanis, A[END_REF]. Tisserand frames are ideal terrestrial frames since they fulfill this condition (Kovalevsky and Mueller, 1984). They have an origin at the Earth's center of mass and an orientation which is such that the contribution of the whole Earth's deformations to the angular momentum -the relative angular momentum -is equal to zero. Among all existing frames, these frames minimize the kinematic energy and thus involve less coordinate variations [START_REF] Boucher | Current intercomparisons between conventional terrestrial systems[END_REF].

In order to apply these concepts, the distribution of masses as well as the associated deformations need to be known at any point of the Earth's surface and of the Earth's interior. Whereas space geodesy is able to determine the Earth's center of mass position without this knowledge, this knowledge is required to define the Tisserand frame orientation time evolution. Also because geodetic stations are located at the Earth's surface, a modification of the Tisserand conditions has been adopted by the International Earth Rotation and Reference Systems Service. The International Terrestrial Reference System (ITRS) is specified as follows (Petit and Luzum, 2010, Chapter 4):

"1. It is geocentric, its origin being the center of mass for the whole Earth, including oceans and atmosphere; 2. The unit of length is the meter (SI). The scale is consistent with the TCG time coordinate for a geocentric local frame, in agreement with IAU and IUGG (1991) resolutions. This is obtained by appropriate relativistic modeling; 3. Its orientation was initially given by the BIH orientation at 1984.0; 4. The time evolution of the orientation is ensured by using a no-net-rotation condition with regards to horizontal tectonic motions over the whole Earth "

The last point is equivalent to nullify the relative angular momentum of the Earth's crust only. As discussed by Altamimi and Dermanis (2012b), this can be addressed by a kinematic approach or by an algebraic approach. The algebraic approach is the one that has been used to define the orientation of the ITRF. It relies on an a priori knowledge of the Earth's surface motion through a geophysical model and on the use of inner constraints for ITRF2000 (Altamimi et al., 2001). It also requires that the kinematic model for coordinate variations, namely piecewise linear functions, is known a priori. Conversely, Athanasios Dermanis has extensively developed the concept of the kinematic approach [START_REF] Dermanis | Global reference frames: connecting observation to theory and geodesy to geophysics Dermanis, A[END_REF][START_REF] Dermanis | Global reference frames: connecting observation to theory and geodesy to geophysics Dermanis, A[END_REF]Dermanis, 2003). This consists in the discrete computation of the relative angular momentum integral over the station network. In particular, he has developed the datum constraint equations by anticipating future changes in the ITRF kinematic model, considering that coordinate variations should not be linear anymore [START_REF] Dermanis | The ITRF Beyond the "Linear'' Model[END_REF]. Thus, polynomial functions, Fourier series and splines were suggested.

We propose in this paper to review the adopted model of ITRF2014 for coordinate variations and explain why some of [START_REF] Dermanis | The ITRF Beyond the "Linear'' Model[END_REF] suggestions are not yet implemented. In section I, ITRF previous releases are shortly described and the ITRF model is presented. Altamimi and Dermanis (2012a;2012b) theoretical work on datum constraint is recalled and a new interpretation of partial inner constraints involving transformation parameters is provided. In section II, the ITRF2014 model is discussed and finally new seasonal parameters are discussed in section III.

I ITRF as a parametric frame

I.1 Regularized frame

The instantaneous coordinates of a point in a terrestrial reference frame can be described as follows according to the IERS conventions (Petit and Luzum, 2010, Chapter 4):

x(t) = x r (t) + ΣΔx(t) (1)
where x r (t) are the coordinates at epoch t in the terrestrial reference frame, called regularized coordinates, and ΣΔx(t) is the sum of conventional corrections.

Up to now, the conventional corrections include displacements at rather short periods from solid Earth tides, ocean tide loading, S1 and S2 atmospheric loading and pole tides [START_REF] Petit | IERS Technical Note 36)[END_REF]IERS, 2015). Any other phenomenon that causes long-term crust motion will then be reported in the reference coordinates x r (t). Thus they reflect tectonic motions, post-glacial rebound or other effects.

ITRF coordinates are regularized coordinates described by simple mathematical functions as described in the next section. They are the result of an adjustment of space geodetic data including Global Navigation Satellite Systems (GNSS), Very long Baseline Interferometry (VLBI), Satellite and Lunar Laser Ranging (SLR and LLR) and Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS) [START_REF] Boucher | IERS Technical Note 27 -The 1997 International Terrestrial Reference Frame (ITRF97)[END_REF][START_REF] Altamimi | ITRF2000: A New Release of the International Terrestrial Reference Frame for Earth Science Applications[END_REF][START_REF] Altamimi | ITRF2005: A new release of the International Terrestrial Reference Frame based on time series of station positions and Earth Orientation Parameters[END_REF][START_REF] Altamimi | ITRF2008: an improved solution of the International Terrestrial Reference Frame[END_REF][START_REF] Altamimi | ITRF2014: A new release of the International Terrestrial Reference Frame modeling nonlinear station motions[END_REF]. Thus, a priori conventional models ΣΔx(t) have been included in the processing of the data which means that reference coordinates implicitly include conventions with respect to some adopted models. A relevant example is the handling of permanent deformations related to solid Earth tides (Petit and Luzum, 2010, chapter 1). A similar situation occurs with technique-specific corrections such as antenna phase center corrections for GNSS [START_REF] Schmid | Absolute IGS antenna phase center model igs08.atx: status and potential improvements[END_REF]. Any change in the used corrections impacts station coordinates, that's why the International GNSS Service (IGS) publishes updates of ITRF such as IGS08 for ITRF2008 every time a new correction file is released [START_REF] Rebischung | IGS08: the IGS realization of ITRF2008[END_REF].

From a user point of view, it means that terrestrial reference frame coordinates are related to certain conventions so that users should take care of the conventional corrections they add to them or use with them. For this reason, there is no attempt to correct technique-specific coordinate offsets in the published coordinates x r (t) if the conventional models are not changed.

II.2 History of ITRF products

There have been 13 ITRS realizations published since 1988 (Petit and Luzum, 2010, Chapter 4). Figure 1 provides a summary of the input data used in terms of space geodetic techniques and data span for the most recent published frames. As can be noticed, the input data span increases for each realization making the most recent ITRF superior to its predecessors.

fig.1 Time spans of recent ITRF input data

The kinematic model of coordinate variations x r (t) has evolved with time starting with constant coordinates in the earlier period. Although a geophysical velocity model was advised to be used to predict coordinates in the past and future, this method implied regular updates of the terrestrial frame since the geophysical model used differs from observations. From ITRF91 [START_REF] Boucher | ITRF 91 and its associated velocity field[END_REF], velocities were estimated and provided to the users. The following kinematic model for reference coordinates was adopted:

x r (t) = x c0 + Σδx i H(t-t x,i ) + (t -t 0 ) v c0 +Σδv i H(t-t v,i ) (2)
where H(t) is the Heaviside function (equal to 0 when t<0; equal to 1 when t≥0), t 0 is the reference epoch of the coordinates and t x,i t v,i are the epochs of discontinuities. In ITRF tables, coordinates are supplied to users as (x c (t), v c (t)) such as:

x c (t) = x c0 + Σδx i H(t-tx x,i ) v c (t)= v c0 + Σδv i H(t-tx x,i ) x r (t) = x c (t) + (t -t 0 ) v c (t) (3)
The dependency of the ITRF coordinates on the epoch is handled by means of the solution number, the so-called "soln". The soln is an integer which is associated to a time interval, so that any linear segment of the piece-wise linear function is identified by a unique integer.

Although computation methods have evolved, including the use of variance/covariance information since ITRF94 [START_REF] Boucher | IERS Technical Note 27 -The 1997 International Terrestrial Reference Frame (ITRF97)[END_REF], the use of time series of frames as input data since ITRF2005 [START_REF] Altamimi | ITRF2005: A new release of the International Terrestrial Reference Frame based on time series of station positions and Earth Orientation Parameters[END_REF], there has been no change in the description of station coordinate time evolutions up to ITRF2014.

II.3 Estimation model

II.3.1 Generalities

The ITRF estimation model for recent releases has been published in [START_REF] Altamimi | ITRF2000: A New Release of the International Terrestrial Reference Frame for Earth Science Applications[END_REF]2007) and more deeply conceptualized in Altamimi and Dermanis (2012b). In this section, we complete Altamimi and Dermanis (2012b) discussion by a complementary physical interpretation of the partial inner constraints also called internal constraints.

Since ITRF2005, time series of coordinates for a set of stations were provided for each measurement technique. Whereas all data could be in theory combined in one step (Equation 1of [START_REF] Altamimi | ITRF2005: A new release of the International Terrestrial Reference Frame based on time series of station positions and Earth Orientation Parameters[END_REF] to provide coordinates and velocities for each station in a well-defined reference frame, a twostep approach was adopted for better data editing and to facilitate the computation.

In the first step, positions, velocities and discontinuities in positions and velocities, hereafter called coordinates, are computed for every station for each technique independently. This first step is called stacking.

The second step, called combination, makes a proper average of all available data to provide coordinates in the same reference frame for all techniques making use of:

 coordinates of stations with their associated variance/covariance matrices from each technique as the results of the stacking step,  the relative positions of instruments from different techniques located in the same observatory -the so-called local tie vectors determined by topometry, levelling and GPS technique,  pseudo-observations to define the reference frame of the output coordinates.

II.3.2 Equation of stacking Equation of stacking

The observation equation of the stacking is given below for station i in the input station coordinate set s, following Altamimi and Dermanis (2012a) notations:

(4) where t k , s k , θ k are the transformation parameters, respectively a translation vector, a scale factor offset (with respect to 1) and a rotation vector -vectors are written in bold. Transformation parameters are assumed to be small so that the similarity mentioned in the introduction has been linearized around a priori position x i ap for station i. Thus, the rotation matrix becomes I+ [x ap x] where I is the identity matrix. t 0 is the reference epoch of the estimated coordinates x c,i and v c,i . In equation 4, the dependency of x c,i and v c,I on the epoch has been omitted but the appropriate coordinate parameters should be chosen in case of position and velocity discontinuities. The simultaneous estimation of transformation parameters together with station positions and velocities makes the system rank deficient (Altamimi et al., 2000) where the rank deficiency corresponds to the number of parameters that are necessary for the definition of the output reference frame. Constraints can be either added to parameters x c,i and v c,i -minimum constraints also called partial inner constraints involving coordinates -or to parameters t k , s k , θ k -internal constraints also called partial inner constraints involving transformation parameters.

In the current ITRF construction, both types of constraints are used, depending on the datum parameter to be defined. For example, internal constraints are used for translation and scale parameters and minimum constraints for orientation parameters (see [START_REF] Altamimi | ITRF2005: A new release of the International Terrestrial Reference Frame based on time series of station positions and Earth Orientation Parameters[END_REF]2011;2016).

II.3.3 Partial inner constraints involving coordinates

Minimum constraints as defined by Altamimi et al. (2001) are algebraic constraints that are imposed on coordinates to force them to have zero transformation parameters with respect to a reference coordinate dataset.

They can be derived from inner constraints that involve the a priori coordinates as reference coordinates of the constraints. The latter read as follows [START_REF] Dermanis | Global reference frames: connecting observation to theory and geodesy to geophysics Dermanis, A[END_REF]:

E T (x c -x ap ) = 0 (5) E T (v c -v ap ) = 0
Where E is the matrix of partial derivatives of coordinates with respect to transformation parameters t k , s k , θ k . x ap and v ap are the vectors of a priori coordinates. In the scope of Altamimi and Dermanis (2012b) developments, this constraint is a partial inner constraint involving only coordinates since the problem to be solved here also involves transformation parameters. However, given a constraint equation, any new constraint equation can be derived by multiplying it by an invertible matrix (Vanícek and [START_REF] Krakiwsky | Geodesy The Concepts[END_REF]. Minimum constraints as defined by Altamimi et al. (2001) are given as follows:

(E T E) -1 E T (x c -x ap ) = (E T E) -1 E T (x ex -x ap ) (6) (E T E) -1 E T (v c -v ap ) = (E T E) -1 E T (v ex -v ap )
where (x ex , v ex ) is an external coordinate set used as reference for the constraints.

II.3.4 Partial inner constraints involving transformation parameters

An alternative type of constraints involves the transformation parameters, in the case of stacking of station position time series. Since any input coordinate dataset at epoch k requires transformation parameters (t k , s k , θ k ) at the same epoch, time series of transformation parameters are actually estimated. If one wants to obtain a time series of transformation parameters with zero mean and zero drift, the following equation must hold for a parameter time series p:

(7)

This equation can be obtained by nullifying the least squares estimates of the constant and drift of a linear regression. By performing the computation, the first equation can be simplified as follows:

(8

)
This equation is the partial inner constraint involving only transformation parameters [START_REF] Altamimi | ITRF2005: A new release of the International Terrestrial Reference Frame based on time series of station positions and Earth Orientation Parameters[END_REF], Dermanis 2012b).

For a more physical interpretation of this type of inner constraints, we suggest working on a simplification of equation ( 4) where only translations would be involved:

The problem of solution numbers is easily handled since a new linear segment can be handled similarly by considering it as a new station. We provide an analytical solution of equation ( 9) under conditional least squares, considering K sets of coordinates as observations and additionally assuming that  t 0 is the average time over all t k ,  the condition equation ( 8) applies to the three components of the time series of vectors t k ,

 the input sets of coordinates have identity variance-covariance matrices.

It can be shown that the stacked coordinates are in this case the following [START_REF] Collilieux | Comparison of very long baseline interferometry, GPS, and satellite laser ranging height residuals from ITRF2005 using spectral and correlation methods[END_REF]):

(10) Thus, the estimated coordinates are those obtained from a standard linear regression without any translation parameters. The origin of the obtained frame has no drift and no mean offset with respect to the input frames due to equation 8. This result reinforces that adding constraints as in equation ( 8) allows defining a mean frame. Moreover, as proposed by [START_REF] Altamimi | ITRF2008: an improved solution of the International Terrestrial Reference Frame[END_REF]2016), it is possible to apply the constraint equation 8 on a subset of time indices t k corresponding to more reliable input coordinate sets, which allows to introduce all the other coordinate sets without affecting the frame definition of the combined coordinate set. In practice, those types of constraints are chosen to define the frame of the technique stacked solutions that contribute to the ITRF frame definition.

II.3.5 Combination

The second-step of the ITRF processing is the combination. The stacked coordinates from the different techniques are combined with local ties. In this process, transformation parameters between each coordinate dataset and the coordinates to be computed are estimated (Equation 1 of [START_REF] Altamimi | ITRF2005: A new release of the International Terrestrial Reference Frame based on time series of station positions and Earth Orientation Parameters[END_REF]. 14 constraint equations have to be added to define the frame of the output coordinates since positions at the reference epoch and velocities are linearly related to transformation parameters in the combination model. Because certain space geodetic techniques are sensitive to the origin and scale, defining the origin and scale of the combined coordinates consists in constraining some estimated transformation parameters to zero (Altamimi et al., 2001).

Algebraic constraints (see II.3.3) are added to x c (t) to define the orientation of the combined frame at the reference epoch and to v c (t) to define the orientation time evolution. Since ITRF2005, the conventional time evolution of the frame orientation has been defined through inner constraints with respect to the previously released velocity field. For ITRF2000, the geophysical model NNR-NUVEL1A, that verifies the no-net-rotation condition as specified in the ITRS system definition, has been used as reference velocity field for the orientation time evolution constraint. The choice of the orientation time evolution impacts studies that interpret polar motion drift. For users interested in that particular aspect, it is still possible to apply afterward a rotation contribution as done in [START_REF] Argus | An estimate of motion between the spin axis and the hotspots over the past century[END_REF] in order to change the orientation of the ITRF, or to use the kinematic methods suggested by [START_REF] Dermanis | Global reference frames: connecting observation to theory and geodesy to geophysics Dermanis, A[END_REF]2008).

II.4 Deficiency of the kinematic model in previous ITRF

There have been intensive researches on the analysis of position time series of space geodetic stations. For example, for GNSS, studies by [START_REF] Dong | Anatomy of apparent seasonal variations from GPS-derived site position time series[END_REF], [START_REF] Williams | Error analysis of continuous GPS position time series[END_REF] and [START_REF] Ray | Anomalous harmonics in the spectra of GPS position estimates[END_REF][START_REF] Ray | Subseasonal GNSS Positionning Errors[END_REF] are worth mentioning. In these studies, regularized positions and velocities are usually computed and removed from the determined coordinates to study their non-linear variations.

The following conclusions are worth reporting:

 Seasonal signals are evidenced in all technique estimated coordinates [START_REF] Collilieux | Comparison of very long baseline interferometry, GPS, and satellite laser ranging height residuals from ITRF2005 using spectral and correlation methods[END_REF][START_REF] Altamimi | Quality assessment of the IDS contribution to ITRF2008[END_REF].

 A significant part of seasonal signals is related to real crust motions, namely loading effects that are the elastic deformations of the Earth due to mass transports in its fluid layers (e.g. [START_REF] Mangiarotti | Annual vertical crustal motions predicted from surface mass redistribution and observed by space geodesy[END_REF]Mangiarotti et al., , van dam et al., 2001;;[START_REF] Dong | Anatomy of apparent seasonal variations from GPS-derived site position time series[END_REF][START_REF] Petrov | Study of the atmospheric pressure loading signal in very long baseline interferometry observations[END_REF][START_REF] Collilieux | Impact of loading effects on determination of the International Terrestrial Reference Frame[END_REF][START_REF] Fritsche | Global-scale validation of model-based load deformation of the Earth's crust from continental watermass and atmospheric pressure variations using GPS[END_REF], see section III.1.

 Loading effects, especially due to ice melting, can cause large non-stationary signals as found by [START_REF] Khan | Geodetic measurements of postglacial adjustments in Greenland[END_REF] in Greenland. See figure 2b) for an updated plot, where a seasonal signal is also visible.

 Non-tidal loading corrections are currently not recommended by the IERS (Petit and Luzum, 2010, Chapter 7).

 Technique specific errors are likely to remain in space geodetic technique coordinates [START_REF] Altamimi | Quality assessment of the IDS contribution to ITRF2008[END_REF][START_REF] Sarti | Height bias and scale effect induced by antenna gravitational deformations in geodetic VLBI data analysis[END_REF][START_REF] Ray | Subseasonal GNSS Positionning Errors[END_REF][START_REF] Appleby | Assessment of the accuracy of global geodetic satellite laser ranging observations and estimated impact on ITRF scale: estimation of systematic errors in LAGEOS observations 1993-2014[END_REF]  Co-seismic signals affect a large number of stations [START_REF] Tregoning | A decade of horizontal deformation from great earthquakes[END_REF][START_REF] Métivier | Global co-seismic deformations, GNSS time series analysis and earthquake scaling laws[END_REF].



III Station motion modeling in ITRF2014

III.1 Non-tidal loading effects

III.1.1 Background on non-tidal loading effects

The elastic deformation of the Earth in response to a surface load can be computed using either Green's functions or a spherical harmonic expansion [START_REF] Farrell | Deformation of the Earth by Surface Loads[END_REF]. In order to predict deformations at a point of coordinates Ω=(λ, φ), both the elastic properties of the Earth (described by Load Love numbers for a spherically symmetric, non-rotating, perfectly elastic and isotropic [START_REF] Dahlen | Theoretical global seismology[END_REF]) and the distribution of the load, i.e. its surface density σ(Ω), have to be known. If σ(Ω) is provided as a spherical harmonic expansion, it can be written as follows:

(11

)
where R n,m (Ω) and S n,m (Ω) are not normalized and related to associated Legendre polynomials by R n,m (Ω)=P nm (sin φ)cos(mλ) and S n,m (Ω)=P nm (sin φ)sin(mλ).The horizontal and vertical displacements can be computed by [START_REF] Blewitt | A New Global Mode of Earth Deformation: Seasonal Cycle Detected[END_REF][START_REF] Blewitt | Self-consistency in reference frames, geocenter definition, and surface loading of the solid Earth[END_REF]:

(12) (13)
with l' n and h' n the load Love numbers of degree n, R the mean radius, M the mass of the Earth and is the unit surface gradient operator (Blewitt and Clarke, 2003).

Thus, forward non-tidal loading models can be derived using numerical models of the mass distributions. Most of the time, forward models are generated by separating fluid layers [START_REF] Van Dam | Detection of atmospheric pressure loading using very long baseline interferometry measurements[END_REF][START_REF] Van Dam | Crustal displacements due to continental water loading[END_REF][START_REF] Van Dam | Nontidal ocean loading: amplitudes and potential effects in GPS height time series[END_REF] such as the atmosphere, ocean non-tidal mass transport, continental hydrology and ice sheet mass variations. However, care should be taken when summing models from various sources [START_REF] Clarke | Effect of gravitational consistency and mass conservation on seasonal surface mass loading models[END_REF] since interactions exist between the different fluid layers (exchange of water). The different models are available with various time samplings, from 3 hours to 1 month, and various spatial resolutions.

As discussed by [START_REF] Blewitt | Self-consistency in reference frames, geocenter definition, and surface loading of the solid Earth[END_REF], load Love numbers of degree-1 are dependent on the reference frame of the computed displacement field which are most often the center of mass of the solid Earth (CE), the geometric center of the solid Earth's surface (CF) or the center of mass of the whole Earth system (CM). CM is the natural frame for satellite techniques. The variations of the vector CM with respect to CF are usually called geocenter motion and the contribution of the gravitational attraction of the mass variations at the Earth's surface and associated elastic deformations can also be computed by using only degree-1 coefficients of σ(Ω):

(14)
The variations of CF with respect to CE are very small, about 2% of the geocenter motion magnitude for elastic deformations [START_REF] Dong | Geocenter variations caused by atmosphere, ocean and surface ground water[END_REF]. As stations are attached to the Earth's surface and because the analysis of translation parameters partial derivatives makes appear the barycenter of the station network [START_REF] Dong | Origin of the International Terrestrial Reference Frame[END_REF]Collilieux et al., 2009), net translations of a geodetic network with respect to mean coordinates are close to the opposite of geocenter motion. It is fundamental to understand that the geocenter motion is a part of the loading effects themselves and should be ideally treated simultaneously. Recent annual geocenter motion estimates predict translations with amplitudes at the level of 3-4 mm along the X, Y and Z components (Collilieux et al., 2009). Thus, mean coordinate variations in the CM frames are close to the geocenter motion magnitude [START_REF] Dong | Origin of the International Terrestrial Reference Frame[END_REF].

III.1.2 Handling of non-tidal loading effects

In the ITRF releases, non-tidal loading effects have been ignored since they were expected to be averaged out when computing positions and velocities from long enough datasets (Blewitt and Lavallée, 2002). However, for stations with a few observations, station positions and velocities have been reported to be biased by loading deformations [START_REF] Collilieux | Impact of loading effects on determination of the International Terrestrial Reference Frame[END_REF], affecting station coordinates from different techniques at co-location sites. This also means that if one wants to sum ITRF coordinates with a non-tidal loading displacement model, small biases may be generated since ITRF coordinates may have already captured a possible drift related to non-tidal loading effects.

There are mainly two ways to handle non-tidal loading effects according to equation (1):

 to correct for a non-tidal loading model, i.e. to model the non-tidal loading effect in the term ΣΔx(t),

 to add new parameters in the regularized coordinate kinematic model x r (t).

The first alternative is in theory the most rigorous. Indeed, it allows accounting for the effects over the whole power spectrum from sub-daily periods to long periods. When processing data for determining a reference frame, loading corrections should be ideally made at the observation level, which in addition allows modeling the induced effect on the geopotential for satellite techniques.

The inclusion of such models has been shown to reduce observation residuals although non-tidal loading models are not free of errors, for example evaluated at the level of 15% of the effect for nontidal atmospheric loading [START_REF] Petrov | Study of the atmospheric pressure loading signal in very long baseline interferometry observations[END_REF]. Hydrology models, for instance, show large discrepancies, particularly in the treatment of glaciers, ice sheets and lakes (e.g. [START_REF] Rodell | The global land data assimilation system[END_REF]. Due to the data processing time that can be important especially for GNSS and because they are quicker and more flexible, a posteriori corrections which consist in correcting coordinates by the mean load effect over the whole data integration time have been also studied (Tregoning and Watson, 2009). However, [START_REF] Dach | Evaluation of the impact of atmospheric pressure loading modeling on GNSS data analysis[END_REF] have suggested a relevant way to add geometric loading model corrections in the software packages while providing a rigorous way to remove them. They have added parameters to scale the non-tidal loading model corrections in the normal equations which allow applying the loading model or not by fixing them to zero or one when inverting the normal system.

The second alternative consists in estimating the loading displacement from station coordinates based on equation ( 12) and ( 13) as proposed by [START_REF] Blewitt | A New Global Mode of Earth Deformation: Seasonal Cycle Detected[END_REF]. [START_REF] Rülke | Realization of the Terrestrial Reference System by a reprocessed global GPS network[END_REF] were the first to estimate station positions and velocities simultaneously with load spherical harmonic coefficients. They estimated them on a monthly basis up to degree 6. The method allows modeling the geocenter motion effect with its associated degree-1 deformation with a relevant time sampling but the truncation degree used was too small to capture all non-tidal loading signals. Increasing the truncation degree would require to use additional datasets such as GRACE data [START_REF] Wu | Seasonal and interannual global surface mass variations from multisatellite geodetic data[END_REF], but the time span of the mission starting in 2002 does not cover the whole ITRF data period (see figure 1). [START_REF] Dermanis | The ITRF Beyond the "Linear'' Model[END_REF] suggestion of extending the station coordinate kinematic model to other function types is a relevant alternative since the main observed coordinate variations can be captured by parametric functions. Polynomials, splines and harmonic functions were suggested. Polynomial functions could for example capture the accelerations observed in figure 2 for station THU3 in Greenland. However, it is not certain that coordinates will follow the same pattern after the publication of the regularized coordinates. Indeed, the kinematic coordinate model is used to extrapolate coordinates in the future for routine precise orbit determination. Spline functions also fail to provide reliable extrapolations of coordinates although they provide an excellent fit to the data. Periodic functions are more relevant and were already used in reference frame analyses [START_REF] Meisel | Influence of Time Variable Effects in Station Positions on the Terrestrial Reference Frame[END_REF]. They have the advantage of providing a more reliable estimation of station velocities [START_REF] Blewitt | A New Global Mode of Earth Deformation: Seasonal Cycle Detected[END_REF]. In addition, as many position discontinuities are evidenced in particular in GNSS position time series, the addition of periodic parameters in the functional model improves the estimation of offsets. The main drawback is that periodic terms may also absorb technique-specific errors. In addition, they are not constant in time [START_REF] Chen | Singular spectrum analysis for modeling seasonal signals from GPS time series[END_REF], but their amplitudes vary moderately. Some cases of drastic seasonal signal changes in GPS coordinate series have been reported [START_REF] Ray | Systematic errors in GPS position estimates[END_REF] but they were related to equipment changes.

III.2 ITRF2014 coordinate modeling and estimation model

For ITRF2014, the ITRS product center explicitly asked not to correct for non-tidal models at the observation level. But tests have been carried out to correct for non-tidal atmospheric loading a posteriori [START_REF] Altamimi | ITRF2014: A new release of the International Terrestrial Reference Frame modeling nonlinear station motions[END_REF] during the ITRF preparation. However, the adoption of periodic parameters was preferred since they impact most significantly the estimated velocities. That's also why only annual and semi-annual periods have been considered, since other frequencies have a smaller impact on the estimated reference positions and velocities [START_REF] Altamimi | ITRF2014: A new release of the International Terrestrial Reference Frame modeling nonlinear station motions[END_REF]. However, seasonal parameters were not published as part of the coordinate kinematic model, see discussion in section III.3.

New functions have been used to model post-seismic deformations at stations where it was needed. The reader is referred to the appendix C of [START_REF] Altamimi | ITRF2014: A new release of the International Terrestrial Reference Frame modeling nonlinear station motions[END_REF] for a complete description of the post-seismic models that were adopted. A model, hereafter written δPSD(t), is made of a combination of exponential and logarithmic functions with various relaxation times. Not more than two functions are used per earthquake and coordinate component in a frame directed along the parallel, the meridian and the normal at the ellipsoid at the station location. For stations far from the earthquake epicenter, δPSD(t) is set to zero since the post-seismic motion can be represented by position offsets and constant velocity changes, as done in previous ITRF releases. Figure 2b) shows a station affected by significant post-seismic signal. The piece-wise linear function part from the kinematic model is shown in green and the complete model including the correction δPSD(t) in red.

As evidenced in this figure, users need to take the post-seismic deformation models δPSD(t) into account when using ITRF2014 positions and velocities.

Equation ( 4) of the stacking has then been modified as follows:

(15)

where two frequencies ω j /(2π) introduced (1 and 2 cycles per year). At a specific frequency, the signal is parameterized by sine and cosine amplitudes making the model linear. The term δPSD(t) has been estimated beforehand from GNSS input coordinate series, the best model being selected by a statistical method. Input series have then been corrected for that model which was assumed identical for co-located stations.

Even if the post-seismic models δPSD(t) were set up and the associated parameters were estimated in the stacking, the addition of these parameters would not cause an additional rank deficiency in the normal equation of the stacking. Indeed, post-seismic models only affect a part of the whole station network. Conversely, as periodic parameters are added for all stations, supplementary rank deficiencies are added to the normal equation for each frequency [START_REF] Petrov | Study of harmonic site position variations determined by very long baseline interferometry[END_REF] due to the simultaneous estimation of transformation parameters. For example, any periodic signal (a j ,b j ) added to all stations can be canceled by subtracting the same signal from the estimated translations. . As a consequence, condition equations must be added for every frequency. It is either possible to add constraints on the station periodic signals or on the transformation parameters.

As done for equation ( 7), it is possible to first cancel periodic signals at every frequency ω /(2π) in the transformation parameters, which gives for a transformation parameter p :

(16)

This condition can be used for canceling the annual signal in the scale parameter. Indeed, it is likely that most of the seasonal variations in the scale parameters be related to loading signals (Collilieux et al., 2010). If applied to the translation parameters, this constraint transfers any seasonal nettranslation of the network to station seasonal parameters.

The alternative approach consists in constraining the seasonal parameters a ] T , minimum constraints can be added as follows:

(E T E) -1 E T (α j -α j,ap ) = (E T E) -1 E T (α ex j -α j,ap ) (17) (E T E) -1 E T (β j -β j,ap ) = (E T E) -1 E T (β ex j -β j,ap )
where α ex j and β ex j have been fitted from an existing loading model or from GRACE observations. This algebraic constraint provides an alternative to the equations derived by [START_REF] Dermanis | The ITRF Beyond the "Linear'' Model[END_REF] to explicitly cancel the relative angular momentum of the network accounting for seasonal signals. In order to be relevant, the constraint of equation ( 17) should be applied with respect to a model which verifies a no-net-rotation condition. In the loading theory that has been briefly introduced in section III.1, a rotation is described by the toroidal term of degree-1. But no degree-1 toroidal deformation is possible when the external forcing is a load on the spherical surface [START_REF] Métivier | Mantle lateral variations and elastogravitational deformations -I. Numerical modelling[END_REF]. Thus, the no-net-rotation condition of the loading model is transferred to the estimated seasonal terms as currently done for the orientation linear time evolution of the reference frame.

III.3 Estimated seasonal parameters

The estimated annual and semi-annual periodic signals as the result of the stacking are dependent on the constraints that have been applied as described in the previous section. In ITRF2014, the following options have been used in the stacking of the four techniques:

 Equation ( 16) was applied to cancel annual and semi-annual signals in the translation and scale parameters.

 Equation ( 17) was applied to define the orientation of the seasonal parameters over well distributed sets of stations. No external loading model was used as reference for the constraints, i.e. α ex j and β ex j were set to zero.

The first constraint impacts the origin of seasonal periodic signals which is theoretically CM for satellite techniques. However, as already found in various studies, the origin of DORIS, SLR and GNSS frames differ at seasonal time-scale (see figure 9 from [START_REF] Rebischung | The IGS contribution to ITRF2014[END_REF][START_REF] Altamimi | ITRF2014: A new release of the International Terrestrial Reference Frame modeling nonlinear station motions[END_REF][START_REF] Ray | IERS Analysis Campaign to Investigate Motions of the Geocenter[END_REF]. In addition, for ITRF2004, geocenter parameters have been introduced in all input GNSS frames while GNSS station coordinates have been aligned to IGb08 frame [START_REF] Rebischung | The IGS contribution to ITRF2014[END_REF] which makes the estimated ITRF2014 GNSS seasonal parameters more consistent with CF frame. The obtained frame origin of VLBI seasonal signals is also arbitrary because the origins of the input VLBI coordinate datasets have been constrained to an a priori coordinate frame beforehand. It is thus likely that mean seasonal biases exist between technique-specific estimated seasonal signals.

The second constraint affects the estimated polar motion coordinates. Due to the limited number of stations and their distribution in space, the station displacements due to loading generate a small residual rotation. According to [START_REF] Collilieux | Strategies to mitigate aliasing of loading signals while estimating GPS frame parameters[END_REF], this residual rotation can be mitigated when a well-distributed network of station is used, which explains the constraint chosen to define the orientation of the seasonal signals in the ITRF2014 computation.

Because the seasonal signals of the different techniques are not provided in the same reference frame as discussed above and because they were not combined in ITRF2014, they have not been provided to enhance the ITRF kinematic coordinate model. The question of delivering a unique seasonal displacement for all stations at the same co-location site needs also to be addressed. Indeed, whereas we could in principle expect the averaged seasonal signal to be more reliable than any individual one at co-location sites, it cannot be ensured as long as the origin of their discrepancies is not understood. In addition, seasonal signals should be provided in the CM frame according to the ITRS specifications. It would require that the SLR origin of the seasonal parameters be transferred to other technique seasonal parameters. And this could be done reliably only if technique specific seasonal parameters are similar (in the sense of a 6-parameter similarity) over the set of common stations. As a conclusion, seasonal parameters at co-located stations should be carefully analyzed before attempting their combination.

  fig.2 a) GNSS position time series at THU3 (Greenland). b) GNSS position time series at SAMP (Indonesia). Source: http://itrf.ign.fr

  i j and b i j . By writing α j = [a 1 j ... a i j ... a n j ] T and β j = [b 1 j ... b i j ... b n j


Assistant Professor - B Lossouarn

J.-F De

Professor G Kerschen

Multimodal Damping of a Nonlinear Structure with a Passive Piezoelectric Network

Keywords: Nonlinear absorber, Multimodal damping, Piezoelectric coupling, Resonant network, Passive control

This work proposes the implementation of a multimodal and fully passive piezoelectric tuned vibration absorber that mitigates several resonances of a nonlinear structure. By extending a principle of similarity, an analogous electrical network is designed in order to reproduce the dynamics of the mechanical structure. Several electrical resonances are simultaneously tuned to the mechanical resonances, thus providing the equivalent of a multimodal vibration absorber from electromechanical coupling through an array of piezoelectric patches. Furthermore, the use of a nonlinear capacitor in the analogous network generates an autonomous adjustment of the electrical resonance frequencies when the structure reaches the nonlinear domain. The interest of this method is proved experimentally by mitigating vibration over a wide frequency range that covers the first three modes of a beam with cubic nonlinearity.

Introduction

Considering piezoelectric tuned vibration absorbers, mechanical nonlinearities generate a detuning that seriously affects the damping performance. This limitation of the classical resonant piezoelectric shunt [START_REF] Hagood | Damping of structural vibrations with piezoelectric materials and passive electrical networks[END_REF] was analyzed by Soltani at al. [START_REF] Soltani | The nonlinear piezoelectric tuned vibration absorber[END_REF], who extended results obtained by Habib et al. in the case of a purely mechanical absorber [START_REF] Habib | A principle of similarity for nonlinear vibration absorbers[END_REF]. In order to retrieve a correct tuning of the absorber, a solution consists in the introduction of an additional nonlinearity in the electrical circuit. The interest of such a nonlinear piezoelectric shunt has already been proven numericaly and experimentally for vibration mitigation of a single nonlinear resonance [START_REF] Soltani | The nonlinear piezoelectric tuned vibration absorber[END_REF][START_REF] Lossouarn | Passive realization of a nonlinear piezoelectric tuned vibration absorber with a saturable inductor[END_REF]. The objective is now to extend the method to multimodal structures by implementing a multi-resonant piezoelectric network [START_REF] Lossouarn | Multimodal vibration damping of a beam with a periodic array of piezoelectric patches connected to a passive electrical network[END_REF] that is able to damp several resonances of a nonlinear structure.

Nonlinear piezoelectric tuned vibration absorber

The experimental setup is based on a cantilever beam involving a cubic nonlinearity. As shown in Fig. 1, a thin lamina is clamped at the left end of the beam, which generates the hardening nonlinearity. The structure is covered with an array of piezoelectric patches in order to mitigate vibration with a piezoelectric tuned vibration absorber [START_REF] Hagood | Damping of structural vibrations with piezoelectric materials and passive electrical networks[END_REF]. A passive inductor is made by winding copper wire around a magnetic circuit in ferrite material and connected to the piezoelectric patches. A linear inductor allows reducing the vibration amplitude of the beam at low excitation but Fig. 2(a) shows that an increasing forcing amplitude leads to a complete detuning of the resonant shunt because of the mechanical nonlinearity. A solution to overcome this effect is to introduce in the absorber a nonlinearity similar to that of the primary structure [START_REF] Habib | A principle of similarity for nonlinear vibration absorbers[END_REF]. Practically, it can be implemented with a nonlinear capacitor [START_REF] Soltani | The nonlinear piezoelectric tuned vibration absorber[END_REF] or a nonlinear inductor [START_REF] Lossouarn | Passive realization of a nonlinear piezoelectric tuned vibration absorber with a saturable inductor[END_REF]. Both solutions provide equivalent adaptive tuning by maintaining two equal peaks in the frequency response function. Such a result represented in Fig. 2(b) has already been proven experimentally by using a physical inductor subjected to magnetic saturation [START_REF] Lossouarn | Passive realization of a nonlinear piezoelectric tuned vibration absorber with a saturable inductor[END_REF].

Multimodal damping with an analogous piezoelectric network

As the objective of the present work is to achieve vibration damping of several resonances, the array of piezoelectric patches is interconnected with electrical components in order to implement a multi-resonant electrical network. Actually, coupling a mechanical structure to its analogous electrical network generates electromechanical energy transfers over a broad frequency range. The electrical analogue of a beam [START_REF] Lossouarn | Multimodal vibration damping of a beam with a periodic array of piezoelectric patches connected to a passive electrical network[END_REF] is considered for the present application. It is a fourth order transmission line made of inductors and transformers, as shown in Fig. 3. Capacitors are not explicitly represented because we directly benefit from the piezoelectric capacitance of the patches that couple the beam to the electrical network. Coupling the mechanical structure to its analogous electrical network gives the frequency response in Fig. 4, which shows the simultaneous vibration mitigation of the first three modes of the beam at low forcing amplitudes. Yet, higher amplitudes still lead to the undesired behavior illustrated in Fig. 2(a) if linear components are used in the electrical network. Contrary to a previous experimental valitation based on a nonlinear inductor [START_REF] Lossouarn | Passive realization of a nonlinear piezoelectric tuned vibration absorber with a saturable inductor[END_REF], a nonlinear capacitor is incorporated in the piezoelectric electrical network to ensure an equivalence with the nonlinear stiffness at the end of the beam. The analogy is not only valid in the linear domain but also in the nonlinear domain thanks to the principle of similarity developed in previous studies on nonlinear vibration absorbers [START_REF] Habib | A principle of similarity for nonlinear vibration absorbers[END_REF][START_REF] Soltani | The nonlinear piezoelectric tuned vibration absorber[END_REF]. The damping performance are thus maintained, providing similar frequency response functions for higher forcing amplitudes.

Conclusions

This work aims to show that the principle of nonlinear similarity can be extended to multimodal damping with piezoelectric networks. To this end, an analogous electrical network is first assembled in order to reproduce the dynamics of the considered structure in the linear domain. A multimodal coupling offers vibration mitigation over a wide frequency range. For higher forcing amplitudes, it is shown experimentally that the mechanical nonlinearity induces a serious mistuning, especially for the first mode of the beam. A nonlinear capacitor is thus used in the electrical network in order to generate variable electrical resonances, which maintain a multimodal vibration reduction. This proves that broadband damping of a nonlinear structure can be achieved with a fully passive vibration absorber.

Fig. 1 Fig. 2

 12 Fig. 1 Cantilever beam with additional clamping through a thin lamina.

Fig. 3 Fig. 4

 34 Fig.[START_REF] Habib | A principle of similarity for nonlinear vibration absorbers[END_REF] Nonlinear beam coupled to its analogous electrical network through an array of piezoelectric patches.


Ibrahim Khalil

Hamzaoui

Mohammed Berrajaa
email: berrajaamo@yahoo.fr

Mostafa Azizi
email: azizi.mos@ump.ma

Giuseppe Lipari
email: giuseppe.lipari@univ-lille.fr

Pierre Boulet
email: pierre.boulet@univ-lille.fr

Measurement-based Methodology for Modeling the Energy Consumption of Mobile Devices 1 Measurement-based Methodology for Modeling the Energy Consumption of Mobile Devices

Keywords: Mobile computing, Operating system, Energy consumption modeling

Energy consumption is the result of interactions between hardware, software, users, and the application environment. Optimization of energy consumption has become crucial, and energy is considered a critical resource, so it is important to know and understand both how energy is measured and consumed on mobile devices. An accurate knowledge will allow us to develop efficient solutions to reduce energy consumption in order to improve the user experience. In this paper we propose an experimental methodology to build a model of the energy consumption of mobile applications. Based on precise measurements, we elaborate predictive models of energy consumption for both unconnected and connected applications.

Introduction

The energy consumption of a smartphone is the amount of energy used to operate its services [START_REF] Khalil | Survey on adaptation techniques of energy consumption within a smartphone[END_REF][START_REF] Song | Reducing energy consumption of smartphones using user-perceived response time analysis[END_REF] by the software running on the hardware. Thus, the hardware components used by the applications determine the resulting energy consumption. To know how to measure and to understand how the energy is consumed on these mobile devices remains an important objective of every developer [START_REF] Khalil | Estimation and optimization of energy consumption on smartphones[END_REF].

The main research presented in this paper consists of developing and practicing a new experimental methodology for modeling the energy consumption in mobile environments and more specifically on the Android platform. The objective of this work is to model the energy consumption of a particular application running on a mobile device.

In this paper, we propose a model and we describe a methodology to identify the parameters of this model (processor frequency, dissipated power and initial battery level). To this end, we analyzed a collection of experimental data collected during a 33 day long "Tour de France" in a wheelchair [START_REF] Khalil | An experimental methodology for modeling the energy consumption of mobile devices[END_REF]. The objectives of our work is to propose a methodology that a developer can use to understand the energy consumption behavior of a device in a specific applicative context. This methodology is based on measurements done on the device under study and leads to the proposal of a mathematical model with identified parameters of the energy consumption of the device. The resulting model can then be used to make decisions to balance energy efficiency and quality of service.

We present the methodology in Section 2, and we validate this methodology on experimental data in Section 3. Then we present related works in Section 4. We conclude our paper in Section 5.

Methodology

The research presented here proposes a methodology to model and evaluate the energy costs in mobile environments. The proposed methodology is developed in four stages:

• Data collection;

• Data preparation;

• Modeling for fixed frequencies;

• Modeling for variable frequency.

The aim of this work is to monitor energy consumption by acting on the following parameters:

• Frequency of processors;

• Initial level of the battery;

• Energy dissipated by clock cycle.

The model can be used to define an optimal frequency in terms of energy consumption for specific situations without degrading too much the quality of service desired by the user.

Data collection

Scientific data was collected on the basis of several scenarios in different places. These measures were operated out in a real environment and confronted with the consumption measures carried out beforehand in a controlled environment.

The tools used to develop our study are:

Trepn Profiler is a diagnostic tool for profiling performance and power consumption of Android applications. All tests of this experimentation were processed by version V6.2s. Trepn profiler provides information on system status, network status, graph performance, speed, processor frequency etc.

Cronoid is an automation tool that allows performing tasks on a regular basis (like cron). It also enables automatic task running when the status of the terminal changes. The version used for this work is Cronoid-3.5.1.

CPU Frequency is a tool that allows the user to change the CPU frequency setting to save energy or to achieve better performance.

The data collection stage of our methodology is based on the following steps:

• Preparation of the test platform (CPU frequency management based on the governor) in order to have the rights to fix the frequencies of the CPU with the CPU Frequency tool.

• The role of the Cronoid tool is to automate tasks in order to minimize interaction with the user.

Data preparation

After each experiment, the data obtained is stored in a file. Given the large number of parameters and data retrieved, the measurements are organized into two tables. Two examples of such tables are shown in Table 1 and Table 2.

For the implementation of our methodology, we selected the following parameters, obtained by Trepn Profiler:

• Total load per CPU;

• Memory usage;

• CPU frequency;

• Battery level;

• Battery power.

All these measures are regularly recorded in a table at a fixed measuring interval.

The methodology used to identify the parameters of the model is based on the determination of the type of correlation between the number of total operations Op and the dissipated elementary energy E under a previously set frequency F which will be expressed as E(F, Op).

The total number of operations, Op, is calculated by summing the active clock cycles of each processor per measurement interval. This active clock cycle number is the load of the CPU multiplied by the measurement interval duration and by the CPU frequency. The dissipated energy during a measurement interval is the product of the battery power and the measurement interval duration.

Modeling

The third step of the methodology is to compute a regression of the energy per measurement interval as a function of Op separately for some fixed (using the CPU Frequency tool) frequencies. More details on the modeling stages will be given in Section 3 on several case studies. In the first case study, a linear regression is enough to have a good model. Then we compute a regression of the coefficient of the fixed frequency models to obtain a model E(F, Op) of the energy consumption per measurement interval as a function of the frequency F and of Op, the number of operations during this measurement interval. In the case of offline mode study, a linear regression is enough to build a good model. The connected mode case is better managed by a polynomial regression.

In the case of variable frequencies, the same methodology can be applied on a variable frequency experiment by filtering the rows of the recorded measurement table by frequency and dealing separately with each frequency. The measuring intervals when the frequency changes can simply be discarded.

Experiments

Experimental context

Each service activated in a smartphone induces an energy consumption relative to an amount of dissipated energy. The purpose of this study is to provide a detailed monitoring of the main sources of energy consumption. In the first case study, we will use the GSM network, which will be the only one to be in active state (Only the GSM network will be active but unused by the application under test). The other communication networks (3G, 4G, Wi-Fi and the GPS) remain inactive throughout the various experiments. The second case study aims at evaluating and measuring the energy consumed in connected use cases, sources related to the data communications are activated completely (3G / 4G or Wi-Fi) or partially (GPS) depending on the scenarios of the experiment.

Each experiment was carried out over a period of 45 minutes with a rate of ten measurements per second, leading to around 27000 measurements per test case. For the Local Video with Fixed Frequency (LVFF) scenario, twenty experiments were carried out with different frequencies, initial battery levels and locations totaling 540000 measurements. We believe that 10 measures per second is a good compromise between the ability to capture DVFS induced frequency changes to get precise measures and the overhead induced by the measuring process. Indeed there are on the order of 100,000 processor cycles per measuring interval (at 1 GHz).

To collect the data, a tour of France was carried out using an electrical wheelchair to perform measurements in a diversified environment based on several scenarios. Figure 1 details the route of the measurement sites (33 days and 3006 km).

All measurements were performed on an HP Pro Slate 8 tablet that has the characteristics listed in Table 3.

Modeling of energy consumption

In this section, we present two cases. The first case corresponds to the disconnected mode, where we study local video playing with a fixed frequency while driving to test the impact of the variation of the GSM network on the energy consumption in disconnected mode.

Modeling in offline mode

The disconnected mode disables all data connection components such as Wi-Fi, Bluetooth and GPS. In this mode, we used local data for energy modeling.

First scenario, Local Video with Fixed Frequency (LVFF): in this scenario the CPU is used to monitor the system state by initiating and measuring the overall system power consumption for a local video with a fixed in advance frequency.

Second scenario, Local Video with Variable Frequency (LVVF): in this scenario the CPU is used to monitor the system state by initiating and measuring the overall system power consumption for a local video with variable frequency. This scenario is used for the validation of the proposed methodology trained with the scenario with fixed frequencies.

Modeling in connected mode

The second case corresponds to the connected mode. On the one hand, we study the case of a remote video with the fixed and variable frequency in motion to test the impact of the variation of the GSM network 3G / 4G or Wi-Fi on the energy consumption in connected mode. On the other hand we study the case of navigation where, in addition to activating the GSM network 3G / 4G or Wi-Fi and Bluetooth, we added the GPS for testing the impact of navigation on the energy consumption.

The active applications during the measurements are: Facebook, GsmaService, Google App, Clock, Avast Mobile Security, Cronoid, CPU Frequency, Google Play Store, Launcher3, Power Battery, Google services, Play, Hangouts, Multimedia storage, Google Play Music, Play-Fi, Keyboard, Google, MusicFX, SmartcardService, System interface, Messenger , VLC, Cron Tasker Free, and Google Partner Configuration.

In all the experiments of the video scenario, we worked on a video of 720 × 302 resolution with 25000 frames using the MPEG-4 video codec over a period of 45 min, the same experiment has been repeated with different processor frequencies (previously set) and with different initial battery levels.

The aim of this part is to evaluate the dissipated energy per clock cycle which will be called "dissipated elementary energy".

To qualify the quality of the model, we systematically show the value of the coefficient of determination which represents the square of the linear correlation coefficient R which is calculated as follows:

R 2 = n i=1 (ŷ i -ȳ) 2 n i=1 (y i -ȳ) 2 = 1 - n i=1 (y i -ŷi ) 2 n i=1 (y i -ȳ) 2 ∈ [0, 1]
y i are the values of the measurements, ŷi the predicted values and ȳ the mean of the measurements.

The average consumption per clock cycle can be expressed as follows:

E clock cycle = Regression( 4 i=1 Load i , E dissipated )
E clock cycle shows the total energy per clock cycle, Load i the load of CPU i and E dissipated corresponds to the dissipated energy measured experimentally.

Results of measurements

Disconnected mode case

The various measures relating to the energy dissipated were carried out in relation to the number of total operations at a fixed frequency. Figure 2 shows an example of experimental data obtained by fixing the frequency of the processor at 1.5 GHz.

The graphs corresponding to the other frequencies (1 GHz, 1.7 GHz and 2.2 GHz) are presented in the same way in Figure 3.

The respective equations of the regressions obtained are: The value of the determination coefficient is equal to 1 in the case of measurement with the frequencies 1.5 GHz, 1.7 GHz and 2.2 GHz, which proves the efficiency of the proposed methodology for the local video. In the case of a frequency fixed at 1 GHz, the coefficient of determination is equal to 0.830. The relevance of the linearity hypothesis was affected by the beginning of the degradation of the quality of the video (following the frequency imposed in the experimentation) which explains the obtained result. The low R 2 is due to the fact that the required computing power to decode the video is too high for a 1 GHz frequency, thus the QoS is degraded and it shows on the R 2 coefficient.

E(1.
On the basis of the experimental results obtained, we can deduce the energy discharge per clock cycle while playing a local video with a fixed frequency. The dynamic consumption corresponds to the energy consumption of the processor due to the fact that the transistors change state. The dynamic energy consumption per cycle for the measured frequencies is presented in Figure 4.

The static consumption corresponds to the energy dissipated when the transistors remain in the same state. The static energy consumption for the measured frequencies is presented in Figure 5.

The respective equations of dynamic energy and static energy dissipated for our model obtained by linear regression are:

E dynamic (F ) = 4.008.e -11 × F -2.101.e -11 , R 2 = 1.000 E static (F ) = 0.03 × F -0, 009, R 2 = 1.000
From these measurements we propose the following model

E(F, Op) = E dynamic (F ) × Op + E static (F )
Where E is the total energy, F the selected frequency, Op is the number of operations per measurement interval (in our case 100 ms).

In this study, the coefficients related to the LVFF scenario are thus defined by the following equation:

E(F, Op) = (4.008.e -11 × F -2.101.e -11 ) × Op + (0.03 × F -0, 009)

Validation of the mathematical model

In order to validate our model, we use an experiment that has not been used for the development of the model, we select the results relative to the dissipated energy (measured experimentally) and we compare them with the energy consumption predicted by the proposed LVFF model.

The experiment that is used for the test consists of 45,945 measurements with a variable frequency. Thus for each measurement interval, we use the recorded frequency as input to the model. The comparison of the predicted quantities by the above mathematical model and the measurements obtained experimentally make it possible to define the precision of the predicted model. After comparison, there is a strong coherency with a relatively low error rate.

The average relative error between our model and the experimental measurements is 1.767 % with a standard deviation of 3.652 %. Further more, 93.9 % of the measurements obtained have an error rate of less than 8 % and 64.73 % of the measurements obtained have an error rate of less than 5 %.

To illustrate this coherency, we have plotted the distribution of the cumulative sum of the difference between the measurements obtained experimentally and the values predicted by our model in order to evaluate the cumulative error rate in all the experimentation. The error rates are grouped by error classes (we have chosen 1000 classes) in Figure 6. The vertical axis was standardized according to the number of existing classes, 1 corresponds to the largest error class (the one with the highest number of members).

The results obtained correspond to a Gaussian distribution, the class which contains the greatest number of errors corresponds at the error interval 0.1 %, a value which consolidates the choice of parameters of model LVFF.

Connected mode case

The purpose of this Section is to study energy consumption while activating Wi-Fi. To stay in the same context, the remote video was identical to the one used in disconnected mode. We determine the type of correlation between the energy consumption and total operations, and we crosscheck with the measurements obtained in the offline mode. Figure 7 shows an example of experimentation by fixing the frequency of the processor at 1.5 GHz for a remote video with fixed frequency (RVFF).

Figures 8 and9 show respectively that the static energy and dynamic energy of the RVFF model vary linearly with respect to the frequency.

The respective equations of dynamic energy and static energy dissipated for our model obtained by linear regression are:

E Dynamic (F ) = 5.002.e -11 × F -2.104.e -11 , R 2 = 0.958 E Static (F ) = 0.03 × F + 0, 03, R 2 = 0.780
For the case of our study, the coefficients related to the RVFF scenario are thus defined by the following equation: E(F, Op) = (5.002.e -11 × F -2.104.e -11 ) × Op + (0.03 × F + 0, 03)

Comparative study between LVFF and RVFF scenarios

In this part we discuss the energy gap existing between the LVFF scenario and the RVFF scenario under similar conditions (with the same frequency). To compare the energy variation between the two scenarios, we compute the difference between the measurements obtained experimentally with the RVFF scenario and consumptions predicted by the LVFF model. This difference corresponds the communication cost and is shown in Figure 10. The energy consumption related to the communication is well modeled according to a quadratic regression which gives a better result, more relevance and more simplicity than a linear regression.

Navigation case

In the navigation case we activate both data connection and the GPS to monitor the energy consumption behavior.

Table 5 shows that the navigation uses a lot of CPU resources (19.14 %) and thus dominates the total energy consumption of the system.

The first case is to set the frequency and to track the total dissipated energy of the active applications separately at first, then the overall consumption of the Navigation with Fixed Frequency (NFF) scenario. Figure 11 shows an example of experimentation relating to the NFF scenario along with the regression. Figures 12 and13 show the linear regressions leading to the variable frequency NFF model.

Navigation with variable frequency

This part deals with the case of navigation with variable frequency (NVF), the frequency being assigned by the operating system. Figure 14 presents the results of NVF experiments. We note the existence of a data flow that strongly depends on the frequency at each measurement. The large number of operations and the rate of charge explain the shape of the obtained graph.

The respective equations of dynamic energy and static energy dissipated for our model obtained by linear regression are:

E Dynamic (F ) = 4.002.e -11 × F -2.103.e -11 , R 2 = 0.999 E Static (F ) = 0.098 × F -0, 0391, R 2 = 0.940
For our study, the coefficients related to the NFF scenario are thus defined by the following equation: E(F, Op) = (4.002.e -11 × F -2.103.e -11 ) × Op + (0.098 × F + 0, 0391)

Comparative study between NFF and NVF scenarios

Here we discuss the energy gap between the NFF and NVF scenarios. To compare the energy variation between the two scenarios, we compute the difference between the measurements obtained experimentally with the NVF scenario and the values predicted by the NFF model to determine the cost due to the variable frequency. This difference is shown in Figure 15.

The result obtained mathematically based on the difference on the respective equations of the NVF and NFF scenarios is very close to that obtained experimentally, the results obtained are as follows:

N F F : 0.05 + 5.e -11 x

(1)

N F V : 0.07 + 2.e -11 x + 5.e -21 x 2

N F V -N F F : 0.02 -3.e -11 x + 5.e -21 x 2

The static energy remains constant in the different experiments, the dynamic energy difference measured experimentally for the NFF and NFV scenarios was compared with the energy obtained by the proposed model.

The equation obtained by the model is expressed as:

0.03 + 2.e -11 x + 2.e -21 x 2 (4)
The value of the dynamic energy obtained by the model (Equation 4) is close to that obtained experimentally (Equation 3). The difference obtained is due to the variability of the frequency (dynamic energy) and presence of tail energy which is relative to the devices (Wi-Fi, GPS . . . ) which are not necessarily permanently active during the experiments.

In this study, we treated the different types of dependencies between energy dissipation and frequency of multiple scenarios in offline mode (LVFF) and the connected mode (RVFF, NFF, NVF). The difference between the scenarios can be explained on the one hand by the variability of the received signal strength for GPS and Wi-Fi (for connected mode) which forces the activation / deactivation of the services that will generate an additional consumption due to the tail energy. On the other hand the scenarios based on the connected mode require the activation of additional devices such as the GPS, which increases the rate of the inputs / outputs which are partially responsible for additional energy cost.

Comparison with the state of the art

Several research studies have also proposed CPU frequency-scaling techniques by taking into account the quality of user experience [START_REF] Li | Smartcap: user experience-oriented power adaptation for smartphone's application processor[END_REF][START_REF] Khalil | Survey on adaptation techniques of energy consumption within a smartphone[END_REF].

Related Modeling Works

TailEnder

Balasubramanian, Balasubramanian and Venkataramani designed TailEnder [START_REF] Balasubramanian | Energy consumption in mobile phones: a measurement study and implications for network applications[END_REF] for measuring the characteristics of the energy consumption of various smartphones. TailEnder aims to minimize energy consumption for applications that can tolerate delay such as e-mails. In particular, the study of 3G's energy consumption characteristics has important and nonintuitive implications for the design of efficient energy application. Previous studies such as [START_REF] Lee | Impact of inactivity timer on energy consumption in wcdma and cdma2000[END_REF][START_REF] Yeh | Comparative analysis of energy-saving techniques in 3gpp and 3gpp2 systems[END_REF][START_REF] Chuah | Impacts of inactivity timer values on umts system capacity[END_REF] have studied the impact of different energy-saving techniques in 3G networks using analytical models [START_REF] Khalil | Survey on adaptation techniques of energy consumption within a smartphone[END_REF].

Cinder

Roy, Rumble, Stutsman, Levis, Mazieres and Zeldovich show in [START_REF] Roy | Energy management in mobile devices with the cinder operating system[END_REF] that the Cinder method, that is designed for phones and mobile devices, allows users and applications to control and manage the limited resources of the device such as energy.

Cinder's mechanism [START_REF] Roy | Energy management in mobile devices with the cinder operating system[END_REF] consists of detecting the main energy-consuming tasks and then assigning this information to the guiding applications in order to allocate the energy resources to the highest priority tasks [START_REF] Zeng | Ecosystem: Managing energy as a first class operating system resource[END_REF][START_REF] Zeng | Currentcy: A unifying abstraction for expressing energy management policies[END_REF][START_REF] Fonseca | Quanto: Tracking energy in networked embedded systems[END_REF][START_REF] Flinn | Powerscope: A tool for profiling the energy usage of mobile applications[END_REF].

BLAST

Mercati, Hanumaiah and Bloch have shown in BLAST (Battery Lifetime-constrained Adaptation with Selected Target) [START_REF] Mercati | Blast: Battery lifetime-constrained adaptation with selected target in mobile devices[END_REF] that performance can be maximized if the user can select a target battery lifetime.

This framework is an application-aware power management framework for mobile devices which controls operating conditions in order to meet a predefined battery lifetime.

Eprof

Pathak, Hu and Zhang have shown in Eprof [START_REF] Pathak | Where is the energy spent inside my app?: fine grained energy accounting on smartphones with eprof[END_REF] that some of the components such as Wi-Fi, 3G GPS have a tail behavior [START_REF] Pathak | Fine-grained power modeling for smartphones using system call tracing[END_REF] in which a component can enter a state of high power and remain in that state of energy beyond the end of the routine.

This approach evaluates the energy consumed by the I/O and the energy of the tail [START_REF] Pathak | Fine-grained power modeling for smartphones using system call tracing[END_REF][START_REF] Balasubramanian | Energy consumption in mobile phones: a measurement study and implications for network applications[END_REF] as well as the energy consumed because of wake-up caused by the wake locks.

RAPL

Hahnel, Dobel, Volp and Hartig, developed the RAPL (Running Average Power Limit) approach that allows to compare applications of two sectors and to show that they have different energy consumption while offering similar services [START_REF] Hähnel | Measuring energy consumption for short code paths using rapl[END_REF].

RAPL is based on energy sensors available for measuring the energy consumption, and takes into account the energy consumption in the software components.

The aim of the RAPL methodology is to measure the energy consumed by the device and that consumed by its pilot.

AppScope

AppScope is a system that automatically evaluates the energy consumption of applications running on Android smartphones [START_REF] Yoon | Appscope: Application energy metering framework for android smartphone using kernel activity monitoring[END_REF]. Its design is based on monitoring the Android kernel at a microscopic level. The objective of the AppScope model is to measure the energy consumption of applications using power equipment models and usage statistics for each hardware component.

Synthesis

We have previously published a more detailed study on the state of art in [START_REF] Khalil | Survey on adaptation techniques of energy consumption within a smartphone[END_REF].

Table 5 compare these energy consumption methodologies and the one we propose in this article on different criteria. Here are the meanings of the columns of the Table :   Model quality quality of the energy model, this criterion shows the size of the work devoted to providing an energy model.

Hard Optim. taking into account hardware optimizations, this criterion shows the quality of the energy model and the hardware energy consumption.

Soft Optim. taking into account software optimizations, shows the quality of the energy model according to the implementation of software.

QoS impact on the quality of service (QoS), a management concept that consists of giving priorities to a few applications to keep a desired level of performance. This criterion details the effect on QoS during the reduction of energy consumption.

Security taking security into account. The security criterion role is to prevent malwares to reduce battery life.

Dev developer help, this criterion is used to reduce the energy consumption of a software during its design.

Eprof and our model HBALB (HBALB corresponds to the initial letters of the names of the authors of this work) achieve better results for quality of the energy model and performance optimization. Eprof makes it possible to analyze the state of the energy demands and to estimate the rate of tail energy of some components. HBALB allows to model the energy consumption so that it is optimized for better performances. TailEnder, Cinder, RAPL satisfy the "Hard opti." criterion. TailEnder performs a minimization of the energy consumption for the application which tolerate a small delay by quantifying the use of the energy, and Cinder by allocating resources to the guiding applications by visualizing the rate of limited resources. RAPL makes an energetic comparison of a similar service and shows that the consumption is not identical.

Appscope, BLAST and HBALB enable optimization of existing applications. Appscope evaluates automatically the power consumption of applications at a microscopic level. It allows to detect events relevant to the operation of a component for better optimization. BLAST maximize performance while letting the device battery to last at least for a certain required lifetime. Our model HBALB allows to define the optimal parameters for a better performance which justifies their relevance in the "Soft optim." criterion.

The quality of service is well managed by Cinder because it automatically manages the allocation of energy resources, the HBALB model does not have much impact on quality of service because user interaction is minimal due to automation of scenario tasks with the Cronoid tool.

Cinder prioritizes the security system that manages the emergency call in case of scarcity of energy resources.

Eprof satisfies the "Dev" criterion because it can help the developers by analyzing the different energy states of the device (characteristics, tail energy etc).

Our proposal provides precise models with respect to the active applications, which is not treated by the other case studies.

Advantages of the proposed methodology

• Our methodology is based on the most elementary entity which allows a high degree of precision for the evaluation of the results.

• The identification of the model parameters is relevant because it is based on measurements carried out in a real and diversified environment.

• The proposed methodology allows for optimization by tracking sources of energy consumption.

• The proposed methodology can be used for decision support for a desired performance.

Limits of the proposed methodology

• The precision of the models depends on the electrical characteristics of the device.

• The proposed methodology is not automated, it requires manual installation of the tools used to developing the model.

• The regression type is not automatically detected.

• The methodology currently lacks GPU and modem energy consumption modeling. It relates the global energy consumption of the whole device to only the activity of the CPU, and so lacks precision when the energy consumption is heavily influenced by other components.

Conclusion

We have presented in this paper a methodology to build a model of the energy consumption of applications on mobile devices. This methodology starts by recording accurate measures of the energy consumption during the use of selected applications. Using these measures, we build a model of the energy consumption as a function of the activity of the processing cores for fixed frequencies by (linear) regression. The last step in our methodology consists in elaborating a model of the energy consumption as a function of the operating frequency and the activity of the processor by a regression on the parameters of the fixed frequency models. Experiments related to this methodology have been carried out on different scenarios:

• Local video case;

• Remote video;

• Navigation.

We have demonstrated this methodology on an experimental test case over local video playing. As this test case does not involve communication, the prediction capabilities of the model are quite good.

For connected use cases, we notice that the prediction is less precise given the number of parameters involved and the diversity of situations encountered. It seems necessary in these cases to take into account more precisely the communication costs. The question on how to do that is left for future work.

The main quality of our methodology comparing with related works is the relevance of the recorded data that leads to an accurate model. The proposed solution can be used to define an optimal frequency for one or more applications in order to offer a better user experience with a reduced energy consumption.

As future work, it would be useful to extend this methodology to deal with other use cases that exercise other components of the complex mobile hardware platform. Checking the contribution of the GPU to the overall energy consumption for gaming use cases would be an interesting study. Such a methodology can help developers to understand the energy consumption behavior of their applications while optimizing the quality of service.

  0, Op) = 0.0512 + 5.995.e -11 Op, R 2 = 0.830 E(1.5, Op) = 0.0321 + 4.226.e -11 Op, R 2 = 1.000 E(1.7, Op) = 0.0337 + 5.073.e -11 Op, R 2 = 1.000 E(2.2, Op) = 0.0309 + 7.210.e -11 Op, R 2 = 1.000

Figure 1 Figure 2 FFigure 3

 123 Figure 1 Route of the scientific tour of France

Figure 5 Figure 6

 56 Figure 5 Static energy in function of the frequency for the LVFF scenario

Figure 7

 7 Figure 7 Energy dissipation for RVFF scenario, Frequency=1.5GHz

Figure 8

 8 Figure 8 Static energy in function of the frequency for the RVFF scenario

Figure 10

 10 Figure 10 Difference of energy consumption between the measures of the RVFF scenario and the predicted values of the LVFF model.

Figure 11 Figure 12

 1112 Figure 11 Energy consumption by total operations for the NFF scenario with a fixed 1.7 GHz frequency

Figure 14 Figure 15

 1415 Figure 14 Energy consumption by total operations for the NVF scenario

Table 3

 3 Characteristics of the experimental equipment

	Type of Product	Tablet
	Operating system	Android 4.4.4 (KitKat)
	Processor	QUALCOMM Snapdragon 800
	Max. CPU frequency	2.3 GHz
	Number of cores	4
	Sensors	Accelerometer, ambient light,
		proximity, compass, barometer,
		gyroscope, Hall effect
	Number of batteries	1
	Technology	Lithium polymer
	Autonomy	Up to 13.75 h
	Operating temp mini	0°C
	Operating temp maxi 40°C
	Storage temp mini	-20°C
	Storage temp maxi	60°C

Table 4

 4 Distribution of system resources for the NFF scenario

	Applications	CPU %
	Maps	19.14
	Parameters	3.93
	Messenger	0.75
	Interface	0.45
	Power Battery	0.4
	Cronoid	0.38
	Configuration of	0.24
	Google Play Music	0.23
	Facebook	0.09
	Gmail	0.072
	Cron Tasker Free	0.03
	Hangouts	0.022
	Drive	0.005
	Application Launcher 1.48E-4
	CPU Frequency	5.51E-5
	Firefox	5.38E-5

Table 5

 5 Comparison of energy consumption methodologies

		Model	Hard	Soft	QoS Security Dev.
		quality optim. optim.
	TailEnder	+	++	NA	+	NA	NA
	Cinder	+	++	NA	++	+++	NA
	Eprof	++	+	+	+	NA	++
	RAPL	+	++	NA	+	+	+
	AppScope	+	NA	+++	+	NA	+
	BLAST	+	NA	++	+	NA	NA
	HBALB	++	+	++	+	NA	+