Daniele Serra
email: daniele.serra@dm.unipi.it

Federica Spoto

Andrea Milani

A multi-arc approach for chaotic orbit determination problems

Keywords: chaos, orbit determination, multi-arc strategies, computability horizon, chirikov standard map

)) to the same dynamical system. These strategies entail the determination of an orbit when observations are grouped in separate observed arcs. For each arc a set of initial conditions is determined and, in the case of the constrained multi-arc method, all subsequent arcs are constrained to belong to the same trajectory. We show that the use of these techniques in place of the standard least squares method has significant advantages: not only can we perform accurate numerical calculations well beyond the computability horizon, in particular the constrained multi-arc strategy improves considerably the determination of the dynamical parameter.

Introduction

The classical techniques of orbit determination were introduced between the XVIII and the XIX century for the solution of practical problems of celestial mechanics, like the determination of a future position of a newly-discovered small celestial body so as to be able to observe it again. Nowadays, orbit determination methods are also essential for the Radio Science experiments of space missions: from measurements of the range and range-rate of a spacecraft orbiting a planet we can infer the trajectory of the orbiter as well as the value of some parameters characterizing the planet, e.g. the gravity field coefficients, the tidal parameters, the inclination of its rotation axis, etc. In abstract terms, orbit determination can be applied to reconstruct an orbit of a dynamical system, starting from a set of observations. The analytical expression of the dynamical system is supposed to be known, therefore in order to reconstruct the orbit it is sufficient to determine the initial state at an initial time instant, the so-called initial conditions. If the system depends on some other free parameters, the dynamical parameters, they can be recovered with orbit determination as well.

Many of the space-related applications of orbit determination are characterized by chaotic phenomena. For example, chaos appears when an asteroid has a deep close encounter with a planet (in the case of Near Earth Asteroids, especially the Earth) and its main effect is to change the orbit of the object and to increase the uncertainty of its state. Space missions grand tours of planetary systems represent another critical example of arising chaos. As in the case of the NASA/ESA/ASI Cassini-Huygens mission, the many close encounters of the orbiter with Saturn and its moons produced great uncertainty in the determination of the entire orbit, possibly affecting also the determination of some parameters regarding the planet.

Generally speaking, chaotic dynamical systems are characterized by a predictability horizon (cf. [START_REF] Lighthill | The recently recognized failure of predictability in Newtonian Dynamics [and discussion[END_REF]), a time limit beyond which it is impossible to forecast the state of the system within a desired accuracy. In the case of orbit determination of a chaotic system, this means that some parameters characterizing the system cannot be determined.

In an attempt of investigating the problem of the determination of chaotic orbits from a mathematically rigorous and abstract point of view, [START_REF] Spoto | Shadowing Lemma and chaotic orbit determination[END_REF] studied a relatively simple model problem in such a way that the results are representative of the real behaviour of more complicated systems, like the ones mentioned above. They considered the dynamical system known as Chirikov standard map (cf. [START_REF] Chirikov | A universal instability of many-dimensional oscillator systems[END_REF]), a discrete system of the two-dimensional torus which has both chaotic and ordered orbits 1 . They simulated a full-cycle orbit de-termination process: the observations were iterations of Chirikov standard map plus a Gaussian random error and used the non-linear least squares method, described in Sect. 2.1, to solve the problem of orbit determination of that dynamical system, determining a set of initial conditions and a dynamical parameter. Their goal was to understand the behaviour of the uncertainty in the determination of these parameters, as a function of the number of observations processed. For both chaotic and ordered orbits, they obtained that the uncertainties of the parameters decrease as a polynomial function of the number of observations. However, in the chaotic case they found that there exists a maximum number of observations that can be processed before the computation becomes numerically unstable. Such time limit is called computability horizon. It strongly depends on the chaoticity of the system, and it is a consequence of the numerical instability of the orbit computation. The computability horizon is therefore a realization of the predictability horizon when dealing with numerical experiments: its value, which depends on the Lyapunov exponent of the orbit, can be computed by a simple empirical formula. No computability horizon was found in the ordered case.

The first goal of this paper is to explore a different numerical technique to overcome the limitation coming from the computability horizon. To be able to compare our outcomes with the previous results, we consider the same problem model studied in [START_REF] Spoto | Shadowing Lemma and chaotic orbit determination[END_REF]. The novelty is the introduction of a multi-arc approach in the orbit determination. The idea of the classical pure multi-arc method (Sect. 2.2) is to consider observations grouped in separate and relatively short arcs, and determine initial conditions for each arc. If we also impose the constraint that all arcs belong to the same trajectory, we obtain the constrained multi-arc strategy (introduced by [START_REF] Alessi | Desaturation manoeuvres and precise orbit determination for the BepiColombo mission[END_REF] for the orbit determination of the ESA/JAXA BepiColombo mission to Mercury). We analyze the case when, in addition to the initial conditions, we also solve for a global dynamical parameter, shared by all arcs. The main advantage of multi-arc methods is that numerical computations are performed arc by arc, therefore we can choose the number of observations of each arc in such a way that the computability horizon is never met. We also aim to study the behaviour of the formal uncertainty associated with the solve-for parameters. We will be particularly interested in assessing whether the constrained multi-arc method provides an improvement in the determination of the dynamical parameter with respect to the pure multi-arc method.

This paper is organized as follows. In Sect. 2 we recall classical notions of orbit determination, multi-arc techniques and introduce the notation used throughout (cf. [START_REF] Celletti | Nearly-integrable dissipative systems and celestial mechanics[END_REF]). Moreover, it locally approximates the Kepler map, an area-preserving map describing the motion of a celestial body on a nearly-parabolic orbit (cf. [START_REF] Petrosky | Chaos and cometary clouds in the solar system[END_REF], [START_REF] Chirikov | Chaotic dynamics of comet Halley[END_REF] and [START_REF] Schevchenko | The Kepler map in the three-body problem[END_REF]).

the article. Sect. 3 describes our framework and presents the mathematical formulation of the pure and constrained multi-arc techniques in the case of Chirikov standard map. In Sect. 4 we set up the numerical experiment, show the results obtained, and discuss the differences between one chaotic example and one belonging to an ordered region. Finally, in Sect. 5 we draw our conclusions and hint at possible applications of this work.

Orbit determination

The scope of this section is to provide the essential notions of orbit determination, non-linear least-squares and multi-arc strategies, as well as to introduce the terminology used in the rest of the article. The reader may refer to [START_REF] Milani | Theory of Orbit Determination[END_REF] for a comprehensive discussion on the subject.

Non-linear least squares

Orbit determination indicates the process of reconstructing the orbit of a dynamical system (possibly dependent on some dynamical parameters µ)

ẋ = f (t, x, µ) ẋ(t 0 ) = x 0 ,
starting from a set of observations r 1 , . . . , r m at times t 1 , . . . , t m . Since the functional expression f on the right-hand side is usually known, the orbit x(t) is uniquely determined by x 0 and µ. Observations are not necessarily made up of direct measurements of the state of the system, but can be a function of it. For instance, in the case of a planetary orbiter the goal is to determine the planetocentric position and velocity of the spacecraft, but the observations are most commonly range (radial distance from the Earth) or range-rate (radial velocity w.r.t. the Earth) measurements.

To model the observations we use a prediction function r(t, x 0 , µ), depending on the time t, the unknown initial conditions x 0 and the dynamical parameters µ. Note that the distinction between initial conditions and dynamical parameters is purely formal: for this reason, in the following we will use the term "parameters" to refer to either x 0 or µ. The residuals are defined as

ξ i = r i -r(t i , x 0 , µ), i = 1, . . . , m,
thus the vector of the residuals ξ = (ξ i ) is a function of x 0 and µ. The least-squares method allows to compute a solution for a desired subset u of the parameters (the solve-for parameters), along with some statistical information associated with it.

The solution of the orbit determination problem should -in some sense -minimize the residuals. More rigorously, let us define the weighted target function

Q(u) = 1 m m i=1 ξ 2 i (u) σ 2 i ,
where σ -1 i is the weight2 associated to observation r i . The least-squares method states that the solution u is a local minimum point of Q. This leads to the differential correction algorithm, that is an iterative procedure to compute u: we start from an initial guess u 0 and generate a sequence (u j ), where u j+1 is obtained solving the normal equations:

C(u j+1 -u j ) = D, j ≥ 0 (1)
where

B = ∂ξ/∂u is called design matrix, C = B T W B is called normal matrix, D = -B T W ξ and W = diag(σ -2 1 , . . . , σ -2 m
) is the weight matrix. At convergence, the obtained value u * is the solution of the orbit determination problem. The statistical information about u * are contained in the covariance matrix Γ = C -1 : the formal uncertainty of parameter u i is √ Γ ii , and the correlation between two parameters u i and u j is Γ ij /( √ Γ ii Γ jj ). When some information on one or more solve-for parameters is available, it may be taken into account during the differential corrections process. Not only is this useful to stabilize the convergence of the algorithm, but it is also the only way to cure a possible rank deficiency of that specific orbit determination problem. Let u P be the apriori value of the N parameters that we want to include in the fit. With each apriori observation u i = u P,i , i = 1, . . . , N is associated an apriori standard deviation σ P,i , i = 1, . . . , N . If C P = diag(σ -2 P,1 , . . . , σ -2 P,N ), adding apriori observations u = u P can be realized modifying the target function Q as follows:

Q(u) = 1 m + N [ξ T W ξ + (u -u P ) T C P (u -u P )].
As a consequence, the normal equations considering the apriori observations are:

(C + C P )(u j+1 -u j ) = D + C P (u j -u P ).

Multi-arc strategies

In the case of planetary space missions with an orbiter, the spacecraft is never continuously tracked from Earth due to visibility conditions, mission design, body occultations of the probe, etc., therefore observations are naturally divided into arcs, separated by intervals in the dark. Dark periods are extended in time (days or even weeks) and scheduled manoeuvres often occur in that time. The uncertainty introduced by such manoeuvres and the difficulty of modeling the long-term effects of non-gravitational perturbations make the orbit determination process a challenging task. One possible solution is to use a classical (in the following, pure) multi-arc method: we can neglect all dynamics happening in the dark periods as long as we determine -in addition to the dynamical parameters -a set of initial conditions of the probe for each observed arc3 . For the practical calculation of the prediction function, this means that numerical propagation of the state of the spacecraft is made arc by arc. On the contrary, in a single-arc strategy the initial conditions have to be propagated for the entire duration of the mission.

Let us shortly describe the mathematical formulation of the pure multi-arc method. Let us suppose that the observations are divided in n disjoint subsets, called arcs, and let x k 0 be the initial conditions vector of arc k, k = 1, . . . , n. Let us divide the set of parameters in two subsets: the global parameters g, which are shared by all arcs, and the local parameters h 1 , . . . , h k , which are characteristic of the single arcs (e.g. h k includes x k 0 ). Moreover, the local parameters from one arc do not influence the dynamics of the other arcs. The residuals vector ξ is divided in k subvectors:

ξ(g, h 1 , . . . , h k ) = (ξ 1 (g, h 1 ), . . . , ξ k (g, h k ), . . . , ξ n (g, h n ))
and consequently ∂ξ i ∂h j = 0, for i = j.

The normal matrix inherits an arrow-like structure:

C =      C gg C gh 1 . . . C ghn C h 1 g C h 1 h 1 0 0 . . . 0 . . . 0 C hng 0 0 C hnhn     
, where

C gg = n k=1 (B k g ) T W k B k g C gh k = (B k g ) T W k B k h k = C T h k g , C h k h k = (B k h k ) T W k B k h k , k = 1, . . . , n,
W k is the weight matrix of arc k, and the design sub-matrices are

B k g = ∂ξ k ∂g , B k h k = ∂ξ k ∂h k , k = 1, . . . , n.
Finally, the right-hand side of the normal equations D is divided into sub-vectors:

D g = - n k=1 (B k g ) T W k ξ k D h k = -(B k h k ) T W k ξ k , k = 1, . . . , n.
This particular structure can be exploited to split the normal equations (1) in smaller systems, which are more convenient to solve from a numerical point of view. [START_REF] Alessi | Desaturation manoeuvres and precise orbit determination for the BepiColombo mission[END_REF] presented a new version of the multi-arc method for the Radio Science Experiment of the ESA/JAXA BepiColombo mission (cf. [START_REF] Benkhoff | BepiColombo-Comprehensive exploration of Mercury: Mission overview and science goals[END_REF]), called constrained multi-arc strategy. This is based on the idea that all arcs belong to the same orbit, therefore they should be connected in a smooth way. The connection is established as follows: for every pair of subsequent arcs, their initial conditions are propagated to the central time of the dark interval between them, and the difference between the two resulting states is constrained to be smaller than the observation noise (see later, Fig. 1). In the case of interplanetary missions, such strategy is more difficult to implement since the dynamical model of the spacecraft has to take into account the perturbations occurring in the dark periods. On the other hand, using such approach improves considerably the stability of the entire orbit determination process. The rigorous mathematical formulation of the constrained multi-arc method will be given in the next section for the model problem we are considering.

Multi-arc strategies for Chirikov standard map

In this section we describe the mathematical formulation of the pure and constrained multi-arc methods for the orbit determination problem of Chirikov standard map.

Let us start recalling the definition of the dynamical system we are dealing with. Chirikov standard map is a conservative discrete dynamical system, defined on a two-dimensional torus by the following formula:

S i+1 µ 0 [(x 0 , y 0 )] = x i+1 = x i + y i+1 y i+1 = y i -µ 0 sin(x i ), i ∈ N, (2)
where S 0 µ 0 [(x 0 , y 0 )] =: (x 0 , y 0 ) are the initial conditions and µ 0 is a parameter. The system could have both ordered and chaotic orbits, depending on the choice of (x 0 , y 0 ) and µ 0 . There are more regular orbits for smaller values of µ 0 , whereas larger values of this parameter result in more chaotic orbits.

Pure multi-arc method

The first element we need for setting up an orbit determination process are the observations. Since we are aiming at using a multi-arc method, we need to reproduce the same natural observational framework where multi-arc techniques are used in practice, thus we need observations organized in arcs and dark periods between them. As done in [START_REF] Spoto | Shadowing Lemma and chaotic orbit determination[END_REF], the observations are simply iterations of the dynamical system and are obtained by means of the following procedure. Firstly, we compute forward and backward iterations of (2) starting from the initial conditions (x 0 , y 0 ). Secondly, we add a random noise, modeled as a Gaussian distribution N (0, σ), obtaining (x i , ȳi ) = (N (x i , σ), N (y i , σ)). From this sequence of modified iterates, we extract observations grouped in 2n + 1 observed arcs, numbered from -n to n, in such a way that:

• each arc k contains 2m + 1 observations (x j , ȳj ) k , j = -m, . . . , m, and (x 0 , ȳ0 ) k denotes the central point of each arc;

• (x 0 , ȳ0 ) := (x 0 , ȳ0 ) 0 is the central point of arc 0;

• any two subsequent arcs are separated by a dark interval composed of 2m d +1 iterates, numbered from -m d to m d .

Consequently, the total number of observations is (2n + 1)(2m + 1). Since the observations are practically iterations of the map with a noise, in the following the word "iteration" will be used as a synonym of "observation". The prediction function is given by the standard map itself, depending on the dynamical parameter µ 0 and the central point of each arc (x 0 , y 0 ) k , which we will consider from now on initial conditions of the arc. The set of the solve-for parameters in the orbit determination process is {µ 0 , (x 0 , y 0 ) k , k = -n, . . . , n}. As already stated before, the parameters can be divided into global and local parameters. In our case, the vector of the global parameters g is only the dynamical parameter µ 0 , whereas the vector of the local parameters h is the vector of all the initial conditions {(x 0 , y 0 ) k , k = -n, . . . , n}, which can be split in subvectors

h k = (x 0 , y 0 ) k , k = -n, . . . , n.
The core of the least squares fit is the minimization of the target function 5) is the difference between the backward-propagated current initial conditions of arc k + 1, namely (x, y) k+1 , and the forwardpropagated current initial conditions of arc k, namely (x, y) k , to the central point of the dark interval between the two arcs.

Q = 1 (2n + 1)(2m + 1) n k=-n (ξ k ) T W k ξ k , (3) • (x, y) k • (x, y) k+1 S m+m d µ [(x, y) k ] • • S -(m+m d ) µ [(x, y) k+1 ] d k
where W k is the weight matrix of arc k and ξ k = (ξ k -m , . . . , ξ k m ) are the residuals relative to arc k, which read

ξ k j = (x j , ȳj ) k -S j µ [(x, y) k ], j = -m, . . . , m.
Here (x, y) k = (x 0 + δx, y 0 + δy) k and µ = µ 0 + δµ are the current guess of the parameters. The crucial property characterizing the multi-arc strategies is described by

∂ξ k 1 ∂(x 0 , y 0 ) k 2 = 0, for all k 1 = k 2 , k 1 , k 2 = -n, . . . , n, (4)
which expresses the fact that the residuals from one arc do not depend on the local parameters of another arc. We do not report the expressions of the normal equations and the matrices involved as they are formally identical to those in Sect. 2.2, except for summations, now running over k = -n, . . . , n. Note that these matrices have explicit formulae and they can be computed as shown in [START_REF] Spoto | Shadowing Lemma and chaotic orbit determination[END_REF], Sect. 2.

Constrained multi-arc

Let us now show the mathematical formulation of the constrained multi-arc method for Chirikov standard map. The observations, the prediction function and the residuals are the same as in the classical multi-arc strategy. The main idea is to impose the constraint that two subsequent arcs belong to the same orbit.

We define the jump vector between arcs k and k + 1 as (see Fig. 1)

d k := S -(m+m d ) µ [(x, y) k+1 ] -S m+m d µ [(x, y) k ], k = -n, . . . , n -1 ( 5
)
where S is given by (2), (x, y) k and (x, y) k+1 are the current guess of the initial conditions of arc k and k + 1 respectively, and µ is the current guess of the dynamical parameter. Constraining two consecutive arcs to belong to the same orbit is equivalent to including the following apriori observations to the fit:

d k = 0, (6)
with associated apriori standard deviation σ P . This corresponds to adding 2n quadratic forms to the target function (3), which now reads

Q = 1 (2n + 1)(2m + 1) + 4n n k=-n (ξ k ) T W k ξ k + 1 σ 2 P n-1 k=-n (d k ) T d k .
The constrained multi-arc design matrix contains also the derivative of the jump vectors with respect to the parameters:

∂d k ∂g =      ∂S -(m+m d ) µ [(x, y) k+1 ] ∂µ 0 - ∂S m+m d µ [(x, y) k ] ∂µ 0 k = -n, . . . , n -1 0 k = n ∂d k-1 ∂h k =      ∂S -(m+m d ) µ [(x, y) k ] ∂(x 0 , y 0 ) k k = -n + 1, . . . , n 0 k = -n ∂d k ∂h k =    - ∂S m+m d µ [(x, y) k ] ∂(x 0 , y 0 ) k k = -n, . . . , n -1 0 k = n.
The constrained multi-arc normal matrix C results in a block tridiagonal arrowshaped matrix:

C =         C gg C gh -n . . . . . . C ghn C h -n g C h -n h -n C h -n h -n+1 0 0 . . . C h -n+1 h -n . . . . . . . . . . . . 0 . . . . . . C h n-1 hn C hng 0 0 C hnh n-1 C hnhn         ;
where

C gg = C gg + 1 σ 2 P n-1 k=-n ∂d k ∂g T ∂d k ∂g C gh k = C gh k + 1 σ 2 P ∂d k ∂g T ∂d k ∂h k + 1 σ 2 P ∂d k-1 ∂g T ∂d k-1 ∂h k = C T h k g , k = -n, . . . , n C h k h k = C h k h k + 1 σ 2 P ∂d k ∂h k T ∂d k ∂h k + 1 σ 2 P ∂d k-1 ∂h k T ∂d k-1 ∂h k , k = -n, . . . , n C h k-1 h k = 1 σ 2 P ∂d k-1 ∂h k-1 T ∂d k-1 ∂h k = C T h k+1 h k , k = -n + 1, . . . , n.
Finally, the constrained multi-arc right-hand side is

D g = D g - 1 σ 2 P n-1 k=-n ∂d k ∂g d k D h k = D h k - 1 σ 2 P ∂d k ∂h k T d k - 1 σ 2 P ∂d k-1 ∂h k T d k-1 , k = -n, . . . , n.

Numerical experiments

Spoto [START_REF] Spoto | Shadowing Lemma and chaotic orbit determination[END_REF] showed that the determination of a Chirikov standard map's orbit with initial conditions (x 0 , y 0 ) = (3, 0), using the classical singlearc method and determining (x 0 , y 0 ) and the dynamical parameter µ 0 , is limited by the computability horizon. For that orbit and in double precision, after ∼ 200 iterations of the map (i.e, when processing more than ∼ 200 observations), computation becomes numerically unstable. Since the Lyapunov time4 of that orbit is about 11 iterations, the limit in this case corresponds to about 200/11 ∼ 18 Lyapunov times. Using a multi-arc strategy is a way to overcome this limitation: since such technique entails that the computation of the prediction function (which also requires numerical propagations of the dynamical system) is performed arc by arc, if the number of observations belonging to each arc is small enough, the computability horizon is never reached and numerical instability is thus avoided. In our experiments with the same dynamical system and the same chaotic orbit, we push orbit determination up to more than 60 Lyapunov times in double precision, as described in the following.

Setup of the experiment

We create the observations as explained in Sec. 3. We use σ = 10 -8 as the Gaussian noise standard deviation, and we choose m = 5, m d = 1, n = 50. In this way each arc is composed of 11 observations -approximately 1 Lyapunov time -, the dark intervals are composed of 3 iterations, and the total number of arcs is 101. The number of iterations required for simulating the observables is ∼ 700, forward and backward, which correspond to more than 60 Lyapunov times of the orbit with initial conditions (3, 0). With this choice of m and m d we are able to avoid the computability horizon.

We apply the least squares fit and the differential corrections algorithm via pure and constrained multi-arc approaches. As already described in Sect. 2, the differential corrections algorithm consists of solving iteratively the normal equations until a given convergence criterion is satisfied. We will discuss such criterion in Sec. 4.2. The solve-for parameters are the initial conditions (x 0 , y 0 ) k , k = -n, . . . , n of the arcs and the dynamical parameter µ 0 .

In order to ensure stability of the process, we tackle orbit determination in a progressive fashion. At step 0 we process only arc 0, solving for the initial conditions of this arc and µ 0 . For each parameter we use as first guess its true value (i.e. the value used to create the observables) incremented by 10σ. We do this to simulate a real orbit determination case, where the first guess of the parameters is usually very different from the value at convergence. Once convergence is reached, we have corrected values for such parameters. We use them in step 1 as new first guess in an orbit determination process in which we now process, along with arc 0, also arcs 1 and -1, and we solve for the initial conditions of the three arcs and µ 0 . The first guess for the initial conditions of the newly added arcs is given by their simulation value plus 10σ, as in the previous steps. In general, at step k + 1 we add arcs -k -1 and k + 1 to the set of the previously processed 2k + 1 arcs (from -k to k) and solve for the initial conditions of all arcs and µ 0 . The first guess values of the parameters are chosen as done in the previous steps. We repeat such scheme until the total number of arcs is reached.

Convergence conditions

The differential corrections entail the iterative solution of the normal equations, therefore it is crucial to establish a criterion for the termination of the algorithm. The principle of convergence is based on two conditions: on the one hand, we want to make sure that the last correction of the parameters (i.e. the difference between the value of the parameters at the present iteration and at the previous one) is "small", in a sense we are going to specify; on the other hand, when coping with a constrained multi-arc strategy, we need to ensure that the jump vectors are not larger than the observation noise, say at least ten times smaller. These two specifications translate into the following mathematical conditions: 1) if ∆u is the last correction to the solve-for parameters, we require that its norm ∆u C (C is the normal matrix of the adopted multi-arc strategy) satisfies

∆u C := ∆u T C∆u (2n + 1)(2m + 1) 1,
meaning that all following iterations will not improve significantly the solution; 2) in order to have a smooth orbit, the norm of the jump vectors should be negligible with respect to the observation noise σ, resulting in

d RMS := n-1 k=-n d k 2 4n ≤ σ * , ( 7
)
where d k is the euclidean norm of d k and σ * is a fraction of σ. The same value σ * is considered for the apriori standard deviation σ P associated with the contraints. We adopt a value of σ P which varies along with the iteration of the fit:

σ P = max d RMS 100 , σ * .

Results: chaotic case

The results presented in this section refer to the determination of an orbit of Chirikov standard map simulated with initial conditions (x 0 , y 0 ) = (3, 0) and µ 0 = 0.5, with the setup described in Sect. 4.1. Such an orbit belongs to a chaotic region of the phase space of the considered dynamical system. We have already observed that the use of a multi-arc strategy allows to avoid the computability horizon, therefore in the following we will focus on the uncertainties attainable for the solve-for parameters.

Pure multi-arc

Let us first present the results obtained when we apply the pure multi-arc method. Figure 2 shows the cumulative formal uncertainties on the initial conditions of arc 0 (x 0 , y 0 ) and the dynamical parameter µ 0 as a function of the total number of arcs k processed, k = 1, . . . , 2n + 1, in a log-log plot.

The uncertainty on µ 0 decreases as a polynomial function k a of the number of arcs, with a ∼ -0.5. This is evident if we observe that by definition of pure multi-arc strategy the observed arcs are independent, that is the information on the dynamical parameter contained in each arc is the same in all arcs. As a consequence, as we add more arcs to the fit, the normal matrix C grows as k, and the covariance decreases as k -0.5 , as we found. Figure 2: Formal uncertainties of the initial conditions of arc 0 (x 0 , y 0 ) = (3, 0) and the dynamical parameter µ 0 = 0.5 as a function of the total number of processed arcs, using a pure multi-arc strategy, chaotic case (log-log plot). The slope of the µ 0 formal uncertainty curve is ∼ -0.5.

The uncertainty on the initial conditions has a different behaviour. Component x 0 first decreases polynomially (roughly with the same exponent as µ 0 ), then it reaches a point where no further improvement is gained, therefore remaining almost constant. Component y 0 shows no improvement whatsoever. A similar behaviour was expected since by definition in a pure multi-arc strategy the information carried by each observed arc is relevant only for the global parameters and the local parameters of that arc. In particular, it does not affect the local parameters of other arcs. The moderate initial decreasing trend of the x 0 component is due to correlation with µ 0 : improving the knowledge of the latter helps better determine the former. After a number of arcs the correlation is broken and no further improvement for x 0 is obtainable.

Constrained multi-arc

Figure 3 shows the results obtained using the constrained multi-arc method in four different scenarios, corresponding to four different choices of the convergence parameter σ * , namely σ * = σ/10 h , for h = 1, 2, 3, 4 respectively. The smaller σ * , the tighter the constraints between the arcs.

As far as the uncertainties of the initial conditions are concerned, the constrained multi-arc does not change the qualitative behaviour observed in the pure multi-arc case. Both uncertainties decrease at first, then they tend to an asymptotic value. As σ * decreases, we observe that the asymptotic value is smaller, thus the accuracies attainable for the initial conditions are higher. This is in agreement with how the constrained multi-arc works: smaller values of σ * correspond to tighter constraints between subsequent arcs, thus better accuracy in the determination of the orbit. Of course σ * cannot be chosen infinitely small, or we would encounter numerical instability. We found indeed that for σ * < σ/10 4 the iterative procedure of differential corrections fails. It is worth remarking that for large values of σ * we find similar results to those obtained in the pure multi-arc case, since imposing a loose constraint between two subsequent arcs is in practice equivalent to not constraining them at all.

As regards the dynamical parameter µ 0 , in all cases its formal uncertainty is a polynomial function of the number of processed arcs, k a , where a < -0.5. The smaller σ * is, the smaller a: in other words, tighter constraints between arcs result in better determination of the dynamical parameter. Note that a is smaller than in the pure multi-arc case for every choice of σ * , therefore the determination of µ 0 benefits from the adoption of a constrained multi-arc technique.

Numerical experiments conducted on other chaotic orbits of the same dynamical system show that such behaviour is in fact general. Although the value a of the slope changes as we change the test orbit, we always find that a < -0.5. This result shows that chaotic orbits provide favourable opportunities for estimating free parameters of a dynamical system, and the only way not to deal with numerical computation issues is the use of a constrained multi-arc method.

Results: ordered case

In this section we investigate the results of orbit determination of Chirikov standard map using a multi-arc strategy (pure and constrained) in an ordered case. We chose (x 0 , y 0 ) = (2, 2), µ 0 = 0.5.

The results obtained using a pure multi-arc strategy are qualitatively similar to the ones obtained with the same strategy in the chaotic case (cf. Fig. 4). The uncertainty in the determination of µ 0 decreases indeed as k a , a ∼ -0.5, where k is the number of arcs, and the initial conditions of arc zero (x 0 , y 0 ) do not improve much, for the same reasons explained in Sect. 4.3.1.

Let us go on with the constrained multi-arc strategy. For comparison with the chaotic case, we examined the ordered case similarly to what was done in Sect. 4.3.2: we repeated the same experiment four times, each time with a smaller value of the convergence parameter σ * . As shown in Fig. 5, in all cases the formal uncertainty of the dynamical parameter µ 0 is a polynomial function of the number of arcs, k a , a ∼ -0.5. With regard to the initial conditions, their formal uncertainties strongly depend on the choice of σ * . For larger values of σ * they initially tend to decrease, then remain constant, analogously to the chaotic case. The reason is the same as in the previous case: when constraints are rather loose, adding observed arcs to the fit provides little information on other arcs, so the uncertainties 2, 2) and the dynamical parameter µ 0 = 0.5 as a function of the total number of processed arcs, using a pure multi-arc strategy, ordered case (log-log plot). The slope of the µ 0 uncertainty is ∼ -0.5.

of the initial conditions do not decrease after they reach a certain minimum value.

On the other hand, if σ * is small enough, the formal uncertainties of the initial condition follow the same law as the dynamical parameter (Fig. 5(d)). This is what happens if we attempt a single-arc solution (cf. [START_REF] Spoto | Shadowing Lemma and chaotic orbit determination[END_REF], Fig. 10), since in the limit for σ * → 0 the constrained multi-arc strategy tends to the single-arc one. Note that in contrast with the chaotic case, in the ordered case no instability is encountered when choosing σ * < σ/10 4 . We have studied indeed the case when σ * = σ/10 5 , one order of magnitude less than what is possible to do in the chaotic case. As a final remark, we observe that the determination of µ 0 in the ordered case is not affected by the choice of the particular multi-arc strategy used in orbit determination process. As in either case the formal uncertainty of the dynamical parameter is decreasing like k -0.5 , the constrained multi-arc method does not improve the determination of the parameters when dealing with non-chaotic problems. Nevertheless, it should be adopted in all practical cases when the stability of the computation is compromised by dynamical effects that are difficult to model, for instance in the BepiColombo MORE experiment [START_REF] Alessi | Desaturation manoeuvres and precise orbit determination for the BepiColombo mission[END_REF]). Figure 5: Formal uncertainties of the initial conditions of arc 0 (x 0 , y 0 ) and the dynamical parameter µ 0 as a function of the total number of processed arcs, using a constrained multi-arc strategy (log-log plot), ordered case. The four figures refer to four different experiments repeated with the same setting, (x 0 , y 0 ) = (2, 2), µ 0 = 0.5, except for different values of the convergence parameter σ * (see formula 7).

Conclusions

We have applied the classical multi-arc method and the recently-introduced constrained multi-arc method for orbit determination of Chirikov standard map. We simulated observations organized in disjoint arcs and we have solved for a pair of initial conditions for each arc and a global dynamical parameter. We have compared the results obtained with the analysis carried out by [START_REF] Spoto | Shadowing Lemma and chaotic orbit determination[END_REF], where they considered a unique observed arc and solved for a single pair of initial conditions and the same dynamical parameter. Since in the multi-arc strategies the computation of the prediction function is performed arc by arc, by choosing relatively short arcs (about one Lyapunov time of the orbit) we have avoided the computability horizon found in the single-arc strategy when trying to compute a prediction for more than 18 Lyapunov times. We have indeed shown that we are able to compute an orbit of the system for at least 60 Lyapunov times, and the limit is far from being reached. We went even beyond the computability horizon encountered in quadrupole precision, occurring after ∼ 39.2 Lyapunov times (cf. [START_REF] Spoto | Shadowing Lemma and chaotic orbit determination[END_REF]), which confirms the remarkable strength of this method. Since the uncertainties of the solve-for parameters decrease as the number of observations increases, this allows to achieve higher accuracies in the orbit determination of a chaotic system.

For the determination of a chaotic orbit, in the case of a pure multi-arc strategy, we found that the knowledge of the initial conditions cannot be significantly improved increasing the number of observed arcs. Although we have shown such result for the initial conditions (x 0 , y 0 ) 0 of the central arc, it is also valid for the initial conditions of all the other arcs. On the other hand, the constrained multiarc strategy helps to improve the accuracy on the initial conditions as we add more arcs, provided that we constrain all subsequent arcs to belong to the same trajectory with very small apriori uncertainty σ * (as large as a fraction of the observation noise). Moreover, we observed that the smaller σ * is, the larger the improvement we obtain. However, we also noted that we are not allowed to choose σ * arbitrarily small. This choise would correspond to impose that the arcs are connected with infinite precision; were we allowed to do so, we would also be able to use a single-arc method, as in [START_REF] Spoto | Shadowing Lemma and chaotic orbit determination[END_REF], to process as many observations as desired, which is not the case, as it has already been stated. In fact, there exists a value of σ * , such that the differential corrections algorithm cannot find a solution. This is the counterpart of the computability horizon for the constrained multi-arc: we can process a number of observations much larger than in the single-arc case, but we are only able to reconstruct a piecewise orbit, whose components are connected within a certain approximation. As regards the determination of the dynamical parameter, we found that its formal uncertainty decreases as a polynomial function of the number of arcs k. Such functional law is k -0.5 when applying a pure multi-arc technique and k a , a < -0.5 for the constrained one, a being smaller when tighter constraints between the arcs are applied. This implies that, using the same number of observations, the constrained multiarc method helps to better determine the dynamical parameter with respect to the pure one.

We have also studied an ordered case. We noted that as opposed to the chaotic case we can push orbit determination with constrained multi-arc method to smaller values of the apriori uncertainty σ * . However, the resulting uncertainty on the dynamical parameter does not depend on the multi-arc technique used: in both cases it decreases as k -0.5 . On the other hand, we observed that using a constrained multi-arc method improves the determination of the initial conditions of arc 0, provided that we use a small value of σ * . For σ * = σ/10 5 such uncertainties behave like k -0.5 , as in the single arc case.

The product of orbit determination with constrained multi-arc strategy is a collection of orbits, each one relative to an arc, which connect within a certain precision at the central time of the dark period between one arc and the following. One might question the validity of determining such a piecewise orbit in contrast to a single one. Is this "patchwork" of orbits a good approximation of the real one?

The answer is positive and comes from an application of the Shadowing Lemma [START_REF] Anosov | Geodesic flows on closed Riemann manifolds with negative curvature[END_REF], [START_REF] Bowen | ω-limit sets for axiom A diffeomorphisms[END_REF], [START_REF] Pilyugin | Shadowing in Dynamical Systems[END_REF] and [START_REF] Spoto | Shadowing Lemma and chaotic orbit determination[END_REF]). Using the standard notation for this theorem, once we fix , the result of the orbit determination via constrained multi-arc is a δ-psudotrajectory, δ = δ(σ 2 P , σ, σ * ) > 0 and the Shadowing Lemma ensures that there exists an -shadowing, that is a unique orbit, which approximates the result of the least-squares.

Although we only dealt with a model problem, we believe that the results of this article are of particular relevance in practical applications, especially in the case of space missions exploring the outer planets of the Solar systems and their satellites, like the NASA/ESA/ASI Cassini-Huygens mission to Saturn or the ESA JUICE mission to the jovian system. These kinds of missions are characterized by many close approaches with the planet and its moons 5 , and therefore the orbit of the spacecraft is chaotic (the Lyapunov time being of the same order as the time between two subsequent close encounters). To the best of the authors' knowledge, no single-arc solution has been published yet for any past mission. As already suggested in [START_REF] Dirkx | On the contribution of pride-juice to jovian system ephemerides[END_REF], the use of a constrained multi-arc method would allow to control possible numerical instability arising when the number of observations grows, and at the same time it would ensure significant benefits for the determination of the dynamical parameters, e.g. the spherical harmonics coefficients of the gravity field of the planets/moons, tidal parameters, etc., whose measurement is the ultimate goal of the Radio Science experiments on board of space missions.

Figure 1 :

 1 Figure 1: The jump vector d k between arcs k and k + 1 given by (5) is the difference between the backward-propagated current initial conditions of arc k + 1, namely (x, y) k+1 , and the forwardpropagated current initial conditions of arc k, namely (x, y) k , to the central point of the dark interval between the two arcs.

Figure 3 :

 3 Figure3: Formal uncertainties of the initial conditions of arc 0 (x 0 , y 0 ) and the dynamical parameter µ 0 as a function of the total number of processed arcs, using a constrained multi-arc strategy (log-log plot), chaotic case. The four figures refer to four different experiments repeated with the same setting (x 0 , y 0 ) = (3, 0), µ 0 = 0.5, except for different values of the convergence parameter σ * (see formula 7). The formal uncertainty of µ 0 is a polynomial function of the number of processed arcs, k a ,.

Figure 4 :

 4 Figure4: Formal uncertainties of the initial conditions of arc 0 (x 0 , y 0 ) = (2, 2) and the dynamical parameter µ 0 = 0.5 as a function of the total number of processed arcs, using a pure multi-arc strategy, ordered case (log-log plot). The slope of the µ 0 uncertainty is ∼ -0.5.

  σ * = σ/10 5

Since its introduction, Chirikov standard map has been used in Celestial Mechanics as a test problem to understand features of more complicated nearly-integrable and dissipative models

Weighing observations is a common technique used to take into account the possibility that some observations are more accurate than the others, for instance because they were realized with better instrumentation.

From another point of view, this is equivalent to having a different spacecraft for each arc.

The Lyapunov time of an orbit is the inverse of the largest Lyapunov exponent of that orbit and indicates the time after which the distance between two nearby trajectories increase by a factor e. The computation of the Lyapunov times of the orbits of Chirikov standard map is described in[START_REF] Spoto | Shadowing Lemma and chaotic orbit determination[END_REF].

For instance, the Cassini mission had 127 fly-bys with Titan alone.


François Nicot

Luc Sibille

Félix Darve

Oana Cazacu

Multiscale Modeling of Heterogenous Materials: From Microstructure to Macro-scale Properties zy

Bifurcation in Granular Materials:

A Multiscale Approach

Introduction

The notion of failure can be encountered in many fields, irrespective of the scale considered. This notion is essential in materials sciences where failure can be investigated on the specimen (the material point) scale. It is also meaningful in civil engineering to prevent or to predict the occurrence of failure on a large scale.

For geomaterials, known as non-associate materials, several failure modes can be encountered strictly within the plastic surface. From a mathematical point of view, this feature is essentially related to the non-symmetry of the tangent constitutive tensor. Whereas the localized mode describes a failure corresponding to a discontinuous displacement field, the diffuse mode is associated with a homogenous kinematic field with no localization pattern. It has been shown that this failure mode can be predicted by the vanishing of the second-order work [NIC 07a and b]. Introduced by Hill [HIL 58], this quantity (hereafter denoted by 2 W ) is defined from the inner product of the incremental first Piola-Kirchoff stress tensor ( Π ) with the incremental displacement gradient tensor:

() ∫ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ Π = o V o j i ij dV X u W δ δ
The interest of this semi-Lagrangian formulation lies in that all variables are reported to the fixed initial configuration defined by the volume o V and the coordinates i X . For a material point corresponding to a representative volume element (RVE) of a granular material, equation [11.1] simplifies into the following expression:

( ) ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ Π = j i ij o X u V W δ δ 2 [11.2]
The second-order work can also be expressed under an Eulerian formulation, introducing the Cauchy stress tensor σ :

() ( ) t L t L u div V W t δ δ σ σ δ σ δ : 2 ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - + = [11.3]
where ij L , such as According to the small strain approximation, Dt δδ ε

( ) i ij j u Lt x δ δ ∂ = ∂ ,

=

, where 1 2

t DL L ⎛⎞ =+ ⎜⎟ ⎝⎠
denotes the strain rate tensor, and δε is the incremental small strain tensor. From the symmetry of this tensor, it follows that:

() 2 2 2 : : : t L V V V W δ σ ε δ σ δ ε δ σ δ - + = [11.4]
The different formulations of the second-order work introduce macroscopic tensorial variables that represent both the complex force and displacement distributions within the granular specimen. The vanishing of the second-order work stems therefore from microstructural origins for which the local variables (contact forces and relative displacements between adjoining particles) are relevant. As a consequence, since the vanishing of the second-order work is a proper criterion for detecting the occurrence of a certain failure mode in geomaterials, it makes sense to track the microstructural origin of this macroscopic vanishing. This analysis will be first carried out by considering our micro-directional model [NIC 05], which is basically a micromechanically based constitutive relation. Then, this approach will be extended based on a general micromechanical derivation.

Microstructural origin of the vanishing of the second-order work

The micro-directional model

The micro-directional model is a multi-scale relation between the Cauchy stress tensor dσ and the strain tensor dε by taking micro-mechanical characteristics into account. In this approach, the granular assembly is described as a distribution of contacts within adjoining particles. Each contact is associated with a given direction of the physical space, corresponding to the normal direction to the tangent contact plane. The texture is therefore described by the distribution of contacts along each direction of the physical space. The probability that some contacts exist in a given direction is investigated and local variables are averaged in each direction, so that directional variables are introduced. Fundamentally, this model is based on a homogenization procedure within a representative volume element (RVE) that can be resolved in the three following basic stages (for more details, see [NIC 05]).

∑ = = c N c c j c i ij l F V 1 1 σ [11.5]
where c l is the branch vector joining the centers of particles in contact on contact c, c F is the contact force, and the sum is extended to all the c N contacts occurring in the RVE of volume V . The norm of the branch vector c l is assumed to be a constant parameter (equal to the mean diameter of the grains) whose evolution over the loading programs is ignored. This ensures that the terms c F and c l are uncoupled. The discrete summation given in equation [11.5] can be replaced with a continuous integration over all the contact directions in the physical space. This scheme confers the directional character to the model:

∫∫ Ω = D j i g ij d n F r ω σ 2 [11.6]
where ω is the density of contacts along each space direction n , g r denotes the mean radius of the sphere-shaped grains, F is the average of all contact forces c F associated with contacts oriented in the direction n , and dΩ is the elementary solid angle. After differentiation it follows that:

∫∫ ∫∫ Ω + Ω = D j i g D j i g ij d n F r d n F r δω ω δ δσ 2 2 [11.7]
The kinematic projection relation is given by: ()

j ij g i n r n u δε δ 2 ˆ= [11.8]
where () ûn is the directional kinematic variable linked to ()

F n along the contact direction n .
The local behavior is described by introducing a constitutive relationship between both average normal force ˆn F and tangential force ˆt F and both average relative normal displacement ˆn u and tangential displacement ˆt u . An elastic-plastic model is introduced, and the following local constitutive incremental relations can be inferred:

n n n u k F δ δ = [11.9a] t t t t t t t t F u k F u k F F ˆ - + + = δ δ ξ δ [11.9b] where ( ) { } ˆm in , tan ttt gnnn F ku F k u ξ δ φ δ =+ + , n k is the normal elastic stiffness, t
k is the tangential elastic stiffness, and g ϕ is the local friction angle.

Microstructural expression of the macroscopic second-order work

Starting from equation [11.7], and noting that the density of contact ω along each direction is expressed as

e V ω ω =
, where e ω is the number of contacts along the considered direction, it follows that the differentiation of the Cauchy stress tensor is the sum of three terms:

ij D e j i g D e j i g ij V V d n F V r d n F V r σ δ δω ω δ δσ - Ω + Ω = ∫∫ ∫∫ 2 2 [11.10]
which can also be written as:

∫∫ ∫∫ Ω = Ω - + D e j i g D e j i g ij ij d n F V r d n F V r V V ω δ δω σ δ δσ 2 2 [11.11]
Now, taking advantage of the kinematical projection relation yields: δδ can be interpreted as the microscopic second-order work associated with the contact "c" between two given adjoining particles. As a consequence, equation [11.13] states that the macroscopic second-order work is equal to the sum of the microscopic second-order works associated with all the contacts existing within the assembly.

= Ω - + ∫∫ D e j ij i g ij ij ij ij d d n F r V V ω δε δε σ δ δε δσ 2 ∫∫ Ω D e i i d u F ω δ δ
This basic result was inferred by considering a given constitutive relation, i.e. the micro-directional model. The purpose of the next section consists of generalizing this result without referring to any constitutive model.

From micro to macro second-order work

Let us consider a granular assembly containing N grains "p", with 1 pN ≤≤ . Each grain "p" is in contact with p n other adjoining grains "q", with ( )

1 t r () y r () x r () z r frame ℜ frame ℜ () 2 t r

Figure 11.2. Galilean frame and local frame

The microscopic second-order work attached to the contact "c" between particles "p" and "q" is given by the relation ([NIC 07a and c]):

q p c q p q p u F W , , , 2 ˆ δ δ ⋅ = [11.14]
where , ˆp q F δ denotes the incremental contact force exerted by particle "p" on particle "q", and , ˆp q c u δ is the incremental relative displacement of particle "p" with respect to particle "q". On the granular assembly scale, the macroscopic second-order work can be related to the second-order time derivative of the kinetic energy as:

() t E u F W c V p p p ext 2 , 2 δ δ δ - ⋅ = ∑ ∂ ∈ [11.15]
Taking into account the expression of the kinetic energy, () ( )

∑ = ⋅ + ⋅ = N p p p p p c M u F t E 1 ω δ δ δ [11.16]
we obtain, after some algebra ([NIC 07a]):

( ) ( ) ∑ ∑∑ ∂ ∈ = - = ⋅ - ⋅ = V p p p ext N p p q q p c q p u F u F W 2 , 1 1 1 , , 2 ˆδ δ δ [11.17]
This basic relation indicates that the macroscopic second-order work is the sum of the microscopic second-order works extended to all the contacts of the whole assembly,

( ) ( ) 2 2 2 2 sin cos tan t c t n c t c n g n c n c u k u u k u k W δ α δ δ α ϕ δ + + = [11.18]
where α is the angle between both vectors 1 ( )

t c t c F t F = and
n c t c n g n c n c u u k u k W δ δ ϕ δ

Some remarks on the basic micro-macro relation for the second-order work

Let us come back to further discuss relation [11.17]. This relation was investigated from discrete element simulations. Considering a cubic granular specimen, at rest after an initial axisymmetric drained triaxial loading, a series of stress probes was imposed, along all the directions of the incremental stress space, and both quantities 2 W (macroscopic second-order work) and 2 W (sum of the microscopic second-order works) were compared.

As seen in Figure 11.5, equation [11.17] is perfectly verified within the elastic tensorial zone (the zone gathering loading directions leading to no plastic dissipation.), when contacts behave essentially elastically, and in a part of the plastic tensorial zone where plastic dissipation is related to sliding of contacts (Figure 11.6). As soon as loading directions are characterized by contact opening and/or creation (which corresponds to the central part of the plastic tensorial zone), a significant shift between 2 W and 2 W exists (Figure 11.6). Should the validity of equation [11.17] be queried? It is our conviction that this basic relation is valid, irrespective of the tensorial zone considered. Nevertheless, it is worth noting that equation [11.17] applies to an equilibrium state; on the contrary, discrete element simulations require considering a finite time interval to calculate both quantities 2 W and 2 W . For loading directions belonging to the central part of the plastic tensorial zone (characterized by contact opening and creation), grain rearrangements continuously take place, so that the medium is (at least locally) no longer in equilibrium. As a consequence, for such loading directions, discrete element simulations do not constitute an appropriate way of checking a relation valid at the equilibrium but involving (force and displacement) rates.

Conclusion

This chapter was devoted to the micromechanical investigation of the vanishing of the second-order work. As this quantity was shown to play a fundamental role in detecting the occurrence of a certain failure mode (diffuse failure mode, related to the spontaneous burst of kinetic energy), it is of great interest to understand what the microstructural conditions that lead to the vanishing of the second-order work are.

First, by considering our micro-directional model, it was inferred that the macroscopic second-order work of a given granular assembly is equal to the sum of the microscopic second-order works extended to all the contacts within the assembly. Then, the validity of the relation was extended based on general micromechanical arguments. This relation is fundamental since it bridges both microscopic and macroscopic worlds. The analysis was pursued by introducing an elastoplastic (frictional) model on the contact scale. The conditions for the vanishing of the microscopic second-order work (which is quadratic form with respect to the relative displacement) were examined, and an interpretation of the fact that the "unstable cones" containing the loading directions of the incremental stress space corresponding to negative values of 2 W are contained in the third quadrant ( 1 0 δσ < and 3 0 δσ < ) was provided. The microstructural ingredient of the analysis is essentially related to the local sliding condition. An important aspect remains to be considered in relation with the sudden (and brutal) deletion of contacts on the (mesoscopic) force chain scale. This geometrical aspect should be considered in addition to the former material aspect (sliding condition) considered in this chapter.

  is the general term of the velocity gradient tensor L . It follows that the second-order work is the combination of three terms: : VL t δσ δ is a material term, : VL t δσ δ is related to the change in volume, and σδ is associated with the change in the texture [NIC 07a and c].

  The stress average corresponds to the Love formula ([LOV 27]; [WEB 66]; [CHR 81]; [MEH 82]):

2 :

 2 .12] can interestingly be compared to equation[11.4]. The term () 2 VLt σδ , which is shown to be related to the change in texture [NIC 07c], can be assimilated to the term 2 for textural change. In these conditions, it can be established that the macroscopic second-order work can be expressed in a very straightforward manner with respect to microscopic variables: over all the contacts contained within the assembly. As demonstrated by Nicot and Darve ([NIC 07a and c]), the term cc ii F u

Figure 11 . 1 .

 111 Figure 11.1. Granular assembly: boundary particles and external forces

.W

  This last term seems to be negligible from simulations based on a discrete element method [SIB 06 and 07]. This relation, that relates the macroscopic second-order work to microstructural elements embedded in the term 2 , provides insight into the microstructural origins of the vanishing of the secondorder work. The next section is concerned with examining this feature.11.2.4. Micromechanical analysis of the vanishing of the second-order workLet us consider the term constitutive equations such as those given in [11.9]. Considering any contact "c", the contact behaves in the plastic regime, the microscopic second-order work 2 c W is a quadratic form that can be positive or negative:

  and equation [11.18] yields:

WFδFigure 11 . 3 .Figure 11 . 4 .

 113114 Figure 11.3. Evolution of the contact force for the vanishing of the microscopic second-order work: the contact force descends the Coulomb line

Figure 11 . 5 .

 115 Figure 11.5. Microscopic and macroscopic second-order work density along different stress loading directions (after [NIC 07c])


Karthik Murali

Madhavan Rathai

Olivier Sename
email: olivier.sename@gipsa-lab.grenoble-inp.fr

Mazen Alamir
email: mazen.alamir@gipsa-lab.grenoble-inp.fr

A COMPARATIVE STUDY OF DIFFERENT NMPC SCHEMES FOR CONTROL OF SEMI-ACTIVE SUSPENSION SYSTEM

Keywords: Vehicle dynamics, Non-linear model predictive control, Suspension systems, Optimal control

Control of semi-active suspension system for vertical dynamics of automobiles plays a vital role in guaranteeing comfort and safety for the on-board passengers. The seemingly simple task of control poses to be daunting under the presence of multiple nonlinearities, physical constraints and specification over objective satisfaction of the system and thereby, it is of paramount importance to account for these issues during control system design for better efficiency and prolonged endurance of the suspension system. Amongst the existing several control methodologies, predictive control techniques serves as a promising approach in dealing with the aforementioned issues. In this paper, we present several nonlinear model predictive control (NMPC) schemes and a detailed performance analysis of the methods. The incorporated NMPC schemes are the direct methodssingle shooting, multiple shooting and collocation methods. The different methods were tested in simulation under MATLAB environment.

INTRODUCTION

In the recent years, study on advanced control methods for automotive systems has gained a huge momentum and most of the automotive industries has embarked in research and development and implementation of better control algorithms for improved passenger comfort and safety and as well as vehicle's performance in terms of energy efficiency, pollution reduction etc. This sudden surge in interest is partly in response to the advent of autonomous vehicles and it has also been the key driving factor for automotive companies to strive for optimal performance. Conditioned upon the aforementioned requirements both objectively and subjectively, optimal control methods fares much better than other control methods due to the systematic approach embodied in its design to tackle the above issues. In order to implement optimal control in practice, model predictive control (MPC) formulation of the optimal control problem (OCP) provides the necessary feedback control framework to work seamlessly in real-time and real-world. The crux of the MPC scheme is the receding horizon method, where an OCP is solved online at every sampling instant. Despite the enormous benefits of the method, one of the major downside is the computational requirements for solving the underlying optimization problem of the OCP at every sampling period. However, given the exponential growth of computational power and easy/cheap availability of computing resources, the chasm between the two is closing in and the method seems promising in the near future.

The main ingredients of the MPC problem are a) the objective function, b) system constraints and c) the system dynamics [START_REF] Alamir | Stabilization of nonlinear systems using receding-horizon control schemes: a parametrized approach for fast systems[END_REF] and when any of the ingredients induce any nonlinearity, the problem is known as nonlinear MPC (NMPC). In general, the NMPC problem can be solved in three prongs which are a) Direct methods (discretize then optimize), b) Indirect methods (optimize then discretize) and c) Dynamic programming [START_REF] Rawlings | Model Predictive Control: Theory, Computation, and Design[END_REF]. This paper utilizes the direct methods approach, however irrespective of the method adopted the problem requires a nonlinear programming (NLP) solver to deduce the optimal solution for the NMPC problem. The direct methods can be broadly classified as single shooting, multiple shooting and collocation methods. Despite the fact that all the three methods solve the same problem, the solutions sought varies due to the difference in the problem formulation. The NLP solver that yields the solution might return different local minima solution for different cases. Under special case of convex problems, the solution of all the methods tends to be the same. In this paper, the subject of focus is on a comparative study on NMPC direct methods for control of semi-active suspension system for a quarter car model. The nonlinearity is induced due to the inherent dissipativity characteristics of the semiactive damper system [START_REF] Rathai | A Parameterized NMPC Scheme for Embedded Control of Semiactive Suspension System[END_REF], [START_REF] Savaresi | Semi-active suspension control design for vehicles[END_REF].

The paper is organized as follows. Section 2 describes the mathematical model of the quarter car model equipped with semi-active suspension system. Section 3 and Section 4 details the NMPC design requirements and direct methods formulation. Section 5 expounds the simulation results obtained and finally, the paper is concluded with conclusions and future works in Section 6.

MATHEMATICAL MODELING

The quarter-car model system equipped with semi-active suspension system shown in Fig. 1 are given by the following set of equations [START_REF] Alamir | Stabilization of nonlinear systems using receding-horizon control schemes: a parametrized approach for fast systems[END_REF] Where, , , , are the sprung mass, unsprung mass, stiffness coefficient of the suspension system and stiffness coefficient of the tyre respectively. , , , = , = are the chassis mass displacement, unsprung mass displacement, road profile displacement, deflection velocity and displacement between the chassis and the tyre respectively. The force exerted due to the suspension system is given by the Guo's [START_REF] Guo | Dynamic modeling of magnetorheological damper behaviors[END_REF] nonlinear equation with , , , the appropriate parameters for the force model. The input for the system is , which is the PWM duty cycle (DC) signal that drives the damper system. The nonlinear dynamics of the system is compactly expressed with , where , , the road disturbance acting on the system and is the PWM-DC signal.

Fig. 1 The quarter-car model

The parameters for the mathematical model are utilized from the INOVE test platform. The INVOVE test platform [START_REF] Vivas-Lopez | INOVE: a testbench for the analysis and control of automotive vertical dynamics[END_REF], shown in Fig. 2 is a 1:5-scaled baja style racing car which consists of 4 controllable Electro-Rheological (ER) dampers and 4 DC motors to generate different road profiles for each wheel corner. The numerical values of the parameter are listed in the reference herein [START_REF] Rathai | A Parameterized NMPC Scheme for Embedded Control of Semiactive Suspension System[END_REF].

Objective requirements

In this paper, the chosen objective design for the semi-active suspension pertains to comfort objective. The prime goal of the comfort based objective design is to minimize the vertical acceleration of the chassis [START_REF] Rathai | Reachability based model predictive control for semiactive suspension system[END_REF] i.e. . For a given look ahead period , the objective is defined with

(2)
Where, which is obtained from equation (1).

Constraint requirements

The constraints incorporated into NMPC problem are:  Semi-active damper input constraint: Minimum and maximum saturation force of the semi-active damper system i.e. .  State constraint: Minimum and maximum stroke displacement of the suspension system i.e. .  PWM-DC input constraint: The PWM-DC input is constrained in the set .  Road disturbance: The road profile is treated constant over the horizon i.e.

This can be appended as an additional state variable and the augmented dynamics of the system can be expressed with , where , .

The list of constraints can be compactly expressed with ) , which is the mixed nonlinear input and state constraints of the system. Thus, given the objective and the constraints of the system, the OCP to be solved at every time instant with initial condition for the NMPC problem is defined with

(3)

NMPC DIRECT METHODS FORMULATION

The OCP of the NMPC scheme (3) leads to an infinite dimensional problem and it is cumbersome and in most cases it is impractical to implement. The direct methods, transcribes the problem into a finite dimensional problem, which yields an approximate solution to (3) by virtue of NLP solvers. In order to set the transcription procedure, the following assumptions are taken into consideration [START_REF] Johansen | Introduction to nonlinear model predictive control and moving horizon estimation[END_REF]:

Assumptions

 The input is finitely parameterized by piecewise constant vector at an integer multiple of the sampling period over the horizon. With this representation, the input can be expressed with , which is a piecewise continuous input signal.  The dynamics (defined by the ODE in (3)) of the system is numerically simulated for the given input signal , which is compactly expressed with , which is evaluated at discrete time instants . The time stamps typically corresponds to the integer multiple of the sampling period and it is utilized to discretize the dynamics, objective function and the constraint functions listed in [START_REF] Rathai | A Parameterized NMPC Scheme for Embedded Control of Semiactive Suspension System[END_REF]. The ODE solver utilized in this paper is a 4 th order Runge-Kutta (RK) solver. With the aforementioned assumptions, the generic NLP framework for the direct method is expressed with [START_REF] Savaresi | Semi-active suspension control design for vehicles[END_REF] Where, is the optimization variable which depends upon the direct method formulation utilized.

Remark 1. The ODE solver mentioned in the assumption is to be considered as a computer code and not in terms of algebraic equations. The ODE solver is a simulator which takes the numerical input trajectory and outputs the numerical state trajectory which is utilized in the objective and constraint functions for the NLP problem mentioned in (4). Thus, the objective and the constraint functions is embedded with the ODE solver code and is ought to be deemed as computer codes (function code). Under conditions of twice differentiability of all the functions (codes) listed in (3), the Jacobians and Hessians for the NLP solver are numerically obtained by methods such as finite differences, algorithmic differentiation etc. [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF] (also known as oracles in optimization parlance) and this information aids the optimization procedure. Thanks to MATLAB's "fmincon" routine which automatically computes the derivatives from the objective and constraint functions.

Direct single shooting

The direct single shooting method also known as sequential method eliminates the dynamics equality constraint in (3) by forward simulation and thus, removing the state variables from the OCP NMPC problem. This reduces the optimization problem only to the input variables , which is obtained from the following NLP problem.

(

) 5
The objective is discretized by means of Riemann sum at time stamps . The states are replaced with the ODE simulator evaluated at these time stamps in the objective as well as the constraints. The optimal solution obtained from ( 5) and the first input is injected into the system and the process is repeated in receding horizon manner.

Direct multiple shooting

The direct multiple shooting method also known as simultaneous method retains the state variables as optimization variables and this increases the number of decision variables in the optimization formulation in [START_REF] Rathai | A Parameterized NMPC Scheme for Embedded Control of Semiactive Suspension System[END_REF]. The ODE solver simulates the system over multiple time intervals i.e.

simultaneously and the final state value for the simulation in the each interval is stipulated to obey the dynamics of the system, which is enforced by equality constraints. The NLP optimization problem is formulated as [START_REF] Vivas-Lopez | INOVE: a testbench for the analysis and control of automotive vertical dynamics[END_REF] The optimal solution for the above optimization problem yields both the optimal state trajectory and optimal input sequence. As per the standard receding horizon policy, the first input is applied to the system and repeated in the future. The benefits of utilizing direct multiple shooting methods include a) Better simulator stability with unstable system, b) Parallelizability of ODE simulation, c) Structural properties of the Hessian matrices aid the optimization routine. However the flip side is that the optimization is carried out over an increased number of variables and a good initialization for the NLP solver is required for faster convergence to the optimal/suboptimal solution (this can be ameliorated by warm start procedure).

Direct collocation

Direct collocation methods are extension to the simultaneous methods, where the ODE simulator is expunged from the multiple shooting formulation [START_REF] Vivas-Lopez | INOVE: a testbench for the analysis and control of automotive vertical dynamics[END_REF] and is replaced with algebraic equality constraints enforced at the collocation points. The optimization problem is expressed with [START_REF] Rathai | Reachability based model predictive control for semiactive suspension system[END_REF] The fundamental difference between ( 6) and ( 7) is (Implicit solver), which satisfies the dynamics of the system at the collocation points. Direct collocation methods are typically suited for stiff systems and implicit RK methods are popularly adopted in optimal control literature. This introduces additional algebraic variables which are casted as equality constraints in [START_REF] Rathai | Reachability based model predictive control for semiactive suspension system[END_REF]. In this paper, vis-à-vis to the system considered, trapezoid collocation method is utilized to enforce the dynamics constraints at the collocation points.

SIMULATION RESULTS

The NMPC methods were implemented in MATLAB environment and NLP solver utilized was the "fmincon" routine. The sampling period was chosen to be 5ms and the look ahead period as 15ms. The input parameterization was chosen to be a constant signal over the horizon. The "fmincon" solver was set to sequential quadratic programming (SQP) mode with maximum number of Newton iteration as 3. The road profile utilized was a chirp signal with an amplitude of 2.5mm and a frequency sweep from 1Hz to 8Hz for a duration of 25s. 4 displays the chassis acceleration for the system. From the plot is evident that the solution for the single shooting, multiple shooting and collocation method are nearly equal in performance, however the collocation method tends to be have better greater RMS value in comparison with the other two methods. However, this is subjected to choice of the collocation method utilized in the problem formulation. The RMS values of the acceleration is listed in the table Table I. The dissipativity constraints are illustrated in Fig. 5.

NMPC method RMS value (m/s^2) Single shooting 6.9020 Multiple shooting 6.8808 Collocation 6.9685 Table I Chassis acceleration RMS values

Fig. 2

 2 Fig.2 INOVE test platform

Fig. 3

 3 Fig.3 PWM duty cycle input

Fig. 3

 3 Fig.3displays the PWM-DC computed from the NMPC methods. Clearly, it is evident that the input profiles are not the same for the methods due to the difference in the problem formulation.

Fig. 4

 4 Fig.4Chassis acceleration Fig.4displays the chassis acceleration for the system. From the plot is evident that the solution for the single shooting, multiple shooting and collocation method are nearly equal in performance, however the collocation method tends to be have better greater RMS value in comparison with the other two methods. However, this is subjected to choice of the collocation method utilized in the problem formulation. The RMS values of the acceleration is listed in the table TableI. The dissipativity constraints are illustrated in Fig.5.NMPC method RMS value (m/s^2) Single shooting 6.9020 Multiple shooting 6.8808 Collocation 6.9685 Table I Chassis acceleration RMS values

Fig. 5

 5 Fig.5 Dissipativity constraint

CONCLUSION AND FUTURE WORKS

In this work, a comparative simulation study is conducted on the different popular NMPC methods and its integration with a quarter car model equipped with semi-active suspension system. The work is conducted in the spirit of learning, exploring and investigating different NMPC methods and its suitability for automotive systems. In the future works, real-time implementation of the above methods are to be conducted to validate the real-time feasibility and viability in the INOVE test platform at GIPSAlab, Grenoble.

ACKNOWLEDGEMENT

This work was supported within the ITEA3 European project, 15016 EMPHYSIS (Embedded systems with physical models in the production code software).