Roberta Lima

Christian Soize
email: christian.soize@univ-paris-est.fr

Rubens Sampaio
email: rsampaio@puc-rio.br

Robust design optimization with an uncertain model of a nonlinear vibro-impact electro-mechanical system

Keywords: electro-mechanical systems, vibro-impact, robust design optimization, nonlinear dynamics

In this paper, the robust design with an uncertain model of a vibro-impact electro-mechanical system is done. The electro-mechanical system is composed of a cart, whose motion is excited by a DC motor (motor with continuous current), and an embarked hammer into this cart. The hammer is connected to the cart by a nonlinear spring component and by a linear damper, so that a relative motion exists between them. A linear flexible barrier, placed outside of the cart, constrains the hammer movements. Due to the relative movement between the hammer and the barrier, impacts can occur between these two elements. The developed model of the system takes into account the influence of the DC motor in the dynamic behavior of the system. Some system parameters are uncertain, such as the stiffness and the damping coefficients of the flexible barrier. The objective of the paper is to perform an optimization of this electro-mechanical system with respect to design parameters (spring component, and barrier g) in order to maximize the impact power under the constraint that the electric power consumed by the DC motor is lower than a maximum value. This optimization is formulated in the framework of robust design due to the presence of uncertainties in the model. The set of nonlinear equations are presented, and an adapted time domain solver is developed. The stochastic nonlinear constrained design

Introduction

Electro-mechanical systems are common in actual technologies, and their design is of a great interest in many areas. Many works have been done in this topic, as [START_REF] Zhankui | Nonlinear and chaos control of a micro-electromechanical system by using second-order fast terminal sliding mode control[END_REF][START_REF] Sadeghian | Comparison of generalized differential quadrature and galerkin methods for the analysis of micro-electromechanical coupled systems[END_REF][START_REF] Lee | Design and analysis of electromechanical characteristics of micromachined stainless steel pressure sensor[END_REF], trying to characterize the mutual interaction between electrical and mechanical parts. This interaction leads us to analyze very interesting nonlinear dynamical systems (see for instance [START_REF] Rocard | Dynamique Générale des Vibrations[END_REF][START_REF] Kononenko | Vibrating Systems with a Limited Power Supply[END_REF][START_REF] Nayfeh | Nonlinear Oscillations[END_REF][START_REF] Fidlin | Nonlinear Oscillations in Mechanical Engineering[END_REF][START_REF] Hagedorn | Non-linear Oscillations[END_REF][START_REF] Evan-Iwanowski | Resonance oscillations in mechanical systems[END_REF][START_REF] Cartmell | Introduction to Linear[END_REF]) , in which the nonlinearities vary with the coupling conditions, and also affects the two most important variables used to evaluate the performance of electro-mechanical systems, related to the power consumed by the electrical part, and the power used into the movement of the mechanical part. As the mutual interaction between electrical and mechanical parts affects the two powers used to evaluate the system performance, the coupling effects must be analyzed in the design optimization problem for electro-mechanical systems. The present work deals with the robust design optimization of a vibro-impact electric-mechanical system in order to improve its performance. The electrical part of the system is a DC motor, and the mechanical part is a vibro-impact system. It should be noted that, in [START_REF] Lima | Stochastic analysis of an electromechanical coupled system with embarked mass[END_REF], the equations and the numerical integration were presented for a similar electric-mechanical system for which the embarked mass was replaced by a pendulum and for which there was no impact. This first work has allowed the electro-mechanical coupling to be analyzed as a function of the mass of the mechanical system. The analysis of vibro-impact systems is not a new subject, and is frequently encountered in technical applications of mechanisms. The interest of analyzing the optimization of their performance is reflected by the increasing amount of research in this area (see for instance [START_REF] Luo | Periodic-impact motions and bifurcations of vibro-impact systems near 1:4 strong resonance point[END_REF][START_REF] Ostasevicius | Numerical analysis of dynamic effects of a nonlinear vibro-impact process for enhancing the reliability of contact-type mems devices[END_REF][START_REF] Yue | Global analysis of boundary and interior crises in an elastic impact oscillator[END_REF], and also the book by Ibrahim [START_REF] Ibrahim | Vibro-Impact Dynamics: Modeling, Mapping and Applications[END_REF], which is completely devoted to this problem). Besides the theoretical research in vibro-impact dynamics, the applications have also been developed, such as percussive actions in rotary drilling, vibration hammer, impact damper, and gears. The vibro-impact dynamics appears also in several other situations, as for example in earthquakes, where the interest is the seismic mitigation [START_REF] Nucera | Targeted energy transfers in vibro-impact oscillators for seismic mitigation[END_REF]. The vibro-impact dynamics can be affected by many factors, and the analysis of vibro-impact dynamical systems requires to taken into account uncertainties in the computational models that are used (see for instance [START_REF] Sampaio | On measures of non-linearity effects for uncertain dynamical systems -application to a vibro-impact system[END_REF]).

The objective of the paper is to perform an optimization of this electromechanical system with respect to design parameters that are the stiffness of the spring component, and the gap, g, of the barrier. The optimization consists in maximizing the impact power under the constraint that the electric power consumed by the DC motor is lower than a maximum value. This optimization is formulated in the framework of robust design due to the presence of uncertainties in the computational nonlinear dynamics model of the electro-mechanical system.

This paper is organized as follows. In Section 2 , without introducing any simplifying hypotheses, the equations are constructed. The initial value problem is presented for the vibro-impact electro-mechanical system. In Section 3 we define the variables of interest for the design optimization. A sensitivity analysis with respect the system parameters is performed numerically in Section 4, in order to define the most sensitive parameters, and then to define what will be the design parameters. The construction of the probabilistic models of the uncertain parameters, and the formulation of the robust design optimization problem are given in Sections 5 and 6. The numerical results of the robust design optimization problem are presented in Section 7.

Dynamic of the coupled system

First, the elements (motor, cart, hammer, and barrier) of the electromechanical system are presented. The coupling between the motor and the mechanical system is described, and the time-evolution coupled problem is mathematically formulated as an initial-value problem.

Electrical system: DC motor

The mathematical modeling of DC motors is based on the Kirchhoff law [START_REF] Karnopp | System Dynamics: Modeling and Simulation of Mechatronic Systems[END_REF]. It is written as

l ċ(t) + r c(t) + k e α(t) = ν , j m α(t) + b m α(t) -k e c(t) = -τ (t) , ( 1
)
where t is time, ν is the source voltage, c is the electric current, α is the angular speed of the motor, l is the electric inductance, j m is the inertia moment of the motor, b m is the damping ratio in the transmission of the torque generated by the motor to drive the coupled mechanical system, k e is the constant of the motor electromagnetic force, and r is the electrical resistance. Figure 1 shows a sketch of the DC motor. The available torque delivered to the coupled mechanical system is represented by τ that is the component of the available torque vector τ in the z-direction shown in Fig. 1.

Motor-cart-hammer coupled system

As described in the introduction, the system is composed by a cart whose movement is driven by the DC motor, and by a hammer that is embarked into the cart. The motor is coupled to the cart through a pin that slides into a slot machined in an acrylic plate that is attached to the cart, as shown in Fig. 2. The off-center pin is fixed on the disc at distance Δ of the motor shaft, so that the motor rotational motion is transformed into a cart horizontal movement. To model the coupling between the motor and the mechanical system, the motor shaft is assumed to be rigid. Thus, the available torque vector to the coupled mechanical system, τ , can be written as

τ (t) = Δ(t) × f(t) , (2)
where Δ = (Δ cos α(t), Δ sin α(t), 0) is the vector related to the eccentricity of the pin, and where f is the coupling force between the DC motor and the cart. Assuming that there is a viscous friction between the pin and the slot, the vector f has two components: the horizontal force that the DC motor exerts in the cart, f x , and the vertical force, f y , induced by the viscous friction. The available torque τ and vertical force f y are written as

τ (t) = f y (t) Δ cos α(t) -f x (t) Δ sin α(t) , ( 3
)
f y (t) = c pin Δ α(t) cos α(t) , ( 4
)
where c pin is the viscous friction. The embarked hammer is modeled as a rigid body of mass m h and its relative displacement is h with respect to the cart.

In the adopted model, the constitutive equation of the spring component between the hammer and the cart is written as

f s (t) = k h1 h(t) + k h3 h(t) 3 .
The rate of nonlinearity of the hammer stiffness is defined as

r h = k h3 /k h1 .
We introduce the natural frequency, ω h , of the hammer suspended to the linear spring with constant k h1 such that ω h = k h1 /m h . The horizontal cart displacement is represented by x. Due to constraints, the cart is not allowed to move in the vertical direction. The spring-damper element modeling the medium on which the impacts occur, is constituted of a linear spring with stiffness coefficient k i and a damper with damping coefficient c i . The equations of the cart-hammer-barrier system are

ẍ(t) (m c + m h ) + ḧ(t) m h + c ext ẋ(t) = -f imp (t) + f x (t) , (5) ẍ
(t) m h + ḧ(t) m h + c int ḣ + k h1 h(t) + k h3 h 3 (t) = -f imp (t) , (6)
where, c ext is the viscous friction coefficient between the cart and the rail and

c int = 2ς int √ m h k h1
is the viscous friction coefficient between the cart and the hammer (ς int is the damping ratio). The term f x is the horizontal coupling force between the DC motor and the cart, and f imp is the impact force between the hammer and the barrier, which is written as

f imp (t) = -φ(t) k i (x(t) + h(t) + g) + c i ( ẋ(t) + ḣ(t)) , (7)
where

φ(t) = 1, if x(t) + h(t) + g < 0 and ḣ(t) + ẋ(t) < 0 , 0, in all other cases , ( 8
)
in which g is defined as the horizontal distance from the hammer (when α = π/2 rad) to the equilibrium position of the barrier. In the model defined by Eq. ( 8), an impact starts when x(t)+h(t) is negative and equal to -g and, ḣ(t)ẋ(t) < 0. During an impact, the action of the barrier on the hammer stops as soon as the total velocity ḣ(t) + ẋ(t) becomes positive (the return of the hammer). Due to the system geometry, x(t) and α(t) are related by the following constraint

x(t) = Δ cos (α(t)) . ( 9
)
Substituting Eqs. ( 3) to ( 9) into Eq. ( 1), we obtain the initial value problem for the motor-cart-hammer-barrier coupled system that is written as follows.

Given a constant source voltage ν, find (α, c, h) such that, for all t > 0,

l ċ(t) + rc(t) + k e α = ν , ( 10
) α(t) j m + (m c + m h )Δ 2 sin (α(t)) 2 -ḧ(t) (m h Δ sin (α(t))) -k e c(t) + α(t) b m + α(t)(m c + m h )Δ 2 cos (α(t)) sin (α(t)) +c pin Δ 2 cos (α(t)) 2 -c ext Δ 2 sin (α(t)) 2 = φ k i (Δ cos (α(t)) + h + g) + c i (-d α(t) sin (α(t)) + ḣ(t)) Δ sin (α(t)) , (11) ḧ(t)m h -α(t) (m h Δ sin (α(t))) -α(t) (m h Δ α(t) cos (α(t))) + ḣ(t)c int + k h1 h(t) + k h3 h 3 (t) = φ(t) k i (Δ cos (α(t)) + h + g) + c i (-Δ α(t) sin (α(t)) + ḣ(t)) , (12)
where

φ(t) = 1, if Δ cos α(t) + h(t) + g < 0 and ḣ(t) -Δ α(t) cos (α(t)) < 0 0, in all other cases , ( 13
) with the initial conditions,

α(0) = 0 , α(0) = 0 , c(0) = ν r , h(0) = 0 , ḣ(0) = 0 . (14)

Measure of the system performance

At time t, the electric power introduced by the electrical grid in the motor is

π in (t) = ν c(t) . ( 15
)
Let t j b and t j e be the instants of begin and end of the j-th impact, such that for all t belonging to [t j b , t j e ], we have ẋ(t) + ḣ(t) < 0. At time t, the impact power, π j imp (t), is then written as

π j imp (t) = k i (x(t) + h(t)) ( ẋ(t) + ḣ(t)), t j b ≤ t ≤ t j e . ( 16
)
The time average of the impact power during the j-th impact, π j imp , is written as

π j imp = 1 t j e -t j b t j e t j b π j imp (t) dt . ( 17
)
The sum, π imp , of the averages of the impact powers, which is one of the variable of interest in the design optimization problem, is written as

π imp = N imp j=1 π j imp , (18)
where N imp is the total number of impacts that occur during time interval [0, T ]. The time average of the electric power consumed in this time interval is

π elec = 1 T T 0 π in (t) dt . ( 19
)
These two variables, π imp and π elec , are considered for measuring the system performance. The biggest is π imp and the smaller is π elec , better will be the system performance.

Sensitivity analysis and choice of the design parameters

To understand the role played by each system parameter in π imp and π elec , a sensitivity analysis has been done. The objective was to determine what were the system parameters that had the biggest influence in π imp and π elec , in order to define those that will be the design parameters for the robust design optimization problem. The initial value problem defined by Eqs. [START_REF] Cartmell | Introduction to Linear[END_REF] to [START_REF] Yue | Global analysis of boundary and interior crises in an elastic impact oscillator[END_REF] has been rewritten in a dimensionless form for computation and some dimensionless parameters were defined. However, in the sensitivity analysis, these dimensionless parameters were not considered as varying parameters since they do not have an easy physical interpretation. The varying parameters used for the numerical simulations are related with the design of the cart and the embarked hammer. They are:

• m c /m h , relation between the hammer mass and the cart mass;

• k h1 /m h , relation between the linear stiffness of the spring component and hammer mass (a sort of natural frequency of the hammer);

• g, horizontal distance from the hammer (when α = π/2 rad) to the equilibrium position of the barrier;

• Δ, eccentricity of the pin. This parameter determines the length of the cart path.

The other parameters, related with the motor properties and viscous friction coefficients, are fixed and the values of these fixed parameters are given in Table 1. The output responses are π imp and π elec . For computation, the initial 10) to [START_REF] Yue | Global analysis of boundary and interior crises in an elastic impact oscillator[END_REF] has been rewritten in a dimensionless form. The main objective was to reduce the computation time. Duration is chosen as T = 10.0 s. The 4th-order Runge-Kutta method is used for the time integration scheme for which we have implemented a varying time-step.

The time-step is adapted to the state of the dynamical system according to the occurrence or the non occurrence of impacts. When the hammer is not impacting the barrier, the time step used is 10 -4 s, but when the hammer is approaching the barrier and when it is impacting it, the time step is chosen as the value 10 -5 s. Simulations with different values to the initial conditions, were performed. As it was verified that they do not have a significant influence in π imp and π elec , in all simulations the initial conditions were taken as constant, given by Eq. ( 14). phenomenon does not occur with respect to the parameter m c /m h . Thus, while π imp and π elec are not very sensitive to m c /m h , they are sensitive to k h1 /m h , g and Δ. It is also seen, that two different kinds of sensitivity can be distinguished among these three parameters. For parameters k h1 /m h and g, it can be seen that π imp and π elec reach their maxima when k h1 /m h and g are equal to 1, 580 rad 2 /s 2 and 0.011 m respectively. For parameter Δ varying in its range of values, Figs. 5(b) and 7(b) show that the highest is Δ, the highest are π imp and π elec . It has been considered as not necessarily to verify the behavior of π imp and π elec for a larger range of Δ because the value Δ = 0.013 m is already sufficiently large when compared with the system dimensions and the motor properties. It should be noted that if parameter Δ were increased more then, the nonlinearities would increase more, that is not the objective of the analysis. Considering that m c /m h does not have a significant influence in π imp and π elec , and considering that the sensitivity of the parameter Δ is easily predictable, these two parameters will not be considered as design parameters in the robust design optimization problem. Only parameters g and k h1 /m h will thus be considered as design parameters.

Construction of the probability model

As explained in the introduction, this paper deals with the robust design of the electro-mechanical system in presence of uncertainties in the computational model. The three parameters that are assumed to be uncertain are k h1 , k i and c i , which are modeled by the independent random variables K h1 , K i and C i . The probability distribution of each one is constructed using the Maximum Entropy Principle [START_REF] Jaynes | Information theory and statistical mechanics[END_REF][START_REF] Shannon | A mathematical theory of communication[END_REF]. This principle allows the probability distribution of a random variable to be constructed using only the available information, avoiding the use of any additional information that introduces a bias on the estimation of the probability distribution. If a large amount of experimental data are available, then the nonparametric statistics can be used. If there are no available experimental data, or if there are only a few experimental data, then the Maximum Entropy from Information Theory is the most efficient tool for constructing a prior probability model for which its hyperparameters are used as sensitivity parameters if there are no experimental data, or, are estimated in using the parametric statistics if a few experimental data are available. The Maximum Entropy Principle states: out of all probability distributions consistent with a given set of available information, choose the one that has maximum uncertainty (the Shannon measure of entropy). The available information of the random variables is defined as 1. K h1 , K i and C i are positive-valued random variables, 2. the mean values are known: 3. in order that the response of the dynamical system be a second-order stochastic process, we impose the following conditions:

E{K i } = K i , E{C i } = C i and E{K h1 } = K h1 ,
E{log K i } < ∞, E{log C i } < ∞ and E{log K h1 } < ∞.
Thus, the Maximum Entropy Principle for each random variable K i , C i , and K h1 , yields a Gamma distribution (see [START_REF] Soize | Maximum entropy approach for modeling random uncertainties in transient elastodynamics[END_REF]),

p(a) = 1 [0,+∞) (a) 1 μ 1 δ 2 1 δ 2 1 Γ(1/δ 2 ) a μ 1 δ 2 -1 exp a δ 2 μ , ( 20
)
where 1 [0,+∞) (a) is an indicator function that is equal to 1 for a ∈ [0, +∞) and 0 otherwise, and where

• Γ is the Gamma function: Γ(b) = ∞ 0 t b-1 exp(-t)dt;
• δ = σ μ is the coefficient variation of the random variable, μ is its mean value representing K i , C i , or K h1 , and σ is its standard deviation.

Robust design optimization problem

In order to formulate the robust design problem, the set of all the system parameters is divided into three subsets. The first subset is the family of the fixed parameters that is represented by the vector p fix = { ν, l, r, j m , k e , b m , c pin , c ext , ς int , r h , m c , m h , Δ }. The second one is the family of the design parameters that is represented by the vector p des = {K h1 /m h , g}. The third one is the family of the uncertain parameters that is represented by the random vector P unc = {K i , C i , K h1 }. Since P unc is a random vector, the outputs of the electro-mechanical system are stochastic processes and, consequently, π imp (p des , p unc ) and π elec (p des , p unc ), become random variables Π imp (p des ) = π imp (p des , P unc ) and Π elec (p des ) = π elec (p des , P unc ). The cost function of the robust design optimization problem is defined by

J(p des ) = E{Π imp (p des )} . ( 21
)
The robust design optimization problem is written as

p opt des = arg max p des ∈C ad J(p des ) , (22)
in which C ad = {p des ∈ P des ; E{Π elec (p des )} ≤ c elec }, where P des is the admissible set of the values of p des , and where c elec is an upper bound.

Results of the robust optimization problem

The hyperparameters δ K i and δ C i , which control the level of uncertainties for K i and C i are fixed to 0.1. The robust design optimization problem is then solved for three levels of uncertainties for K h1 , defined by the following values of the hyperparameters δ K h1 = 0, δ K h1 = 0.1, and δ K h1 = 0.4. The optimization problem is also considered whitout uncertainties in the systems parameters, that is, the deterministic case (δ K h1 = δ K i = δ C i = 0). For p des ∈ C ad , the cost function is estimated by the Monte Carlo simulation method using 100 independent realizations of random vector P unc following its probability distribution. The optimization problem (defined by Eq. ( 22)) is solved using the trial method for which the admissible set C ad is meshed as follows: for K h1 /m h , 13 values are nonuniformly selected in the interval [703 , 3 830], and for g, 20 nonuniform values in [0 , 0.038]. Thus, 26, 000 numerical simulations have been carried out to solve optimization problem for each level of uncertainties. For computation, the initial value problem, defined by Eqs. [START_REF] Cartmell | Introduction to Linear[END_REF] to [START_REF] Yue | Global analysis of boundary and interior crises in an elastic impact oscillator[END_REF] and rewritten in a dimensionless form, is solved using the time-integration scheme described in Section 4 with duration T = 10.0 s. The values of the fixed parameters are m c = 0.3 Kg, m h = 0.5 Kg, Δ = 0.01 m, and the others are given in Table 1. Upper bound c elec is fixed to the value 6.00 W. For the deterministic case, the components of the optimal solution p opt des are (K h1 /m h ) opt = 1, 580 rad 2 /s 2 and g opt = 0.011 m. For case with uncertainties, for which δ K i is fixed to 0.1, and δ C i to 0.1, we obtain, for δ K h1 = 0, (K h1 /m h ) opt = 957 rad 2 /s 2 and g opt = 0.018 m, for δ K h1 = 0.1, (K h1 /m h ) opt = 1, 950 rad 2 /s 2 and g opt = 0.008 m, and for δ K h1 = 0.4, (K h1 /m h ) opt = 2, 360 rad 2 /s 2 and g opt = 0.008 m. Figures 8 and9 display the graphs of the cost function defined by Eq. ( 21) as a function of the design parameter for these four cases. These figures show that, for each case, the optimal value of the design parameter correspond to a global maximum in C ad . The role played by uncertainties on the optimal values of the design parameters can be analyzed through Figs. 10 and11, which display the graphs g → E{Π imp ((K h1 /m h ) opt , g)}, K h1 /m h → E{Π imp (K h1 /m h , g opt )}, g → E{Π elec ((K h1 /m h ) opt , g)}, and K h1 /m h → E{Π elec (K h1 /m h , g opt )}. These figures show that the optimal design point strongly depends on the level of uncertainties. In particular, it can be deduced that the mean value of the electric power increases with the increase of the gap. The robustness of the optimal design point, p opt des , can be analyzed in studying the evolution of the coefficient variation, δ Π imp (p opt des ), of random variable Π imp (p opt des ) as a function of the uncertainty level. However, in order to better analyze the sensitivity of the responses with respect to the uncertainty level, we have constructed Fig. 12 that displays the graphs g → δ Π imp ((K h1 /m h ) opt , g) and K h1 /m h → δ Π imp (K h1 /m h , g opt ). For each level of uncertainties, it can be seen that the value δ Π imp (p opt des ) occurs in a region for which the two following functions g → δ Π imp ((K h1 /m h ) opt , g) and K h1 /m h → δ Π imp (K h1 /m h , g opt ) are minima. This means the optimal design point is robust with respect to uncertainties.

Conclusions

In this paper, the formulation and the solution of a robust design optimization problem have been presented for a nonlinear electro-mechanical vibro-impact system in presence of uncertainties in the computational model. Since this nonlinear electro-mechanical system is devoted to the vibro-impact optimization, the time responses exhibit numerous shocks that have to be identified with accuracy, and consequently, a very small time step is required. We have thus chosen an explicit time-integration scheme and not an implicit one. Nevertheless, due to the presence of low-frequency contributions in the time responses, a long time duration is required, which will imply a huge number of integration time step if the time step were chosen constant. This is the reason why we have implemented an adaptive integration time step. It was one of the difficulties encountered for the solver implementation. The use a varying time-step integration scheme was not the only strategy adopted to reduce the computation time. The initial value problem has been rewritten in a dimensionless form, which reduced the computation time of each simulation from 8 minutes to 5 minutes on average. Furthermore, a cluster with 20 computers has been used to to parallelize the simulations carried out in the sensitivity analysis and in the optimization problem. Observing the results of numerical integration, as time histories and phase diagrams, some interesting phenomena were verified, as for example bifurcation. Bifurcation is a typical nonlinear phenomena, and it is frequently discussed in many works (see for instance [START_REF] Ragulskis | The effect of dynamical self-orientation and its applicability for identification of natural frequencies[END_REF]). In the present paper, it appears because depending on the values of the system parameters, the system response will have the occurrence or the non occurrence of impacts. But this topic is an ongoing research that will be object of a future work. The construction of the solution for the design optimization problem, has been prepared by carrying out a sensitivity analysis with respect all the possible design parameters. This pre-analysis has allowed for reducing the number of design parameters to two parameters. Consequently, a random search algorithm or a genetic algorithm was not necessary, and we have thus used a trail method. It should be noted that in the framework of a robust analysis formulated in the context of the probability theory, and taking into account the types of nonlinearities in the dynamical system, the Monte Carlo numerical simulation method has been used, and this introduces a significant increase of the numerical cost. The design optimization problem of the dynamical system without uncertainties yields an optimal design point that differs from the nominal values, and which can not be determined, a priori, without solving the design optimization problem. In addition, the robust analysis that has been presented demonstrates the interest that there is to take into account the uncertainties in the computational model. The optimal design point that has been identified in the robust design framework significantly differs from design point obtained with the computational model without uncertainties. For this electro-mechanical system, it has been seen that, the minimum value of the dispersion of the random output occurs in the region of the optimal design parameters, which means that the optimal design point is robust with respect to uncertainties.

Figure 1 :

 1 Figure 1: Sketch of the DC motor.

Figure 2 :

 2 Figure 2: Motor-cart-hammer coupled system. The nonlinear component spring is drawn as a linear spring with constant k h1 and a nonlinear cubic spring with constant k h3 .

Figure 3 :

 3 Figure 3: For the optimal values (m c /m h ) and Δ : (a) graph of π imp as a function of g and k h1 /m h (varying in all its range of values), (b) graph of π imp as a function of g and k h1 /m h (varying in [0.06 , 0.02] and [1 250 , 1 953] respectively).

Figure 4 :Figure 5 :Figure 6 :Figure 7 :

 4567 Figure 4: (a) Graph of π imp as a function of m c /m h with (k h1 /m h ) , g , and Δ . (b) Graph of π imp as a function of k h1 /m h with (m c /m h ) , g , and Δ .

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: (a) Cost function as function of the design parameters for the deterministic case. (b) Cost function as function of the design parameters for the case in which δ Ki = δ Ci = 0.1 and δ K h1 = 0.

Figure 11 :Figure 12 :

 1112 Figure 11: (a) Mean value of the time average of electric power as function of g with (K h1 /m h ) opt . (b) Mean value of the time average of electric power as function of K h1 /m h with g opt . In both graphs, the E{Π elec (p opt des )} is highlighted for each level of uncertainties with markers.

Table 1 :

 1 Values of the system parameters used in simulations.

	Parameter Value	Parameter Value
	m c	0.50 Kg	ν	2.4 V
	r h	0.30 1/m 2	r	0.307 Ω
	c pin	5.00 Ns/m	l	1.88 × 10 -4 H
	c ext	5.00 Ns/m	j m	1.21 × 10 -4 Kg m 2
	ς int	0.05	b m	1.5452 × 10 -4 Nm/(rad/s)
	k i	10 6 N/m	k e	0.0533 V/(rad/s)
	c i	10 3 Ns/m
	value problem defined by Eqs. (

Acknowledgements

The authors acknowledge the support given by FAPERJ, CNPq and CAPES.


Julien Mille
email: julien.mille@liris.cnrs.fr

Sébastien Bougleux
email: sebastien.bougleux@unicaen.frcohen@ceremade.dauphine.fr

Laurent Cohen

MINIMALLY OVERLAPPING PATHS SETS FOR CLOSED CONTOUR EXTRACTION

Keywords: Segmentation, boundary extraction, minimal path, active contour, overlap

Active contours and minimal paths have been extensively studied theoretical tools for image segmentation. The recent geodesically linked active contour model, which basically consists in a set of vertices connected by paths of minimal cost, blend the benefits of both concepts. This makes up a closed piecewise-defined curve, over which an edge or region energy functional can be formulated. As an important shortcoming, the geodesically linked active contour model in its initial formulation does not guarantee to represent a simple curve, consistent with respect to the purpose of segmentation. In this paper, we propose to extract a similarly piecewise-defined curve from a set of possible paths, such that the resulting structure is guaranteed to represent a relevant closed curve. Toward this goal, we introduce a global constraint penalizing excessive overlap between paths.

Introduction

Methods addressing the problem of two-phase segmentation based on energy minimization techniques and variational principles provide a solid mathematical background, and have proven to find suitable solutions in many practical situations. Among them, active contour models consist in deforming an initial curve until it captures the boundary of the target object. Whether they are implemented in an explicit fashion or using level sets, their evolution is usually driven by gradient descent of the Euler-Lagrange equation, which makes them sensitive to local minima, specifically in the presence of noisy images. Consequently, the quality of the resulting segmentation strongly depends on the initial contour position. Several attempts have been made to reduce this sensitivity, including the addition of terms such as the balloon force [START_REF] Cohen | On active contour models and balloons[END_REF] or the use of discrete optimization heuristics such as dynamic programming [START_REF] Amini | Using dynamic programming for solving variational problems in vision[END_REF] or greedy algorithms [START_REF] Williams | A fast algorithm for active contours and curvature estimation[END_REF]. However, these methods still lead to a local minimum of the energy.

To overcome sensitivity to local minima, [START_REF] Cohen | Global minimum for active contour models: a minimal path approach[END_REF] proposed to find a global minimum of the geodesic active contour functional, provided that one or two points of the target object boundary are initially supplied by the user. The resulting global geodesic curve, which can be respectively closed or open, is efficiently derived from the solution of the Eikonal equation obtained with the Fast Marching method [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF][START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF]. Since the control points are fixed and must be located on the target contour, this latter model does not represent a curve which deforms its shape. Moreover, due to the restricted number of these points, the geodesic may fail to capture a relevant contour if the image is too noisy, not enough contrasted, or if the target contour is too lengthy. While several methods concentrate on avoiding this second drawback [START_REF] Benmansour | Fast object segmentation by growing minimal paths from a single point on 2D or 3D images[END_REF]) [START_REF] Benmansour | Tubular structure segmentation based on minimal path method and anisotropic enhancement[END_REF] [START_REF] Kaul | Detection of curves with unknown endpoints using minimal path techniques[END_REF], the geodesically linked active contour model of [START_REF] Mille | Geodesically linked active contours: evolution strategy based on minimal paths[END_REF] allows to overcome the first one. This latter model combines the advantages of geodesics with the ones of greedy algorithms in order to deform a piecewise geodesic curve. Moreover, it is also able to include region-based energies, such as the minimal variance term proposed by [START_REF] Chan | Active contours without edges[END_REF], or even shape prior terms. Whereas this model is relatively robust to local minima, it can fail to construct a valid closed curve, from the initialization step to the end of the evolution.

To overcome this drawback, we design a new energy functional allowing to find a piecewise smooth curve with minimal overlapping. Given several possible relevant paths, subsequently referred to as admissible paths, the key idea of our contribution is to select the combination of paths generating the most relevant contour. In this extent, we introduce an energy functional, combining contour and region terms with a novel overlapping measure.The construction of admissible paths, the overlapping energy as well as the selection of the optimal combination of paths are described in Section 3. The effectiveness of our extended geodesically linked active contour model, whether given initial points are located on the target contour or far from it, is shown in Section 4 through several experiments. The concepts on which relies the proposed approach are recalled in the following section.

Related concepts

Minimal paths

To extract structures in a given image I : D → R d , [START_REF] Cohen | Global minimum for active contour models: a minimal path approach[END_REF] proposed to find curves of minimal length according to an heterogeneous isotropic metric defined from a potential P :

D → R * + .
This potential, which is chosen to take lower values on the structure of interest, allows to measure the length of piecewise smooth curves γ : [0, 1] → D as fol-lows:

L[γ] = 1 0 P(γ(u)) γ (u) du.
(1)

In the context of contour extraction, curves should be located along edges. The potential is thus defined as P = g + w, where g : D → R + is a decreasing function of the gradient magnitude of the image (usually convolved with the derivatives of a Gaussian with given standard deviation σ),

g(x) = 1 1 + ∇(K σ * I)(x) (2)
and w ∈ R * + is a regularizing constant. The target image structure is then extracted by finding a path of minimal length among all paths connecting two given points x 1 and x 2 located on the structure

argmin γ⊂D {L[γ]} s.t. γ(0) = x 1 γ(1) = x 2 . ( 3
)
Such a globally defined minimal path is called a geodesic. The solution of minimization problem (3) can be obtained by considering the geodesic distance map, also referred to as the minimal action map, U v : D → R + which assigns, to each point x ∈ D, the length of the minimal path connecting x to a given point v ∈ D:

U v (x) = inf γ {L[γ]} s.t. γ(0) = v, γ(1) = x. ( 4
)
This map is the unique viscosity solution of the Eikonal equation

∇U v (x) = P(x), ∀x ∈ D \ {v}, U v (v) = 0, (5)
see for instance [START_REF] Crandall | User's guide to viscosity solutions of second order partial differential equations[END_REF]. This allows to replace optimization problem (4) by a partial differential equation. Its discrete version, on a cartesian grid, can be efficiently solved by the Fast Marching (FM) method in O(N log N) operations, where N is the number of grid points [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF][START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF][START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF]. Once the distance map has been numerically computed, the minimal path between v and any other point x of D can be extracted by a gra-

dient descent on U v    γ (u) = - ∇U v (γ(u)) ∇U v (γ(u)) , γ(0) = x, (6)
that corresponds to a back-propagation starting from x until v is reached. In practice, since the FM is a monotonically front propagation method, finding the minimal path bewteen two points does not require to compute the distance on the whole domain D. Starting from one point, the FM can be stopped when the second point is reached, ensuring that the minimal path can be extracted with (6).

The minimal path approach is not restricted to extract an open curve, provided its endpoints. In particular, in the context of object extraction, it is able to find a closed curve, provided only one point on the target object boundary. The closed curve is obtained by detecting a saddle point of the distance map and then by performing two back-propagations, in opposite directions, starting from this saddle point [START_REF] Cohen | Global minimum for active contour models: a minimal path approach[END_REF]. Whether the curve is closed or open, the minimal path approach can fail to extract the desired curve. As depicted in Fig. 1(d), some portions of the minimal path do not follow the desired curve. This happens for instance when P is too noisy or not enough contrasted, when the length of the target curve is too important, or when the regularization constant w is too high. This undesirable behaviour hides a sampling problem, that is one or two points are usually not enough to capture the whole desired curve.

To overcome this drawback, several approaches aim at finding a piecewise geodesic curve Γ = (V, {γ i } i ), where V = {v i } i is a set of vertices that samples the structure to extract, and {γ i } i is the set of geodesics connecting pairs of succesive vertices (see Fig. 1(e)):

γ i = argmin γ {L[γ]} s.t. γ(0) = v i , γ(1) = v i+1 . (7)
Given an initial vertex set, new vertices can be recursively and efficiently detected in several practical situations such that the resulting piecewise geodesic curve matches the desired structure. See [START_REF] Benmansour | Fast object segmentation by growing minimal paths from a single point on 2D or 3D images[END_REF][START_REF] Peyré | Geodesic methods in computer vision and graphics[END_REF] for a complete survey. In Section 3 we propose an alternative approach to overcome this drawback when a closed curve needs to be extracted. This approach is closely related to the geodesically linked active contour de-scribed in the following section.

The geodesically linked active contour model

In order to extract an object from an image, [START_REF] Mille | Geodesically linked active contours: evolution strategy based on minimal paths[END_REF] proposed an active contour model explicitely represented by a closed piecewise geodesic curve, allowing initialization inside the object or around the object boundary. The optimal contour is defined as the piecewise geodesic closed curve Γ = (V, {γ i } 1≤i≤n ) that minimizes a weighted sum of edge-based and region-based energy functionals:

E[Γ] = ω edge E edge [Γ] + ω region E region [Γ]. (8)
The edge-based energy integrates the edge indicator function g (Eq. ( 2)) along the geodesics

E edge (Γ) = 1 |Γ| n ∑ i=1 1 0 g(γ i (u)) γ i (u) du. (9)
In order not to penalize lengthy contours, it is normalized by the euclidean length

|Γ| = n ∑ i=1 1 0 γ i (u) du.
One may note that the edge indicator g is used instead of the potential P so that the euclidean component of the curve length is not taken into account. This ensures that short curves, which could be undesirable shortcuts, are not preferred over longer ones. This edge term may be associated with a balloon term similar to the one introduced in [START_REF] Cohen | On active contour models and balloons[END_REF] in order to increase the capture range of the contour when located far from the object boundary.

In addition, the region-based energy allows to overcome limitations of the edge-based only model, in particular when dealing with noisy, low-constrasted or textured images. [START_REF] Mille | Geodesically linked active contours: evolution strategy based on minimal paths[END_REF] proposed to use a modified version of the two-phase piecewise constant segmentation model developped by [START_REF] Chan | Active contours without edges[END_REF]. Assuming that curve Γ partitions the image into inner region Ω in and outer region Ω out , the region term is expressed as the sum of inner and outer image variances:

E region [Γ] = 1 |Ω in | Ω in I(x) -µ in 2 dx + 1 |Ω out | Ω out I(x) -µ out 2 dx, ( 10
)
where µ in and µ out are average colors in these regions. Following [START_REF] Mille | Narrow band region-based active contours and surfaces for 2D and 3D segmentation[END_REF], a relaxed image homogeneity term focusing on the vicinity of the curve, referred to as narrow band region term, was also addressed as a possible replacement for the previous region term.

Evolution of active contours, whether they are implemented explicitly or in a level-set fashion, is usually performed with gradient descent of the Euler-Lagrange equation. However, in the present case, the energy cannot be differentiated with respect to a given vertex v i . It depends on geodesics linked to v i , which are not expressed in closed form. Hence, the piecewise geodesic structure is evolved thanks to a greedy algorithm similar in principle to the one proposed in [START_REF] Williams | A fast algorithm for active contours and curvature estimation[END_REF], which is discrete by nature and does not imply differentiation. Basically, vertices are moved in local windows in order to minimize the selected energy. Let W N be a normaloriented window of length m centered at vertex v i :

W N (v i ) = v i + kn i k = - m 2 • • m 2
where n i is the inward unit normal vector, estimated by finite difference on corresponding points on geodesics γ i and γ i+1 , respectively. Greedy evolution is performed by moving vertex v i to the position in the window which corresponding geodesically linked contour has the smallest energy E. Let us consider a test position ṽi belonging to the window, and its associated test geodesics γi-1 and γi linking it to the neighbors of v i-1 and v i+1 , respectively. The evolution scheme for vertex v i is formalized by the iteration:

v (t+1) i = argmin ṽi ∈W N v (t) i E( Γ)
where Γ is the tested piecewise geodesic curve:

Γ = {γ 1 , ..., γ i-2 , γi-1 , γi , γ i+1 , ..., γ n } (11)
The behavior of the geodesically linked active contour is depicted in Fig. 2. We can observe that it is able to capture accurately the object boundaries with a reduced number of vertices.

While the geodesically linked active contour model allows to blend the benefits of minimal paths and region-based terms, it turns out to have a significant drawback, as its initial state is not necessarily a simple closed curve. As depicted in Fig. 3(a), this can occur when the initial vertices are unevenly distributed around the target boundary. In this case, geodesics are very likely to gather on particular sides of the target instead of roughly covering the boundary. The reason is that each geodesic is generated independently of the others, such that the obtained piecewisegeodesic curve does not depend on the visiting order of pairs of adjacent vertices. This undesirable phenomenon may occur either as soon as the geodesically linked contour is initialized, or after several evolution steps on a previously well initialized contour.

As in Section 2.1, this problem can be seen as a sampling one. Intuitively, one could think of imposing evenly spaced vertices, as depicted in Fig. 3(b), or adding vertices near the parts of the target boundary which are not covered by the piecewise geodesic curve, like in Fig. 3(c). In the considered context, such sampling criteria are difficult to express, since the target boundary is unknown and applications usually need minimal user interaction. Otherwise, one could think of imposing hard constraints on the overlapping between paths or penalizing paths enclosing a region with excessively small area, but the independent construction of paths prevents such constraints to be implemented. We address this shortcoming in what follows.

Finding the best path set

To overcome the drawbacks of the geodesically linked active contour model, we focus on determin- Assuming that several possible relevant paths linking successive vertices are available, the idea of our contribution is to select the combination of paths generating the most relevant boundary curve. This piecewise smooth closed curve is built by selecting a single path from each set, related to each pair of successive vertices. The relevancy of the generated contour is measured by an energy functional, combining contour and region terms with an overlapping measure. This last term ensures the resulting curve to minimally overlap.

v 1 v 2 v 3 v 1 v 2 v 3 (a) (b) v 1 v 2 v 3 v 4 v 1 v 2 v 3 (c) (d)

Sets of admissible paths

Let V = {v i } 1≤i≤n be a sequence of n given vertices.

Instead of a single geodesic γ i for each pair of successive vertices v i and v i+1 , we consider a set S i of K admissible paths available for this pair, as exemplified in Fig. 3(d):

S i = {γ i, j } 1≤ j≤K
Paths in S i are sorted by cost in ascending order, so that γ i,1 actually corresponds to the minimal path between v i and v i+1 whereas the remaining curves γ i, j , 2 ≤ j ≤ K, are only short paths of increasing cost. Moreover, all paths in S i are constrained to be pairwise disjoint, except at their endpoints, which is for-mulated as the following condition:

γ i, j 1 (u) = γ i, j 2 (v), ∀( j 1 , j 2 ) ∈ [1..K] 2 , j 1 = j 2 , ∀(u, v) ∈]0, 1[ 2 , u = v.
One may notice that the current approach is more constrained than the so-called K shortest paths problem [START_REF] Yen | Finding the K shortest loopless paths in a network[END_REF][START_REF] Eppstein | Finding the k shortest paths[END_REF], which, in its basic formulation, does not impose paths to be disjoint. In the present case, the non-overlap constraint simplifies the generation of several paths. Intuitively, in a graph, the disjoint paths between a pair of vertices can be found by running several instances of the shortest path algorithm, after removal of vertices and incident edges belonging to already found paths.

In our approach, recall that paths are extracted using the Fast Marching method. Hence, the K admissible paths are built by successive deletion of already existing paths from the potential map. Curve γ i,1 is the minimal path between v i and v i+1 in the space endowed by the initial potential P 1 = P. Once the minimal path γ i,1 has been computed, the second admissible path γ i,2 is sought under the constraint that it should not pass through points belonging to γ i,1 . Hence, γ i,2 is not a geodesic in the space induced by potential P, but in the space induced by a modified potential P 2 . The deletion of γ i,1 in the modified potential map is achieved by setting the potential to +∞ at all points of the geodesic. Extending this principle to the construction of the j th admissible path γ i, j as shown in Fig. 4, a recursive definition of potential functions can be written as:

P j (x) = +∞ if x ∈ γ i, j-1 P j-1 (x) otherwise.
This leads to the following recursive definition of the set of admissible paths:

γ i, j = argmin γ 1 0 P j (γ(u)) γ (u) du s.t. γ i, j (0) = v i and γ i, j (1) = v i+1 .
From a practical point of view, one may note that the gradient descent scheme used to built path is continuous, which generates path points with real coordinates, whereas the potential is implemented on a discrete grid. Hence, in order to set P j (x) to +∞, we actually set the 4 integer points surrounding x to +∞.

Paths configuration of minimal cost

The computation of an admissible closed contour consists in selecting one path out of each set S i , such that the contour resulting from the concatenation of selected paths exhibit desirable properties in the image. One of these properties is that the generated contour should be simple, i.e. non-intersecting. In practice, it is reasonable to allow some overlapping between paths. A natural example arises when vertices are located far from the target boundaries, which might cause several admissible paths to have common sections before splitting up. Hence, the non-overlapping condition should be reformulated as a soft constraint. Towards this purpose, we first introduce the overlap measure O between two curves:

O[C 1 , C 2 ] = max 1 |C 1 | 1 0 ψ[C 1 (u), C 2 ] C 1 du, 1 |C 2 | 1 0 ψ[C 2 (u), C 1 ] C 2 du .
(12) It may be considered as the similarity counterpart of the modified Hausdorff distance [START_REF] Dubuisson | A modified Hausdorff distance for object matching[END_REF], as the integrated quantity is a proximity measure instead of a distance. Penalty functional ψ measures the cost of the proximity of point x to curve C . We chose a truncated linear decreasing function of the euclidean distance between x and its nearest point located on C :

ψ[x, C ] = max 0, 1 -α min v∈[0,1] x -C (v)
, where weight α controls the decreasing slope, which is related to the fuzziness of the overlap cost. Note

that O is symmetrical and O[C , C ] = 1.
The computation of an admissible closed contour can be formulated as determining the sequence of labels {x 1 , x 2 , . . . , x n } ∈ [1..K] n minimizing an energy functional E, where label x i corresponds to the chosen path in set S i : min

{x 1 ,x 2 ,••• ,x n }∈[1..K] n E [Γ (γ 1,x 1 , γ 2,x 2 , • • • , γ n,x n )] ,
where Γ (γ 1,x 1 , γ 2,x 2 , ...., γ n,x n ) is the closed contour built by concatenation of paths γ i,x i . It is subsequently shortened to Γ for simplicity. Energy E is the mathematical formulation of required properties of Γ within the image, extending the energy functional ( 8) involved in the geodesically linked active contour model. It is designed to penalize contours exhibiting strongly overlapping sections, poorly fitting to image edges or enclosing regions with high color disparity:

E[Γ] = E overlap [Γ] + ω edge E edge [Γ] + ω region E region [Γ]. (13)
Weights ω edge and ω region are user-defined parameters controlling the relative significance of the edge and region terms over the overlap term. This last one is defined by applying the overlap measure defined in equation ( 12) over all pairs of paths:

E overlap [Γ] = n-1 ∑ i=1 n ∑ j=i+1 O γ i,x i , γ j,x j .
The edge energy integrates the edge indicator function g along paths normalized by their euclidean length:

E edge [Γ] = n ∑ i=1 1 |γ i,x i | 1 0 g(γ i,x i (u)) γ i,x i (u) du.
Unlike the previous edge term in Eq. ( 9), normalization by euclidean length is performed on each path before summation. This makes E edge a separable sum of path-wise terms, which is an advantageous property for optimization (this point is further discussed in subsection 3.3). As the current curve to be optimized is closed, we propose to use a region term, similar to (10), which combines image color variances of the two regions delimited by the curve, as proposed by [START_REF] Chan | Active contours without edges[END_REF]):

E region [Γ] = λ |Ω in | Ω in I(x) -µ in 2 dx + 1 -λ |Ω out | Ω out I(x) -µ out 2 dx, (14)
where λ ∈ [0, 1] controls the blending of the two terms.

While the overlap and the edge energy functionals constitute the building blocks of the proposed model, the region term can be easily replaced in specific situations, e.g. with piecewise-smooth models [START_REF] Lankton | Localizing regionbased active contours[END_REF][START_REF] Brox | On local region models and a statistical interpretation of the piecewise smooth Mumford-Shah functional[END_REF] or texture features [START_REF] Sagiv | Integrated active contours for texture segmentation[END_REF].

Optimization

The best sequence of labels {x 1 , x 2 , ..., x n } is determined using a brute force search among the K n possible configurations. Note that all energy terms are fully or partially precomputed before testing these configurations. Trivially, the edge term needs to be computed only once for each path. Overlap coefficients O[•,•] are pre-computed between all pairs of path and stored in an upper triangular similarity matrix, allowing straightforward computation of the overlap term.

Regarding the region term, Green's theorem enables to convert region integrals over Ω in (Γ) -and simultaneously over Ω out (Γ) = D\Ω in (Γ) -into a sum of contour integrals over each path, according to the following template formula. For any integrable function f over D, we have:

Ω in (Γ) f (x)dx = n ∑ i=1 1 0 F(γ i,x i ) • γ i,x i ⊥ du ( 15
)
where F is a vector field verifying div F = f . The computation of color means and variances is therefore separable over each path, which allows precomputation. Vector field F is obtained by integrating f along the x and y-dimensions:

F(x, y) = x 0 f (t, y)dt , y 0 f (x,t)dt
T Two vector fields are computed once at every location, for f = I and f = I 2 , which allows color means and variances, over inner and outer regions of a given tested configuration Γ, to be efficiently determined.

Experiments and discussion

We demonstrate the ability of the model to recover closed boundaries of objects in natural color images, given few user-provided points. These points are either located on the target boundary, to assess the relevancy of the proposed approach independently from any deformation algorithm, or far from the boundary, in order to show its benefits when integrated into the deformation process of the geodesically linked active contour model. In all experiments, regularization weight w was set to 0.01, the RGB components being assumed to vary from 0 to 1. Low values of w prevent paths from creating undesirable shortcuts, therefore favouring high gradient edges. Reported execution times were obtained with a C++ implementation running on a standard Intel Core2 Duo 2.8GHz architecture with 4Gb RAM.

From points localized on the contour

Fig. 5 depicts an experiment intended to demonstrate the consistency of contour extraction with respect to various initial locations of vertices. The closed contour was generated given n = 2 vertices and K = 5 admissible paths. On the aforementioned architecture, paths generation took 0.94s and contour selection took 0.14s. As neighboring areas are sufficiently contrasted, the edge map alone turned out to be reliable enough, which allowed not to use the region homogeneity term. Hence, only the overlap and edge terms were used on this particular image.

The proposed approach proves to recover suitable contours regardless of the positions of endpoints. In particular, vertices do not need to be evenly spaced along the actual object boundary.

One of the main benefits of the proposed approach over classical minimal path-based segmentation is the ability to formulate a region-based criterion, as in classical active contours. Figures 6 and7 illustrate the interest of using such criterion, as well as the overlap constraint. The 700 × 529 data in Fig. 6 was processed with n = 3 vertices and K = 5 admissible paths per pair of successive vertices. Paths generation took 13.3s and contour selection 3.9s. The annoying overlapping phenomenon yielded by the basic geodesically linked active contour is shown in Fig. 6(c). Incorporation of the overlap constraint enables to generate a closed contour, which remains nevertheless unrelevant with respect to image partition. This can be explained by the fact that various edges, stronger that the actual boundaries of the target object, can be encountered in neighboring structures. This undesirable feature is addressed by the addition of the region homogeneity term (weight λ balancing the inner and outer region terms was set to 1, as only the inner object is almost homogeneous).

The last experiment shown in Fig. 7 follows a similar principle, with n = 4 and K = 5, on a 800 × 600 image. Reported execution times for paths generation and contour selection are 22.5s and 7.8s, respectively.

For this image, while the proposed approach outperforms the geodesically linked active contour model at its initialization, the resulting curve fails to fully extract the desired contour, due to weak edges near the object center. Since the object of interest is textured, this could be improved by considering a region energy that favors homogeneous textured regions.

(a2) (b1) (b2) (c1) (c2) (d1) (d2)

Integration into deformation process

We report experiments where the proposed approach was integrated into the evolution algorithm of the geodesically linked active contour model, so that con-tour extraction can be performed with initial points located far from the target contour. The integration is as follows: during deformation, when displacing a given vertex v i , tested geodesics γi-1 and γi are built such that they do not overlap with existing geodesics {γ 1 , ..., γ i-2 , γ i+1 , ..., γ n }. This is achieved by deleting these existing geodesics in the potential map, as described in subsection 3.1. The energy of the tested contour is endowed with the overlap term according to Eq. ( 13).

The proposed approach can solve the overlapping problem arising in two different cases. The first case, shown in Fig. 8(a), corresponds to an overlapping present as soon as the contour is initialized and propagated afterwards. On the other hand, the second case, depicted in Fig. 8(c), shows the result of an overlapping occuring during evolution of a well-initialized curve. In both cases, the integration of our approach when updating geodesics during evolution allows to maintain a valid contour (Fig. 8(b) and Fig. 8(d)), at the expense of additional time cost to check paths configurations.

Conclusion and perspectives

By searching the best paths configurations among sets of admissible paths, given an energy functional combining an edge fitting term, a region homogeneity term and a novel overlap-penalizing energy, we aimed at overcoming some important shortcomings arising in geodesic-based segmentation. The introduced constraints allowed to guarantee consistent closed contours, whether given initial points were located on the target boundaries or far from them. Incorporation into the geodesically linked active contour model demonstrated the advantages of the approach. Future work may focus on designing finer search methods to determine the optimal set of paths, since a basic brute force search was implemented so far. A related possible investigation deals with the generation of admissible paths. In this extent, instead of generating all admissible paths per pair of successive vertices at initialization, one could think of an adaptive approach in which only necessary extra admissible paths would be created during the search process.

Figure 1 :

 1 Figure 1: (a) Input image. (b) Potential P. (c) Geodesic between two given points. (d) Undesirable sructure extraction. (e) Piecewise geodesic curve.

Figure 2 :

 2 Figure 2: Evolution of geodesically linked active contour: in the evolution steps, all test geodesics from neighboring vertices to test positions in windows are represented.

Figure 3 :

 3 Figure 3: Towards a relevant initialization of the geodesically linked active contour model: (a) undesirable overlapping with unevenly spaced vertices, (b) improvement by even spacing of vertices, (c) improvement by addition of vertex, (d) admissible paths sets between pairs of vertices with K = 4 paths per pair and high regularization constant w

Figure 4 :

 4 Figure 4: Successive potential maps P j (top row) and corresponding admissible paths γ •, j (bottom row) given two endpoints.

Figure 5 :Figure 6 :

 56 Figure 5: Robustness of contour extraction with respect to vertices locations: all admissible paths (top row) and selected contour (bottom row).

Figure 7 :

 7 Figure 7: Influence of the region homogeneity term: (a) input image, (b) inverted gradient magnitude, (c) initial configuration of the basic geodesically linked active contour (independent minimal paths), (d) all admissible paths, (e) selected contour with overlap and edge terms, (f) selected contour with overlap, edge and region terms

Figure 8 :

 8 Figure 8: Integration of non-overlapping constraint in the geodesically linked active contour model. For each subfigure: (x1) initial state, (x2) final state. (a & c) Basic geodesically linked active contour, (b & d) Geodesically linked active contour model with non-overlapping constraint

a a This work was supported by a grant from the GdR ISIS (french research group on Information, Signal, Image and Vision) as well as the REVES project (funded by the French National Agency for Research with the reference ANR-10-CORD-005)


Alan Chatham

Amani Naseem

Allan Fowler
email: allan.fowler@waiariki.ac.nz

Regina Bernhaupt
email: regina.bernhaupt@irit.fr

Ben Schouten
email: b.a.m.schouten@tue.nl

Florian ' Floyd' Mueller

Menno Deen
email: menno.deen@fontys.nl

Robert Cercos
email: robert@exertiongameslab.org

Keywords: Game Jams, Design, Experimentation, Researchthrough-design

Recent years have witnessed a rise in Game Jamsorganized events to create playable prototypes in a very short time frame. Game Jams offer a unique and quick way to prototype games. Beyond that, we believe Game Jams can also be seen as a design research method, situated in the research-through-design tradition, to create knowledge in a fast-paced, collaborative environment. The goal of this Game Jam is thus twofold: first, participants will use the Game Jam approach to investigate a research question; second, participants can, through actual practice, identify advantages and disadvantages of Game Jams as a research method. Hereby the Game Jam workshop provides a unique opportunity for HCI practitioners and researchers to gain experience in applying gameoriented methods for research.

Introduction

In 2013, a Game Jam workshop introduced the CHI community to a collaborative, method of designing, as currently successfully practiced in game design communities [START_REF] Cowan | Ubiquitous sketching for social media[END_REF]. Inspired by this workshop's reflective session, this year's jam examines the opportunity of utilizing Game Jams as a design research tool. The CHI 2014 Game Jam Workshop invites researchers and practitioners to design a game where the design process and resulting game artifact can serve as a case study. This method will contribute to existing HCI methods, by generating knowledge that benefits from being derived from a range of case studies, practical craft knowledge and collaborative reflection. Furthermore, the method sits within the existing HCI research-through-design tradition, helping the researcher to approach the design of the game both as a researcher and a game designer [START_REF] Preston | Game Jams: Community, Motivations, and Learning among Jammers[END_REF]. As such, we believe Game Jams as a research tool provides opportunities to supplement existing HCI methods [START_REF] Gaydos | Game Jams: Designing Design Spaces[END_REF]. This workshop therefore gives HCI researchers and practitioners an opportunity to explore the value of Game Jams as research tool by practically engaging with the game design process. What's more, jamming may broaden and deepen social interactions [START_REF] Olberding | AugmentedForearm: exploring the design space of a display-enhanced forearm[END_REF] amongst participants.

Background

Recent years have witnessed a rise in Game Jam events. Participants from diverse backgrounds (art, programming, design and scholars) meet to create functional games in a short, predefined time frame. Game Jams have proven to be successful methods for rapidly prototyping innovative games [START_REF] Fowler | The evolution and significance of the Global Game Jam[END_REF]and have educational value for jammers [START_REF] Musil | Synthesized essence: What Game Jams teach prototyping of new software products[END_REF], [START_REF] Preston | Game Jams: Community, Motivations, and Learning among Jammers[END_REF], [START_REF] Zook | Game Conceptualization and Development Processes in the Global Game Jam[END_REF]. As Game Jams are time bound and unconstrained from commercial concerns, they proved a unique method to push the limits of game design and development. For example, the Global Game Jam (GGJ) attracts thousands of developers, and is growing every year. Based on the continued growth of this event, we predict an even larger crowd for the next GGJ [7]. Furthermore, some games developed during the GGJ, amongst others Spelunky and Surgeon Simulator 2013, are now commercially available. Other Game Jams such as the MolyJam, LudumDare and 7FPS add to the number of developed games. In four years combined, these Game Jams produced a total of approximately 17,000 games. [8] shows how a Game Jam can be utilized as a method to explore a specific research question. Initiated by academic research programs, these Game Jams study ways to elicit physical activity amongst youngsters in different situational contexts by a research-through-design method. Explorations in game development for public spaces and swimming pools presented researchers with insights into the design processes of said games. Hereby the design guidelines for games that elicit physical activity were deduced from both theoretical resources as well as actual design experiences.

The annual Games [4Health] Jam

The CHI 2013 Game Jam resulted in a range of divergent games. It also resulted in an interesting reflection exercise at the conclusion of the workshop where participants discussed the possibility of bringing this format into their research labs as a method to accelerate their progress. This year's workshop is an answer to this reflection: proposing to utilize the established rapid-pace format of Game Jams as a research tool to reveal knowledge and create new understandings. The findings will mainly concern design knowledge, but are not limited to this As such, the Game Jam [4Research] the opportunity of utilizing Game Jams as a design research tool. During the Game Jam, researchers team up in an attempt to gain insight into research question by creating a game. These games or design processes will act as case studies that create insights into the questions posed. For example, explorations can be done in making Augmented Forearm interfaces [START_REF] Olberding | AugmentedForearm: exploring the design space of a display-enhanced forearm[END_REF] more playful, or how games like DrawSome or Pictionary can benefit from digital pen input.

The collaborative nature of the Game Jam provides an opportunity for participants from varied backgrounds to network and bridge gaps that exist between institution, research, and industry. Musil et al. [START_REF] Musil | Synthesized essence: What Game Jams teach prototyping of new software products[END_REF] suggest that Game Jams provide an effective and focused experience and that participants gain valuable skills in prototyping and collaboration. Preston et al. [START_REF] Preston | Game Jams: Community, Motivations, and Learning among Jammers[END_REF] advocate the use of Game Jams to nurture innovation, and computational and rule-based thinking. By bringing together developers and researchers in a Game Jam we may build research ventures, networks, and long-term collaborations between participants.

Goals •

Explore jamming as research method • Learn this method by doing it • Identify the advantages of this method (rapid, collaborative)

•

Identify the disadvantages (difficult to generalize results, early bugs in the prototype that may hinder any study) • Develop a playable game

Workshop plan

Participants are encouraged to think about tools and methods to bring to the workshop in order to build games. The workshop will take place over two days.

During the first day, participants will introduce themselves, and pitch ideas. Then teams will be formed and work will begin.

After lunchtime jammers discuss which existing research method would best benefit the research question posed and how the Game Jam can contribute to this method. At the end of the first day, we will close with short presentations from each group, sharing ideas and techniques, followed by an informal dinner to encourage cross-group discussion and collaboration.

The second day will begin with continued work on projects after a short reflection of the day before: "what are the pros and cons of a Game Jam as research tool", and "how can the Game Jam format be improved to contribute insights?" In the afternoon an informal presentations among the participants is scheduled. The workshop will close by bringing the games into the common areas of the conference to expose participants' work to the broader community. Special attention is given by organizers to elicit the 'Game-Jam-feeling' of belonging, which appears to foster co-operation and social negotiations amongst jammers.

Follow-up and dissemination

We will document the attendees' processes and results with photographs and video. An example of this can be found on: http://youtu.be/bcPHa4NfC1tM. We aim to write an article about lessons learned, involving everyone by means of the format of massively collaborative editing [START_REF] Tomlinson | Massively distributed authorship of academic papers[END_REF].


V Robert

D Pascu

J E Arlot

V Lainey

A new reduction of USNO photographic plates of the Martian satellites

Thanks to the new technologies, a new astrometric reduction of old photographic plates can provide a better knowledge of the orbital motion of planetary satellites. In the cadre of the FP7 European project, USNO plates were digitized with the new generation DAMIAN scanning machine. The procedure was applied to various photographic plates and in particular to the USNO photographic plates of the Galilean satellites. The astrometric results are the most precise.

Here we consider a set of a few hundred photographic plates of the Martian satellites, taken at the USNO, and covering the years 1969-1997. A specific procedure was developed to obtain a high precision and we expect an accuracy better than 100 mas in (RA,Dec) positions of each moon. Since the position of Mars may also be deduced from the observed RA and Dec positions of Phobos and Deimos, we can also assess the accuracy of Mars' ephemeris. First astrometric results will be presented.

Extraction and identification

The positions on the images of the stars are extracted by means of the Source Extractor software [START_REF] Bertin | SExtractor : software for source extraction[END_REF], used to create a list of objects detected on the plates. A selection for about circular objects and FWHM and magnitude constraints permit to reduce the list of possible detected sources. The positions of the planet are extracted by an IDL ellipse fitting, the positions of the satellites are extracted by an IDL Gaussian fitting after correcting for the background gradient due to the filters' halo. The objects assumed to be present in the field are identified from existing catalogs to select only real astrometric sources. This process defines areas on the images in which the objects have to be found (Robert, 2011). All the available stars (depending on the catalog used) are identified and more, those that are not visible with the eyes. The method can also be applied with all planetary systems; tests were successfully performed with Jupiter and Saturn USNO images, with Pluto OHP and Bucarest images. The obtained positions are finally corrected for the optical distortion introduced by the objective/camera unit during the digitization and for the instrumental coma-magnitude effect.

Astrometric results

The (O-C)s are for observations of Mars, Phobos and Deimos. Positions are deduced from the measurements. Figures show details of the (RA,Dec) (O-C)s according to INPOP10 [START_REF] Fienga | INPOP10a[END_REF] planetary ephemeris and IMCCE [START_REF] Lainey | First numerical ephemerides of the Martian moons[END_REF]

ESPaCE FP7 program and USNO Martian plates

The European Satellite PArtnership for Computing Ephemerides FP7 program activity is focused on the extraction and analysis of astrometric data from space measurements not yet applied to the dynamics and to combine them with ground-based astrometric data. Part of this project aims at making accurate positional measurements of planets or satellites. The resulting astrometric data will be used to model the motions of the bodies through new SPICE kernels. We focus here on a subset of a few hundreds photographic plates of the Martian satellites, taken with the USNO 61-inch reflector and the USNO 26-inch refractor by D. [START_REF] Pascu | Astrometric techniques for the observation of planetary satellites[END_REF][START_REF] Pascu | The Naval Observatory Program for the Astrometric Observation of Planetary Satellites[END_REF][START_REF] Pascu | An appraisal of the USNO program for photographic astrometry of bright planetary satellites[END_REF] on Kodak 103a-type from 1967 to 1997. Each plate contains from 2 to 3 observations; the exposures are shifted on the RA direction. The exposure time is about 40 sec; the field is 57' on the X-axis. All these photographic plates were recently digitized with the new generation DAMIAN digitizer (Robert et al., 2012). We present here the first results for observations from 1971 to 1997.

Astrometric reduction

Astrometric (RA,Dec) results of the Martian satellites and the planet are determined in an ICRS geocentric frame to be easily compared with the most recent satellite and planetary ephemerides. Because of the few number of available stars, the astrometric reduction is quite different from a common astrometric process: the positions are here calculated after correcting for instrumental and spherical effects that take into account the parallax and aberration effects and the total atmospheric refraction. (RA,Dec) positions are deduced from tangential apparent coordinates; the following equation system present the (x,y) → (X,Y) transformation in the case of the 26-inch observations:

Conclusion

We demonstrated the high interest to continue the analysis of old photographic plates such as USNO's. Thanks to the new technologies, we were able to provide astrometric data with a high accuracy after applying the necessary corrections. We now provide an accuracy about of 70 mas for (RA,Dec) positions of the planet and its moons. Note that the previous accuracy was about of 200 mas and only for relative positions. More important, these results indicate that we now can reach an accuracy better than that of CCD observations [START_REF] Colas | Comparison of observations of the Martian satellites made in 1988 with ephemerides[END_REF] and as good as the old spacecraft measurements of the system. These results encourage us to continue the analysis of old photograhic plates. We have at the present time a set of several hundred USNO plates of the Saturnian system. The digitization and new reduction are also part of the FP7 European project. The step after will be to reduce other relevant old photographic plates such as Yale Southern Station Saturnian observations.

  satellite ephemerides. The first and second tables show details of the (RA,Dec) (O-C)s and rms residuals in mas according to INPOP10 and DE423 (JPL) planetary ephemerides respectively, and IMCCE satellite ephemerides in both cases. INPOP10 and DE423 ephemerides can thus be compared. The first and third tables show details of the (RA,Dec) (O-C)s and rms residuals in mas according to INPOP10 planetary ephemeris in both cases, and IMCCE and mar097 (JPL) satellite ephemerides respectively. The two Martian ephemerides can thus be compared.

  𝑋 = 𝜌 cos 𝜃 𝑥 -(𝜌 + 𝜖 1 sin(𝜖 2 𝑡 𝑚 + 𝜖 3 )) sin 𝜃 𝑦 + Δ 𝑥 + 𝐶 𝑥 𝑥(𝑚 -𝑚 0 ) 𝑌 = 𝜌 sin 𝜃 𝑥 + (𝜌 + 𝜖 1 sin(𝜖 2 𝑡 𝑚 + 𝜖 3 )) cos 𝜃 𝑦 + Δ 𝑦 + 𝐶 𝑦 𝑦(𝑚 -𝑚 0 ) Only 4 parameters modeling the scale, orientation and center of field are fitted for a minimum of 2 reference stars. The contribution of each effect is separated. The USNO plate n°11001 Center of the USNO plate n°11001 The USNO plate n°11014 Center of the USNO plate n°11014


V Robert

D Pascu

V Lainey

J E Arlot

Astrometry of natural satellites Improving the dynamics of planetary systems with old observations

A new astrometric reduction of old photographic plates, benefiting from modern technologies such as submicrometric scanners associated with a reduction using accurate catalogues (UCAC at the present time and Gaia in a near future), provides improved knowledge of the orbital motion of planetary satellites. In the framework of an international collaboration first, and in the FP7 ESPaCE european project afterward, U.S. Naval Observatory plates were digitized with the new generation DAMIAN scanning machine of the Royal Observatory of Belgium. The procedure was applied to a few hundred photographic plates of the Galilean satellites covering the years 1967-1998, and of the Martian satellites covering the years 1967-1997.

We provide results with an accuracy better than 70 mas in (RA,Dec) positions of the Galilean moons, and better than 60 mas in (RA,Dec) positions of the Martian satellites. Since the positions of Jupiter and Mars may be deduced from the observed (RA,Dec) positions of their satellites, we can also assess the accuracy of the ephemerides of Jupiter and Mars.

Extraction and identification

The positions of the images of the planet, satellites and stars are extracted by means of the Source Extractor software (Bertin et al., 1996), used to create a list of objects detected on the plates. A selection for approximately circular objects and FWHM and magnitude constraints permit to reduce the list of possible detected sources. More particular, the positions of the planet Mars are extracted by an IDL ellipse and limb fitting, the positions of the satellites are extracted by an IDL Gaussian fitting after correcting for the background gradient due to the planet's halo. The objects assumed to be present in the field are identified from existing catalogs to select only real astrometric sources. This process defines areas on the images in which the objects have to be found (Robert, 2011). All the available stars (depending on the catalog used) are identified and more, those that are not detected by eye. We demonstrated that the method can be applied with all planetary systems. The obtained positions are finally corrected for the optical distortion introduced by the objective/camera unit during the digitization and for the instrumental coma-magnitude effect.

Galilean astrometric results

The (O-C)s are for observations of Io, Europa, Ganymede and Callisto first, and for Mars, Phobos and Deimos following. Positions are deduced from the measurements. The two first Figures show details of the (RA,Dec) (O-C)s according to the INPOP10 [START_REF] Fienga | IMCCE Memorandum[END_REF] planetary ephemeris and IMCCE [START_REF] Lainey | [END_REF](Lainey et al., , 2009) ) [START_REF] Folkner | [END_REF][START_REF] Folkner | [END_REF])) planetary ephemerides.

ESPaCE FP7 program and USNO photographic plates

The European Satellites PArtnership for Computing Ephemerides ESPaCE project aims at strengthening the collaboration and at developing new knowledge, new technology, and products for the scientific community in the domains of the development of ephemerides and reference systems for natural satellites and spacecraft. The main European research centers involved in space sciences and dynamics contribute by combining their expertise. Their activities are focused on the extraction and analysis of astrometric data from space measurements not yet applied to the dynamics and to combine them with ground-based astrometric data. We focus here on a subset of a few hundred photographic plates of the Galilean and Martian satellites, taken with the USNO 61-inch reflector (37' field) and the USNO 26-inch refractor (57' field) by D. Pascu (1977[START_REF] Pascu | Natural and Artificial Satellite Motion[END_REF][START_REF] Pascu | Galactic and Solar System Optical Astrometry[END_REF] from 1967 to 1998. Each plate contains several exposures shifted on the Dec axis for the Galileans and on the RA axis for the Martians. All these photographic plates were recently digitized with the new generation DAMIAN digitizer (Robert et al., 2011).

Astrometric reduction

Astrometric (RA,Dec) results of the satellites and the planets are determined in an ICRS geocentric frame to be easily compared with the most recent satellite and planetary ephemerides. Because of the few number of available stars, the astrometric reduction is quite different from a common astrometric process: the positions are calculated after correcting for instrumental and spherical effects that take into account the parallax and aberration effects and the total atmospheric refraction. (RA,Dec) positions are deduced from tangential apparent coordinates; only 4 parameters modeling the scale, orientation and center of field are fitted for a minimum of 2 reference stars. The contribution of each effect is separated.

Galilean

Conclusion

We demonstrated the high interest to continue the analysis of old photographic plates such as USNO's. Thanks to the new technologies, we were able to provide astrometric data with a high accuracy after applying the necessary corrections. We now provide an accuracy about of 70 mas (~210 km) for (RA,Dec) positions of the planet Jupiter and its moons. We also provide an accuracy about of 60 mas (~18 km) for (RA,Dec) positions of the planet Mars and its satellites. Note that the previous accuracy was about of 200 mas and only for relative positions. More important, these results indicate that we now can reach an accuracy better than that of CCD observations (Colas et al., 1991) and as good as the old spacecraft measurements of the Martian system for example. These results encourage us to continue the analysis of old photograhic plates. We have at the present time a set of several hundred USNO plates of the Saturnian system. The digitization and new reduction are also part of the FP7 European project. The steps after will be to reduce other relevant old photographic plates such as Yale Southern Station Saturnian observations, and to continue the effort in reducing old observations, participating in new programs such as the NAROO project (New Astrometric Reduction of Old Observations) of the IMCCE / Paris Observatory.

  satellite ephemerides. The next Figure shows the contribution of the temperature on the scale variation, while the last Figure indicates that the planetary halo has negligible effect on the satellite positions. Tables compare the averages of the Galilean and Martian (RA,Dec) (O-C)s and rms residuals respectively, in mas, according to the INPOP08 (Fienga et al., 2008), INPOP10, DE421 and DE430

The USNO Galilean plate n°24014


V Robert

D Pascu

V Lainey

J E Arlot

Astrometry of natural satellites Improving the dynamics of Jupiter, Mars and Saturn systems with old observations

A new astrometric reduction of old photographic plates, benefiting from modern technologies such as submicrometric scanners associated with a reduction using accurate catalogues (UCAC at the present time and Gaia in a near future), provides improved knowledge of the orbital motion of planetary satellites. In the framework of an international collaboration first, and in the FP7 ESPaCE European project afterward, U.S. Naval Observatory plates were digitized with a new generation scanning machine of the Royal Observatory of Belgium. The procedure was applied to a few hundred photographic plates of the Galilean satellites covering the years 1967-1998, of the Martian satellites covering the years 1967-1997, and of the major Saturnian satellites covering the years 1974-1998. We provide results with an accuracy better than 70 mas in (RA,Dec) positions of the Galilean moons, better than 60 mas in (RA,Dec) positions of the Martian satellites, and better than 100 mas in (RA,Dec) positions of the Saturnian moons. Since the positions of the planets may be deduced from the observed (RA,Dec) positions of their satellites, we can also assess the accuracy of the ephemerides of Jupiter, Mars and Saturn.

Analysis and astrometric reduction

The positions of the images of the planet, satellites and stars are extracted by means of a specific process including the use of the Source Extractor software (Bertin et al., 1996) and an IDL ellipse and limb fitting. The objects assumed to be present in the field are identified from existing catalogs to select only real astrometric sources. This process defines areas on the images in which the objects have to be found (Robert, 2011). All the available stars (depending on the catalog used) are identified and more, those that are not detected by eye. The obtained positions are finally corrected for the optical distortion introduced by the objective/camera unit during the digitization and for the instrumental coma-magnitude effect. Astrometric (RA,Dec) results of the satellites and the planets are determined in an ICRS geocentric frame to be easily compared with the most recent satellite and planetary ephemerides. Because of the few number of available stars, the astrometric reduction is quite different from a common astrometric process: the positions are calculated after correcting for instrumental and spherical effects that take into account the parallax and aberration effects and the total atmospheric refraction. (RA,Dec) positions are deduced from tangential apparent coordinates; only 4 parameters modeling the scale, orientation and center of field are fitted for a minimum of 2 reference stars. The contribution of each effect is separated.

ESPaCE FP7 program and USNO photographic plates

The European Satellites Partnership for Computing Ephemerides ESPaCE project aims at strengthening the collaboration and at developing new knowledge, new technology, and products for the scientific community in the domains of the development of ephemerides and reference systems for natural satellites and spacecraft. The main European research centers involved in space sciences and dynamics contribute by combining their expertise. Their activities are focused on the extraction and analysis of astrometric data from space measurements not yet applied to the dynamics and to combine them with ground-based astrometric data. We focus here on a subset of a few hundred photographic plates of the Galilean, Martian and Saturnian satellites, taken with the USNO 61-inch reflector (37' field) and the USNO 26-inch refractor (57' field) by D. [START_REF] Pascu | [END_REF][START_REF] Pascu | Natural and Artificial Satellite Motion[END_REF][START_REF] Pascu | Galactic and Solar System Optical Astrometry[END_REF] from 1967 to 1998. Each plate contains several exposures shifted either on the RA or the Dec axis. All these photographic plates were recently digitized with the new generation ROB digitizer (Robert et al., 2011(Robert et al., , 2014)).

Galilean obs.

(

Conclusion

We demonstrated the high interest to continue the analysis of old photographic plates such as USNO's. Thanks to the new technologies, we were able to provide astrometric data with a high accuracy after applying the necessary corrections. We now provide an accuracy about of 70 mas (~210 km) for (RA,Dec) positions of the planet Jupiter and its moons, 60 mas (~18 km) for (RA,Dec) positions of the planet Mars and its satellites, and 100 mas (~700 km) for (RA,Dec) positions of the planet Saturn and its satellites. Note that the previous accuracy was up to 200 mas and only for relative positions. More important, these results indicate that we now can reach an accuracy better than that of CCD observations (Colas et al., 1991) and as good as the old spacecraft measurements of the Martian system for example. These results encourage us to continue the analysis of old photograhic plates. The steps after will be to reduce other relevant old photographic plates such as Yale Southern Station Saturnian observations, and to continue the effort in reducing old observations. Old CCD observations could also be re-reduced in the same scope, increasing the accuracy of these data.

To fullfill our program, we have started the NAROO project (New Astrometric Reduction of Old Observations) of the IMCCE -Paris Observatory, A new sub-micrometric machine will be set up in France to digitize several thousands photographic plates for astrometry and thus dynamics purposes. Moreover, the Gaia astrometric catalog of reference stars will provide proper motions of sources up to mag 18 with an accuracy better than 6 mas over one century. We look forward its arrival which will be be a revolution in Solar System astrometry: reductions of old observations will yield increase accuracy by eliminating errors due to "old" reference star catalogs. We will be able to observe in the past with today accuracy, that is essential for fast moving objects of the Solar System.

Martian

Astrometric results

The

  (O-C)s are for observations of the Galilean satellites first, then for Mars, Phobos and Deimos, and for the eight major Saturnian satellites. Positions are deduced from the measurements. The plots show details of the (RA,Dec) (O-C)s according to the INPOP10 (Fienga et al., 2010) planetary ephemeris and IMCCE (Lainey et al., 2007, 2009, 2011) satellite ephemerides. The last Figure shows the submicrometric ROB digitizer used for the digitization. Tables compare the averages of the Galilean, Martian and Saturnian (RA,Dec) (O-C)s and rms residuals respectively, in mas, according to the INPOP08 (Fienga et al., 2008), INPOP10, DE421 and DE430 (Folkner et al., 2008, 2014)) planetary ephemerides.

th DPS meeting, Tucson, Arizona, 2014 The research leading to these results has received funding from the European Community's Seventh Framework Programme (FP7/2012-2013) under grant agreement n°263466.