Haoyu Li

Stéphane Derrode
email: stephane.derrode@ec-lyon.fr

Lamia Benyoussef

Wojciech Pieczynski
email: wojciech.pieczynski@telecom-sudparis.eu

Free-walking 3D Pedestrian Large Trajectory Reconstruction from IMU Sensors

Keywords: Pedestrian navigation, Inertial sensor, Hidden Markov model, Kalman filter, Stance detection

This paper presents a pedestrian navigation algorithm based on a foot-mounted 9-DOF Inertial Measurement Unit (IMU), which provides tri-axial accelerations, angular rates and magnetics. Most of algorithms used worldwide employ Zero Velocity Update (ZUPT) to reduce the tremendous error of integration from acceleration to displacement. The crucial part in ZUPT is to detect stance phase precisely, a cyclic left-to-right style Hidden Markov Model is introduced in this work which is able to appropriately model the periodic nature of signals.

Stance detection is then made unsupervised by using a suited learning algorithm. Then orientation estimation is performed independently by a quaternion-based method, a simplified errorstate Extended Kalman Filter (EKF) assists trajectory reconstruction in 3D space, neither extra method nor prior knowledge is needed to estimate the height. Experimental results on large free-walking trajectories show that the proposed algorithm can provide more accurate locations, especially in z-axis compared to competitive algorithms, w.r.t. to a ground-truth obtained using OpenStreetMap.

I. INTRODUCTION

In recent years, Pedestrian Navigation System (PNS) using a foot-mounted IMU sensor has been investigated extensively. IMUs have great advantages as they are small and can be worn on the body. It measures motion kinematic information such as acceleration, angular rate and magnetics, so that it is theoretically possible to transfer the signals from the sensor frame to the earth frame, also called "global frame", based on the sensor orientation, and then to compute the velocity and displacement of motion. Unfortunately, a simple numerical integration of acceleration to obtain the displacement suffers from a tremendous error, which is mainly caused by accelerometer noise and numerical integration. Therefore, ZUPT was then developed and showed great performances in reducing the cumulative double integration errors. The idea behind ZUPT is the following: when a person is walking, his feet alternatively swing in the air and step on the ground, so that his foot velocity is zero when it is attached to the ground. The period during which the velocity is zero is called "stance phase". Foxlin [START_REF] Foxlin | Pedestrian tracking with shoe-mounted inertial sensors[END_REF] firstly proposed an error-state EKFaided ZUPT algorithm for PNS, which exploits the periodical stance phases to reduce both velocity and displacement errors for a much better precision of position. From this work, numerous strategies were then proposed that mainly enrich the measurement vector of EKF, by using for examples angular rate [START_REF]Personal dead reckoning system with shoe mounted inertial sensors[END_REF] or digital compass [START_REF] Norrdine | Step detection for ZUPTaided inertial pedestrian navigation system using foot-mounted permanent magnet[END_REF].

State-of-the-art algorithms on using an EKF-aided ZUPT algorithm allow to precisely reconstruct the trajectory in X-Y plane, also called horizontal plane, but fail in reconstructing the z-axis as they have huge bias in estimating the height. So getting reliable reconstructed trajectory in 3D space requires extra method for robust height estimation. Y. Hsu et al. [START_REF] Hsu | A wearable inertial pedestrian navigation system with quaternion-based ekf for pedestrian localization[END_REF] used a Probabilistic Neural Network (PNN) classifier to identify and count the number of steps climbed. The height can then be estimated assuming that the height of each stair is known, which is limited. Y. Liu et al. [START_REF] Liu | Height error correction based on shoe-mounted pedestrian navigation systems[END_REF] proposed a floor identification method for height estimation, based on a prior knowledge of buildings. SK Park et al. [START_REF] Park | Height compensation using ground inclination estimation in inertial sensor-based pedestrian navigation[END_REF] proposed a height compensation algorithm by estimating road inclination angle. Also, we can also cite fusion algorithms of inertial sensor data with extra sensors (barometer: [START_REF] Shin | Motion recognition-based 3D pedestrian navigation system using smartphone[END_REF], [START_REF] Kim | Performance improvement and height estimation of pedestrian dead-reckoning system using a low-cost MEMS sensor[END_REF], and WI-FI: [START_REF] Chai | 2D/3D indoor navigation based on multi-sensor assisted pedestrian navigation in Wi-Fi environments[END_REF]). Hence, a large majority of works depend on some prior knowledge of the environment or on some additional sensors to assist height estimation.

This paper proposes a new EKF-aided ZUPT based 3D PNS algorithm with only one 9-DOF IMU sensor. By introducing a decomposed Gradient Descent Algorithm (GDA) based quaternion estimation method, the walking orientation of sensor is computed independently from the displacement. Stance phases are then detected by a cyclic Left-to-Right Hidden Markov Model (LR-HMM), so that the EKF-aided ZUPT allows to compensate for the integration error. The precision of the algorithm has been evaluated through long walks down the valley. The remaining of the paper is organized as follows. Section II describes the sensor used and the structure of the proposed PNS algorithm. Section III then explains how the independent quaternion estimation method is developed and works within the PNS. Section IV describes the LR-HMM based stance detection method, and a simplified error-state EFK is introduced. Results of experiments are proposed in Section V.

II. PEDESTRIAN NAVIGATION SYSTEM

The PNS is based on one 9-DOF IMU sensor from Shimmer Sensing 1 with 100Hz sampling rate. The scale of ac- celerometers, gyroscopes and magnetometers are set to ±8g, 1000deg/s and 2.5Ga respectively. The data is stored in the embedded SD card and retrieved after experiments. The sensor is placed on a foot as shown in Fig. 1. The trajectory is represented in the North-West-Up (NWU) earth frame, and the original point of the earth frame is set to the starting point of walking.

The entire PNS algorithm follows the framework of Foxlin [START_REF] Foxlin | Pedestrian tracking with shoe-mounted inertial sensors[END_REF], see Fig. 1. It is mainly composed of two modules, the ZUPT and INS (Inertial Navigation System):

-The ZUPT module takes data collected from the accelerometers and gyroscopes. Stance detection is firstly operated. Then the error-state EKF updates the velocity error δv and displacement error δr according to the result from stance detection. -The INS module has two main steps. First, it computes orientation by a quaternion method, which is based on a decomposed GDA. Second, it transfers the measured acceleration from the sensor frame (s) to the earth frame (e), for the use of integration to obtain velocity and displacement. At last, velocity and displacement are compensated by δv and δr from ZUPT module. Compared to many works of PNS [START_REF] Foxlin | Pedestrian tracking with shoe-mounted inertial sensors[END_REF] [4] [START_REF] Jiménez | Indoor pedestrian navigation using an INS/EKF framework for yaw drift reduction and a foot-mounted IMU[END_REF], the main difference comes from the separation of orientation estimation with the error-state EKF, which will provide more precision to the trajectory reconstruction, as explained in Section III.

III. INERTIAL NAVIGATION SYSTEM

Let us detail some aspects of the INS.

A. Magnetic Disturbance Detect

The magnetic field can be easily affected by other magnetic sources. Therefore, it is necessary to detect magnetic disturbance before estimating orientation. If e m = [ e m x , e m y , e m z ] , the measured magnetic field in the earth frame coordinates, satisfies the two criterions below, then the magnetic measurement cannot be used directly (see below).

| e m -e b | | e b | ≥ ε (1) arctan ( e m x / e m z ) -arctan ( e b x / e b z ) arctan ( e b x / e b z ) ≥ ε, (2)
in which the earth magnetic field e b = [ e b x , 0, e b z ] is supposed to be known. The threshold ε is set to 0.2 in this work.

B. Orientation Estimation

Let s e q = [ q0 , q1 , q2 , q3 ] , | s e q |= 1, denotes an estimated quaternion used to transfer coordinates from the earth frame to the sensor frame. A simple way to compute quaternion is integrating from its differential s e q, which is computed from angular rate. In this work, quaternion from the simple integration is denoted as s e q ω,n , ω and n respectively indicate the angular rate and time epoch. Quaternion is related to Euler angles, ie. yaw ψ, pitch θ and roll φ, [START_REF] Wang | Quaternion-based Kalman filter for AHRS using an adaptive-step gradient descent algorithm[END_REF] introduced the conversion from quaternion to Euler angles, ( ψ, θ, φ) = h ( s e q). Since gyroscopes have drift and noise, the quaternion integration error will inevitably increase as time elapses. So, [START_REF] Madgwick | Estimation of IMU and MARG orientation using a gradient descent algorithm[END_REF] propose a GDA-based quaternion estimation method as a guess to reduce the drift, if the time sampling interval is small enough. Their method is based on two objective functions related to acceleration and magnetic, f g and f b respectively, and on their corresponding Jacobian, J g and J b respectively.

f b and J b use the normalized earth magnetic field as a reference of north, the x-axis of the earth frame. If no disturbance is detected, the local earth magnetic field is referred to as the xaxis of the earth frame. Once a disturbance is confirmed by the method in section III-A at time epoch n, it regards e m n |s e qn-1 as the reference of the x-axis of the earth frame, e m n |s e qn-1 represents the magnetic measurement at time n that transferred into the earth frame by s e qn-1 . It is because that magnetic measurement will consist of extra magnetic resources once e b is interfered, a change in yaw angle is easier to be estimated if the new combined magnetic field is referred to as the x-axis of the earth frame. Therefore, transferring disturbed magnetic measurement back to the earth frame according to s e qn-1 will be an appropriate option.

The basic GDA-based quaternion estimation method merges the two objective functions into one. This may provide good estimation when the sensor is static or moving slowly, but may have bias when there is a large acceleration. It is because that the objective function for acceleration regards the normalized acceleration measurement as a reference of gravity g, different extents of motions may lead to different gradient values even if the objective function for magnetic is the same. It means that f g is influenced by acceleration and angular rate, but f b is only influenced by angular rate. Therefore, the gradient quaternion estimation needs to be computed separately:

s e q∇g,n = s e qn-1 -µ n J T g f g J T g f g , (3)
s e q∇ b ,n = s e qn-1 -µ n J T b f b J T b f b , (4)
where µ n = α s e qn ∆t, α > 1 is a ratio for the gradient, ∆t is the sampling interval. s e q∇g,n and s e q∇ b ,n are quaternion estimations based on GDA, it is reasonable because ∆t is small and quaternion does not change a lot within a tiny period. A data fusion method is then used to combine the estimations from eq. ( 3) and ( 4):

s e qg,n = ρ g,n • s e q∇g,n + (1 -ρ g,n ) s e q ω,n (5)
s e qb,n = ρ b,n • s e q∇ b,n + (1 -ρ b,n ) s e q ω,n (6)
ρ g,n = β g µ n /∆t + β g , ρ b,n = β b µ n /∆t + β b , (7)
s e qg,n and s e qb,n are the estimated quaternions w.r.t acceleration and magnetic, β g and β b represent the divergence rate of s e q ω,n based on the extent of motion (a larger value means a more intensive motion). Hence, they are dependent on the value of acceleration and angular rate. ρ g,n and ρ b,n are the weights yield to µ n and β g , β b .

Since pitch θ and roll φ depend on the estimation of gravity and yaw ψ depends on the estimation of north (the direction of magnetic field), it is possible to extract Euler angles from ( ψg , θg , φg ) = h( s e qg,n ) and ( ψb , θb , φb ) = h( s e qb,n ), then use the pitch and roll represented by s e qg,n and yaw represented by s e qb,n to get the final quaternion estimation:

s e qn = h -1 ψb , θg , φg . (8)

IV. ZERO-VELOCITY UPDATE

In ZUPT, stance detection is firstly performed for finding the stance phase, then the error of acceleration, velocity and displacement is updated by an error-state EKF, and at last the estimated error is used for compensating integration error during every stance.

A. Stance Detection

Generally, one walking step is broken down into four phases [START_REF] Mannini | A hidden Markov model-based technique for gait segmentation using a foot-mounted gyroscope[END_REF]: stance, push-up, swing and step-down. The four phases switch one after one and repeat periodically, thus we introduce a LR-HMM to represent the walking model, in which the state can only stay the same or switch to the next state.

Let's assume a hidden Markov chain model with observations Y = {Y 1 , . . . , Y N }, each Y n ∈ R 3 , and with unknown states X = {X 1 , . . . , X N }, each X n = k ∈ Ω = {1, . . . , 4}. Ω represents the stance, push-up, swing and step-down phases respectively. Assuming a discrete time independent Markov process, X can be parametrized by an initial probabilities vector π = p(x 1 ) and a transition matrix A(x 1 , x 2 ) = p(x 2 |x 1 ). In a cyclic LR-HMM, the transition matrix only allows states to keep the same or to switch to the next:

A(x 1 , x 2 ) =     1 -∆ 2 ∆ 2 0 0 0 1 -∆ 3 ∆ 3 0 0 0 1 -∆ 4 ∆ 4 ∆ 1 0 0 1 -∆ 1     , (9)
where ∆ k = p(x n = k|x n-1 = k -1) denotes the transition probability from state k -1 to state k.

The observation y n = ( s ω xn , s ω yn , s ω zn ) comes from the measurement of gyroscopes, which represents the angular rate along the three axes of the sensor frame. The distributions of observations conditional to states are assumed to be Gaussian

p (y n | x n = k) ∼ N (µ k , Σ k ), (10)
where µ k (3 × 1 vector) and Σ k (3 × 3 matrix) are the mean and co-variance of observations corresponding to state k. So that the LR-HMM model is parametrized by the following set of parameters

Θ = {π k , ∆ k , µ k , Σ k } k∈Ω .
Then the state transition probability p (x n |x n-1 ) and observation probability conditioned on state p (y n |x n ) can be easily derived by eq. ( 9) and [START_REF] Jiménez | Indoor pedestrian navigation using an INS/EKF framework for yaw drift reduction and a foot-mounted IMU[END_REF].

Learning of Θ parameters can be done by Baum-Welch algorithm [START_REF] Baum | A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains[END_REF], which is based on the EM principle (Expectation-Maximization) for finding the maximum likelihood iteratively, starting from an initial guess Θ (0) of parameters and stopping after a criterion or a maximum number of iterations is reached.

B. Error-State Extended Kalman Filter

Once the stance phase is successfully detected in Section IV-A, velocity is assumed to be zero and the velocity integration error can be obtained easily, unlike the displacement integration error which needs to be estimated appropriately. Based on the work of [START_REF] Jiménez | Indoor pedestrian navigation using an INS/EKF framework for yaw drift reduction and a foot-mounted IMU[END_REF], an EKF was used for estimating the integration error of velocity and displacement, so the elements in the state of EKF are the integration error and the EKF used here is also called as error-state EKF. In this work, we presents a simplified error-state EKF, since orientation estimation is performed independently in Section III-B, the error-state only takes into account the acceleration error (δ s a n ), velocity error (δ e v n ) and displacement error (δ e r n ):

δη n = δη n|n = δ e r n , δ e v n , δ s a n . ( 11
)
The error-state model is the one of [START_REF] Jiménez | Indoor pedestrian navigation using an INS/EKF framework for yaw drift reduction and a foot-mounted IMU[END_REF], δη n = Φ n δη n-1 + w n , where w n represents the process noise with covariance matrix Q n = E (w n w n ) and where the error-state transition model Φ n is a 9 × 9 matrix given by

Φ n =   I 3×3 ∆t • I 3×3 0 3×3 0 3×3 I 3×3 ∆t • e s C n 0 3×3 0 3×3 I 3×3   . ( 12
)
Φ n is time-variant and depends on the value of e s C n , which represents the rotation matrix required to convert vectors from the sensor frame to the earth frame at time n [START_REF] Wang | Quaternion-based Kalman filter for AHRS using an adaptive-step gradient descent algorithm[END_REF]. Now, the measurement equation writes z n = Hδη n|n + ν n , where z n is the measurement, H = [0 3×3 , I 3×3 , 0 3×3 ] is the measurement model, and ν n is the measurement noise, assumed as Gaussian with covariance

R n = E (ν k ν k ).
When the stance phase is detected, the error-state measurement in stance phase is z n = e v n -[0, 0, 0] (zero represents the real velocity), then the displacement integration error can be estimated because of the correlation between δ e v n and δ e r n described in the error-state model.

The error-state is predicted at every time, while is only updated during stance phases, since measurements are only available during stance phases. The prediction and update computation is illustrated in [START_REF] Simon | Optimal state estimation: Kalman, H infinity, and nonlinear approaches[END_REF]. Then the velocity and displacement can be compensated through e v n -δ e v n and e r n -δ e r n respectively, δ e r n and δ e v n in error-state should be reset to zero once the compensation is completed.

V. EXPERIMENTAL RESULTS

An experiment was conducted on a road nearby the campus of École Centrale de Lyon, Écully, France. The ground truth were obtained using OpenstreetMap, the total travel distance is 1075m with a walking speed of about 4.2km/s. The magnetic sensor was calibrated manually using the calibration toolbox proposed by Alex Faustino 2 .

Before starting to walk, a short time of standing without motion is necessary for initializing the quaternion corresponding to the earth frame, and the local magnetic field was obtained by transferring the mean value of magnetic measurements to the earth frame. The magnetic declination at Écully is 1.2 • . The LR-HMM method is tested and compared to another threshold based stance detection method detailed in [START_REF] Qiu | Inertial/magnetic sensors based pedestrian dead reckoning by means of multi-sensor fusion[END_REF]. Parameters of both algorithms are learned or tuned to get the best results.

The total steps number in experiment is 1458, so the steps number of one foot is 729. Compared with the threshold based stance detection method (see Fig. 2), LR-HMM obtains a more regular stance pattern, detecting rarely a false negative detection or missing rarely one step. Indeed, the missing and false negative detection number by LR-HMM are 0 and 6 (0.82% in total) respectively, while the numbers by threshold method are 17 and 14 (1.92% in total) respectively. Stance phases detected by LR-HMM have similar durations, which is not the case for the threshold-based method.

The proposed PNS algorithm involves LR-HMM to detect stance phase, and the decomposed GDA-based quaternion estimation (GDA d ) and simplified error-state EKF (EKF s ) to reconstruct trajectory. Comparative algorithms include the threshold stance detection, basic GDA-based quaternion estimation (GDA b ), a full error-state EKF (EKF f ) that contains compass measurement as described in [START_REF] Jiménez | Indoor pedestrian navigation using an INS/EKF framework for yaw drift reduction and a foot-mounted IMU[END_REF], and Zero Angular Rate Update (ZARU). Various combinations of these algorithms are tested and compared to the proposed algorithm, denoted here by LR-HMM+GDA d +EKF s +ZUPT.

All the parameters needed in the algorithms are tuned and learned for the best results, particularly β g and β b are set to 0.0215 and 0.026 respectively. In Fig. 4 and Table I, the proposed algorithm reconstructs a travel distance of 1298.9m, with an end-to-end error of 18.07m, the position error is 1.68%. It can be seen that independent orientation estimation provides a good accuracy in height, the proposed algorithm giving a final height error of 11.41m. GDA d -based quaternion estimation leads to a better result than the GDA b , which means that the decomposed GDA algorithm proposed in this paper can provide a much better dynamic performance than the basic GDA algorithm. The threshold-based stance detection has a good precision at the end point. However the highly enlarged travel distance (1376.9m in X-Y plane and 1464.4m in 3D space) is explained by some smaller stance phases duration than the truth. Thus more displacement compensation are made because of the unexpected extra integration error, which leads to an extra travel distance. The huge height error by EKF f is due to that it only takes the yaw angle into consideration, estimation for pitch or roll is not involved. At last, the DTW distance between the reconstructed trajectory and the groundtruth shows that the proposed algorithm reconstructed the bestfitted 3D trajectory.

VI. CONCLUSION

This paper presents a new strategy for PNS algorithm that uses a foot-mounted IMU sensor. Compared to the works proposed before, the orientation estimation is separated from the error-state EKF and a decomposed GDA-based quaternion estimation algorithm is employed for better dynamic performance. Due to the orientation is estimated independently, the error-state EKF can be simplified for reducing the complexity. Then a LR-HMM is used for stance detection for the use of ZUPT algorithm. Finally 3D trajectory is reconstructed by this new PNS algorithm structure.

The experimental results show that LR-HMM detects stance phase accurately and the proposed PNS algorithm reconstructs 3D trajectory accurately without extra method or prior knowledge. Our future work is to improve the dynamic performance of orientation estimation and to allow the algorithm to deal with different walking speeds during the same walking.

Figure 1 :

 1 Figure 1: Up: Shimmer3 sensor and its placement on a shoe. Down: Structure of PNS algorithm.

2 .Figure 2 :

 22 Figure 2: Step segmentation, jump 1 represents the stance phase. Stance detection by (a) a threshold method, where false negative detection and missing ones can be observed. (b) the proposed LR-HMM method.

Figure 3 :

 3 Figure 3: Reconstructed 3D trajectory by different algorithms.

Figure 4 :

 4 Figure 4: Reconstructed trajectory in X-Y plane by different algorithms.

Table I :

 I PNS Trajectory Error.

			Algorithms a			Ground
		1	2	3	4	5	Truth
	End-to-End Error (m)	18.45 16. 11 19.47 25.06 18.07	0
	Height Error (m)	9.81 241.39 145.10 13.84 11.14	0
	Total Distance in XY Plane (m)	1376.9 1199.6 1180.8 1303.1 1224.5 1075
	Total Distance in 3D (m)	1464.4 1452.3 1346.9 1387.3 1298.9 1097.6
	DTW Distance 1.149E4 16.36E4 9.514E4 1.348E4 1.279E4	-
	a 1 Threshold+GDA d +EKFs+ZUPT,
	2 LR-HMM+EKF f +ZUPT COMPASS(yaw) +ZARU,
	3 LR-HMM+EKF f +ZUPT COMPASS(yaw+inclination) +ZARU,
	4 LR-HMM+GDA b +EKFs+ZUPT,
	5 LR-HMM+GDA d +EKFs+ZUPT.

ACKNOWLEDGEMENT

We would like to thank the Chinese Scholarship Council (CSC) which partially supports this research.


François Castella
email: francois.castella@irisa.fr

Guillaume Dujardin
email: guillaume.dujardin@irisa.fr

Bruno Sericola
email: bruno.sericola@irisa.fr

MOMENTS ANALYSIS IN MARKOV REWARD MODELS

Keywords: Markov models, accumulated reward, performability, uniformization tsvp) Modèles markoviens, récompense cumulée, performabilité, uniformisation

We analyze the moments of the accumulated reward over the interval (0, t) in a continuous-time Markov chain. We develop a numerical procedure to efficiently compute the normalized moments using the uniformization technique. Our algorithm involves auxiliary quantities whose convergence is analyzed, and for which we provide a probabilistic interpretation.

I R I S A 1 Introduction

In the dependability analysis of repairable computing systems, there is an increasing interest in evaluating transient measures, and in particular, the accumulated reward over a given period. This measure is also known as a performability measure that takes into account both the performance of a system and its reliability. The state of the system is modeled by an irreducible continuous-time homogeneous Markov chain X = {X t , t ≥ 0}, over a finite state space S. To each state i ∈ S, is associated a reward rate f (i), which is assumed to be a nonnegative real number, without any loss of generality. The accumulated reward Y (t) over interval (0, t) averaged over t and normalized by the maximum reward rate f = max{f (j), j ∈ S} is defined, for every t > 0, by

Y (t) = 1 f t t 0 f (X s )ds.
The distribution of accumulated reward has been considered in several papers (see, for instance, [START_REF] Meyer | Closed-form solutions for performability[END_REF][START_REF] Iyer | Analysis of performability for stochastic models of fault-tolerant systems[END_REF][START_REF] Kulkarni | Numerical evaluation of performability and job completion time in repairable fault-tolerant systems[END_REF][START_REF] Ciardo | Performability analysis using semi-Markov reward processes[END_REF][START_REF] Smith | Performability analysis: measures, an algorithm, and a case study[END_REF][START_REF] Nabli | Performability analysis: a new algorithm[END_REF][START_REF] Souza E Silva | An algorithm to calculate transient distributions of cumulative rate and impluse based reward[END_REF][START_REF] Meyer | On evaluating the performability of degradable computing systems[END_REF][START_REF] Donatiello | On evaluating the cumulative performance distribution of fault-tolerant computer systems[END_REF][START_REF] Cloth | Five performability algorithms: A comparison[END_REF] and the references therein). This paper is concerned with the study and the evaluation of the moments of accumulated reward. As we shall see there is no need to evaluate the distribution if only some of the moments are required. Ther evaluation has been studied in [START_REF] Rácz | Numerical analysis of large markovian reward models[END_REF] where an algorithm has been proposed to compute them. In that paper the moments considered are not averaged over (0, t) so that they increase to infinity when t goes to infinity. The first moment of Y (t) has been analyzed in [START_REF] Sericola | Availability analysis of repairable computer systems and stationarity detection[END_REF], where a stationary regime detection has been developed to improve its computational time.

In this paper, we consider the moments averaged over (0, t). They are studied in Section 2 where we propose a first computational algorithm. Our approach, based on the uniformization technique, only involves numbers in [0, 1], which provides numerical stability. Besides, our method involves auxiliary quantities, called U(n, r) in the sequel, whose convergence is studied in Section 3 which allows us to develop a new procedure to compute the moments more efficiently. The U(n, r)'s are considered in Section 4 and they are characterized using a probabilistic interpretation. Section 5 is devoted to a numerical example.

Moments analysis

We denote respectively by α and Q the initial distribution and the infinitesimal generator of the irreducible continuous-time homogeneous Markov process X = {X t , t ≥ 0} over the finite state space S. For i ∈ S and r ≥ 0, we denote by m (r) i (t) the r-th moment of the accumulated reward Y (t) over (0, t), given that the initial state of X is equal to i, i.e.

m (r) i (t) = E(Y (t) r |X 0 = i).
We denote by m (r) (t) the column vector containing the m (r) i (t). Clearly, we have m (0) (t) = 1, where 1 is the column vector with all entries equal to 1, its dimension being clear in the context.

The r-th moment of the accumulated reward over interval (0, t) is thus given, for r ≥ 0, by

E(Y (t) r ) = αm (r) (t).
We denote by D the diagonal matrix defined by D = diag(d(i)), where d(i) = f (i)/f . The column vectors m (r) (t) can be obtained recursively from the following result.

PI n˚1869

Theorem 1. For every r ≥ 1 and t > 0, we have

m (r) (t) = r t t 0 e Q(t-u) Dm (r-1) (u)du. (1)
Proof. Let g r (s) be the random variable defined, for r ≥ 0 and 0 ≤ s ≤ t, by

g r (s) = t s f (X u )du r .
Differentiating with respect to s, we get

g r (s) = -rf (X s ) t s f (X u )du r-1 = -rf (X s )g r-1 (s),
which is equivalent to

g r (s) = r t s f (X u )g r-1 (u)du. ( 2
)
Note that from the homogeneity of the Markov chain X, we have, for r ≥ 0 and 0

≤ s ≤ t, E(g r (s)|X s = j) = E t s f (X u )du r |X s = j) = E t-s 0 f (X u )du r |X 0 = j) = f r t r m (r) j (t -s). ( 3
)
Taking the expectation in (2), we get, using the Fubini theorem

E(g r (s)|X 0 = i) = r t s E(f (X u )g r-1 (u)|X 0 = i)du = r t s j∈S E(f (X u )g r-1 (u)1 {Xu=j} |X 0 = i)du = r t s j∈S f (j)E(g r-1 (u)|X u = j, X 0 = i) Pr{X u = j|X 0 = i}du = r t s j∈S f (j)E(g r-1 (u)|X u = j) Pr{X u = j|X 0 = i}du = rf r-1 t r-1 t s j∈S f (j)m (r-1) j (t -u)e Qu (i, j)du,
where the fourth equality follows from the Markov property and the fifth is due to homogeneity property [START_REF] Cloth | Five performability algorithms: A comparison[END_REF]. Taking s = 0, we obtain

m (r) i (t) = r t t 0 j∈S d(j)m (r-1) j (t -u)e Qu (i, j)du,
Irisa which can be written, in matrix notation

m (r) (t) = r t t 0 e Qu Dm (r-1) (t -u)du,
or, by a variable change,

m (r) (t) = r t t 0 e Q(t-u) Dm (r-1) (u)du,
which completes the proof.

Differentiating relation [START_REF] Bowerman | Calculation of the Poisson cumulative distribution function[END_REF] with respect to t, we get

rm (r) (t) + t dm (r) (t) dt = tQm (r) (t) + rDm (r-1) (t). ( 4
)
We now make use of the uniformization technique [START_REF] Ross | Stochastic Processes[END_REF]. We introduce the uniformized Markov chain Z = {Z n , n ≥ 0} associated to the Markov chain X, which is characterized by its uniformization rate ν and by its transition probability matrix P . The uniformization rate ν verifies ν ≥ max{-Q(i, i), i ∈ S} and matrix P is related to the infinitesimal generator Q by P = I + Q/ν, where I denotes the identity matrix. The number of transitions during the interval (0, t), which we denote by N t , is a Poisson process with rate ν. Since Q = -ν(I -P ), relation (4) can be written as

rm (r) (t) + t dm (r) (t) dt = -νtm (r) (t) + νtP m (r) (t) + rDm (r-1) (t). (5)
In the following theorem, we determine the sequence of column vectors U(n, r) so that the solution to [START_REF] Donatiello | On evaluating the cumulative performance distribution of fault-tolerant computer systems[END_REF] has the form of the series (6).

Theorem 2. For every t ≥ 0, we have

m (r) (t) = ∞ n=0 e -νt (νt) n n! U(n, r), ( 6
)
where the column vectors U(n, r) = (U i (n, r), i ∈ S) are given by

U(0, r) = D r 1 for r ≥ 0, (7)
U(n, 0) = 1 for n ≥ 0, (8)
and, for n, r ≥ 1,

U(n, r) = n n + r P U(n -1, r) + r n + r DU(n, r -1). ( 9
)
Proof. Differentiating expression [START_REF] Iyer | Analysis of performability for stochastic models of fault-tolerant systems[END_REF] with respect to t, we get

t dm (r) (t) dt = -νtm (r) (t) + νt ∞ n=0 e -νt (νt) n n! U(n + 1, r) = -νtm (r) (t) + ∞ n=0 e -νt (νt) n n! nU(n, r).

PI n˚1869

François Castella, Guillaume Dujardin & Bruno Sericola

Replacing this expression together with expression (6) in equation ( 5), we obtain, for every t ≥ 0,

∞ n=0 e -νt (νt) n n! (n + r)U(n, r) = ∞ n=0 e -νt (νt) n+1 n! P U(n, r) + ∞ n=0 e -νt (νt) n n! rDU(n, r -1) = ∞ n=0 e -νt (νt) n n! (nP U(n -1, r) + rDU(n, r -1)).
Note that for n = 0, it is not necessary to define U(n -1, r) since, in that case, we have nP U(n -1, r) = 0. The equality holds for every t ≥ 0, so we get, for every n ≥ 0,

U(n, r) = n n + r P U(n -1, r) + r n + r DU(n, r -1).
The initial condition easily follows from the fact that the paths of the Markov chain X are supposed right continuous at t = 0. We thus have, for every r ≥ 0

lim t-→0 m (r) (t) = D r 1.
Taking t = 0 in relation ( 6), we get

U(0, r) = D r 1.
Since m (0) (t) = 1, we obtain U(n, 0) = 1, for every n ≥ 0, which completes the proof.

The uniformization method

The computation of the moments of the normalized accumulated reward E(Y (t) r ) is based on the following relation, which can be easily obtained from relation [START_REF] Iyer | Analysis of performability for stochastic models of fault-tolerant systems[END_REF],

E(Y (t) r ) = ∞ n=0 e -νt (νt) n n! u(n, r),
where u(n, r) = αU(n, r). Note that since matrix P is stochastic and since diagonal matrix D has its entries in the interval [0, 1], we easily obtain from relations ( 7), ( 8) and ( 9) that 0 ≤ U i (n, r) ≤ 1, and thus 0 ≤ u(n, r) ≤ 1.

Let ε > 0 be a given specified error tolerance and N be defined as

N = min n ∈ N n j=0 e -νt (νt) j j! ≥ 1 -ε . ( 10
)
Then we obtain

E(Y (t) r ) = N n=0 e -νt (νt) n n! u(n, r) + e(N),

Irisa

where the remainder of the series e(N) verifies

0 ≤ e(N) = ∞ n=N +1 e -νt (νt) n n! u(n, r) ≤ ∞ n=N +1 e -νt (νt) n n! = 1 - N n=0 e -νt (νt) n n! ≤ ε.
The computation of integer N can be made without any numerical problems even for large values of νt by using the method described in [START_REF] Bowerman | Calculation of the Poisson cumulative distribution function[END_REF].

The truncation level N is in fact a function of t, say N(t). For a fixed value of ε, N(t) is an increasing function of t. It follows that if we want to compute E(Y (t) r ) for J distinct values of t, denoted by t 1 < • • • < t J , we only need to compute u(n, r) for n = 1, . . . , N(t J ) since the values of u(n, r) are independent of the parameter t.

The pseudo code of the uniformization method can then be written as follows.

input : ε, t 1 < • • • < t J , R output : E(Y (t j ) r ) for j = 1, . . . , J and r = 1, . . . , R Compute N from Relation (10) with t = t J for r = 0 to R do U(0, r) = D r 1 u(0, r) = αU(0, r) endfor for n = 1 to N do U(n, 0) = 1 for r = 1 to R do U(n, r) = n n + r P U(n -1, r) + r n + r DU(n, r -1) u(n, r) = αU(n, r) endfor endfor for j = 1 to J do for r = 1 to R do E(Y (t j ) r ) = N n=0 e -νt j (νt j ) n n! u(n, r) endfor endfor Table 1: Algorithm for the computation of E(Y (t) r ).

Stationarity detection

We consider in this section the sequence of column vectors U(n, r) and we show that it converges when n tends to infinity. This allows us to stop the computation of the U(n, r) as soon as they are close enough to their limit. In the following theorem, we express the vectors U(n, r) recursively over index r. This recursive expression will be used in Theorem 4 to prove the convergence of the sequence U(n, r).

PI n˚1869

Theorem 3. For every n ≥ 0 and r ≥ 1,

U(n, r) = 1 n + r r n =0 + r -1 r -1 P n-DU( , r -1). ( 11
)
Proof. We prove this relation by recurrence over index n. For n = 0, this relation gives U(0, r) = DU(0, r -1). This leads to U(0, r) = D r U(0, 0) = D r 1, which is relation [START_REF] Kulkarni | Numerical evaluation of performability and job completion time in repairable fault-tolerant systems[END_REF].

Suppose that relation [START_REF] Rácz | Numerical analysis of large markovian reward models[END_REF] is true for integer n. We have to show that

U(n + 1, r) = 1 n + 1 + r r n+1 =0 + r -1 r -1 P n+1-DU( , r -1).
From relation ( 9) and using the recurrence hypothesis, we have

U(n + 1, r) = n + 1 n + 1 + r P U(n, r) + r n + 1 + r DU(n + 1, r -1) = n + 1 (n + 1 + r) n + r r n =0 + r -1 r -1 P n+1-DU( , r -1) + r n + 1 + r DU(n + 1, r -1) = 1 n + 1 + r r n =0 + r -1 r -1 P n+1-DU( , r -1) + r n + 1 + r DU(n + 1, r -1) = 1 n + 1 + r r n+1 =0 + r -1 r -1 P n+1-DU( , r -1),
which is the desired result.

We denote by π the stationary probability distribution of the Markov chain X. This row vector satisfies πA = 0 or equivalently π = πP . Theorem 4. For every r ≥ 0, we have

lim n-→∞ U(n, r) = (πD1) r 1. ( 12
)
Proof. We proceed by recurrence over integer r. The result is true for r = 0, since we have

U(n, 0) = 1.
Suppose the result is true for integer r -1, i.e. suppose that lim n-→∞ U(n, r -1) = (πD1) r-1 1.

Let us define

β(n, ) = + r -1 r -1 n + r r , H( ) = P D, U( ) = U( , r -1)
, and V (n) = U(n, r).

Irisa

We have β(n, ) ≥ 0 and

n =0 β(n, ) = 1,
and, since the sequence β(n, ) is increasing with ,

lim n-→∞ sup 0≤ ≤n β(n, ) = lim n-→∞ β(n, n) = lim n-→∞ r + 1 n + r + 1 = 0.
To ensure the convergence of the sequence of matrices P n , we require either that the output rates -Q(i, i), i ∈ S, of the Markov chain X are not all equal or that the uniformization rate ν is such that ν > max{-Q(i, i), i ∈ S}. This guarantees that the transition probability matrix P is aperiodic which means that it is ergodic and thus, we have

H = lim -→∞ H( ) = ΠD,
where Π is the matrix with all its lines equal to π. By using the recurrence hypothesis, we get

U = lim -→∞ U( ) = (πD1) r-1 1 = (ΠD) r-1 1.
From relation [START_REF] Rácz | Numerical analysis of large markovian reward models[END_REF] and by using Lemma 10, we get

lim n-→∞ U(n, r) = lim n-→∞ V (n) = ΠD(ΠD) r-1 1 = (ΠD) r 1 = (πD1) r 1,
which completes the proof.

Using this result, the algorithm of Table 1 can be modified as follows to take into account the convergence of the sequence U(n, r) and thus of the sequence u(n, r). In order to do that we consider that we make the following assumption, which is satisfied in practice.

if ∃ K < N s. t. max r≤R |u(K, r) -(πD1) r | ≤ ε then ∀n ≥ K, max r≤R |u(n, r) -(πD1) r | ≤ ε.
If such a K does not exist, there is no stationarity detection and we come back to the previous algorithm. Such a situation means that the value of N is not large and neither is the value of t J . If such a K exists, we have

E(Y (t) r ) = K n=0 e -νt (νt) n n! u(n, r) + (πD1) r 1 - K n=0 e -νt (νt) n n! + e 1 (K),
where the remainder of the series e 1 (K) verifies

|e 1 (K)| = ∞ n=K+1 e -νt (νt) n n! (u(n, r) -(πD1) r ) ≤ ∞ n=K+1 e -νt (νt) n n! |u(n, r) -(πD1) r | ≤ ε.
We obtain the following algorithm, in which we suppose that πD1 has been computed beforehand.

PI n˚1869 input : ε, t 1 < • • • < t J , R
output : E(Y (t j ) r ) for j = 1, . . . , J and r = 1, . . . , R Compute N from Relation [START_REF] Nabli | Performability analysis: a new algorithm[END_REF] with t = t J for r = 0 to R do

U(0, r) = D r 1 u(0, r) = αU(0, r) endfor K = N for n = 1 to N do U(n, 0) = 1 for r = 1 to R do U(n, r) = n n + r P U(n -1, r) + r n + r DU(n, r -1) u(n, r) = αU(n, r) endfor if |u(n, r) -(πD1) r | ≤ ε for every r = 1, . . . , R then K = n break endif endfor for j = 1 to J do if K = N then for r = 1 to R do E(Y (t j ) r ) = N n=0 e -νt j (νt j ) n n! u(n, r) endfor endif if K < N then for r = 1 to R do E(Y (t j ) r ) = K n=0 e -νt j (νt j ) n n! u(n, r) + (πD1) r 1 - K n=0
e -νt j (νt j ) n n! endfor endif endfor Table 2: Algorithm for the computation of E(Y (t) r ) using stationarity detection.

The following section is devoted to probabilistic interpretations of the sequence U(n, r).

Probabilistic interpretation

In the following theorem, we obtain an explicit expression of the column vectors U(n, r) using the matrices P and D. Irisa Theorem 5. For every n ≥ 0 and r ≥ 1, we have

U(n, r) = 1 n + r r n 1 =0 P n-1 D 1 2 =0 P 1 -2 D • • • r-2 r-1 =0 P r-2 -r-1 D r-1 r =0 P r D1. ( 13
)
Proof. For r = 1, the expression (13) reads

U(n, 1) = 1 n + 1 n 1 =0 P 1 D1,
which is obtained by [START_REF] Rácz | Numerical analysis of large markovian reward models[END_REF] and [START_REF] Meyer | On evaluating the performability of degradable computing systems[END_REF].

Suppose that relation ( 13) is true for integer r -1, i.e. suppose that, for every 1 ≥ 0,

U( 1 , r -1) = 1 1 + r -1 1 -1 1 2 =0 P 1 -2 D 2 3 =0 P 2 -3 D • • • r-2 r-1 =0 P r-2 -r-1 D r-1 r =0 P r D1.
Using relation [START_REF] Rácz | Numerical analysis of large markovian reward models[END_REF], we obtain directly [START_REF] Sericola | Availability analysis of repairable computer systems and stationarity detection[END_REF].

The first interpretation of the U(n, r) we obtain is given by the following theorem.

Theorem 6. For every n ≥ 0, r ≥ 1 and i ∈ S, we have

U i (n, r) = 1 n + r r E   n 1 =0 d(Z 1 ) n 2 = 1 d(Z 2 ) • • • n r-1 = r-2 d(Z r-1 ) n r = r-1 d(Z r ) Z 0 = i   . ( 14
)
Proof. The relation is true for r = 1 since the right hand side becomes

1 n + 1 E n 1 =0 d(Z 1 ) Z 0 = i = 1 n + 1 n 1 =0 E (d(Z 1 )|Z 0 = i) = 1 n + 1 n 1 =0 j∈S P 1 i,j d j = 1 n + 1 n 1 =0 (P 1 D1) i = U i (n, 1).
Suppose that the relation is true for integer r -1; i.e. that

U i (n, r -1) = E   n 1 =0 d(Z 1 ) n 2 = 1 d(Z 2 ) • • • n r-2 = r-3 d(Z r-2 ) n r-1 = r-2 d(Z r-1 ) Z 0 = i   n + r -1 r -1 . PI n˚1869
We then have

E   n 1 =0 d(Z 1 ) n 2 = 1 d(Z 2 ) • • • n r-1 = r-2 d(Z r-1 ) n r = r-1 d(Z r ) Z 0 = i   = j∈S E   n 1 =0 d(Z 1 ) n 2 = 1 d(Z 2 ) • • • n r-1 = r-2 d(Z r-1 ) n r = r-1 d(Z r ) Z 1 = j, Z 0 = i   P 1 i,j = n 1 =0 j∈S P 1 i,j d(j)E   n 2 = 1 d(Z 2 ) • • • n r-1 = r-2 d(Z r-1 ) n r = r-1 d(Z r ) Z 1 = j, Z 0 = i   = n 1 =0 j∈S P 1 i,j d(j)E   n 2 = 1 d(Z 2 ) • • • n r-1 = r-2 d(Z r-1 ) n r = r-1 d(Z r ) Z 1 = j   = n 1 =0 j∈S P 1 i,j d(j)E   n-1 2 =0 d(Z 2 ) • • • n-1 r-1 = r-2 d(Z r-1 ) n-1 r = r-1 d(Z r ) Z 0 = j   = n 1 =0 j∈S P 1 i,j d(j) n -1 + r -1 r -1 U j (n -1 , r -1) = n 1 =0 1 + r -1 r -1 j∈S P n-1 i,j d(j)U j ( 1 , r -1) = n 1 =0 1 + r -1 r -1 (P n-1 DU( 1 , r -1)) i = n + r r U i (n, r),
which completes the proof. Note that the third equality follows from the Markov property, the fourth equality follows from the homogeneity of the Markov chain {Z}, the fifth one follows from the recurrence hypothesis and the last one is relation [START_REF] Rácz | Numerical analysis of large markovian reward models[END_REF].

Clearly, for r = 1, U i (n, 1) is the expectation of of the number of visits to each state i during the n first transitions, weighted by the reward d(i), averaged by n + 1. For r ≥ 2, such a simple interpretation of U i (n, r) is not at hand. The remainder part of this section is devoted to establishing another way to obtain a probabilistic interpretation of U i (n, r). The two following lemmas will be used to prove Theorem 9. (n + + 1) n + .

Proof. Let I(n, ) denote the integral in the right hand side. Using an integration by parts, we obtain

I(n, ) = n + 1 I(n -1, + 1).

Irisa

This leads to

I(n, ) = n + I(0, n + ) = n + t n+ +1 n + + 1 ,
which completes the proof.

We define

G i (t, n, r) = E(Y (t) r |N t = n, X 0 = i).
Lemma 8. For every t ≥ 0, n ≥ 1, r ≥ 0 and i ∈ S, we have

G i (t, n, r) = n t n+r r =0 r d(i) j∈S p i,j t 0 (t -u) u n+r--1 G j (u, n -1, r -)du.
Proof. We introduce the quantity

V i (t, n, r) = E(Y (t) r 1 {Nt=n} |X 0 = i).
Let T 1 be the sojourn time in the initial state. Using a renewal argument, we have

V i (t, n, r) = E(Y (t) r 1 {Nt=n} |X 0 = i) = t 0 j∈S p i,j E(Y (t) r 1 {Nt=n} |X u = j, T 1 = u, X 0 = i)νe -νu du = t 0 j∈S p i,j E f i u + t u f (X s )ds r 1 {Nt-Nu=n-1} |X u = j, T 1 = u, X 0 = i νe -νu du = t 0 j∈S p i,j E f i u + t u f (X s )ds r 1 {Nt-Nu=n-1} X u = j νe -νu du = t 0 j∈S p i,j E [f i u + Y (t -u)] r 1 {Nt-Nu=n-1} |X u = j νe -νu du = t 0 j∈S p i,j E [f i u + Y (t -u)] r 1 {N t-u =n-1} |X 0 = j νe -νu du = r =0 r f (i) j∈S p i,j t 0 u E Y (t -u) r-1 {N t-u =n-1} |X 0 = j νe -νu du = r =0 r f (i) j∈S p i,j t 0 u V j (t -u, n -1, r -)νe -νu du.
Unconditioning on the number of transitions during (0, t), we get

G i (t, n, r) = n!V i (t, n, r) f r t r e -νt (νt) n ,
and thus

G i (t, n, r) = n t n+r r =0 r d(i) j∈S p i,j t 0 u (t -u) n+r--1 G j (t -u, n -1, r -)du = n t n+r r =0 r d(i) j∈S p i,j t 0 (t -u) u n+r--1 G j (u, n -1, r -)du,
which completes the proof.

PI n˚1869

The following theorem gives a clearer interpretation of U i (n, r) in terms of the normalized accumulated reward Y (t) and of the Poisson process N t . Theorem 9. For every n ≥ 0, r ≥ 0, i ∈ S and t ≥ 0, we have

U i (n, r) = E(Y (t) r |N t = n, X 0 = i). ( 15
)
Proof. Conditioning on N t , we have

E(Y (t) r |X 0 = i) = ∞ n=0 e -νt (νt) n n! E(Y (t) r |N t = n, X 0 = i),
and from ( 6), we have

E(Y (t) r |X 0 = i) = ∞ n=0 e -νt (νt) n n! U i (n, r).
So, to prove the theorem, we just have to show that E(Y (t) r |N t = n, X 0 = i) is independent of t. We proceed by recurrence over index n.

For n = 0, we have

E(Y (t) r |N t = 0, X 0 = i) = d(i) r since Y (t) = f (i)t when N t = 0 and X 0 = i. Suppose that U i (n -1, r) = E(Y (t) r |N t = n -1, X 0 = i).
Using this hypothesis and Lemma 8, we obtain

G i (t, n, r) = n t n+r r =0 r d(i) j∈S p i,j U j (n -1, r -) t 0 (t -u) u n+r--1 du.
Using now Lemma 7, we obtain

G i (t, n, r) = n t n+r r =0 r d(i) j∈S p i,j U j (n -1, r -) t 0 (t -u) u n+r--1 du = n n + r r =0 r d(i) n -1 + r j∈S p i,j U j (n -1, r -), which is independent of t.
It is worthwhile noting that the moments of the normalized accumulated reward Y (t), given the number of transitions over (0, t) in the uniformized process, are independent of t.

In the next section, we present some numerical experiments to show the importance of the stationarity detection in the reduction of the computational time of the moments E(Y (t) r ).

Numerical example

We consider a fault-tolerant multiprocessor system which consists of n identical processors and b buffer stages. Processors fail independently at rate λ and are repaired singly with rate µ. Buffers stages fail independently at rate γ and are repaired with rate τ . Processor failures cause a graceful degradation of the system and the number of operational processors is decreased by Irisa Moments analysis in Markov reward models 15 one. The system is in a failed state when all the processors have failed or any of the buffer stages has failed. No additional processor failures are assumed to occur when the system is in a failed state. The model is represented by a Markov process with state transition diagram shown in Figure 1. The state space of the system is S = {(i, j); 0 ≤ i ≤ n, j = 0, 1}. The component i of a state (i, j) means that there are i operational processors and the component j is zero if any of the buffer stages has failed, otherwise it is one. It follows that the set U of operational states is U = {(i, 1); 1 ≤ i ≤ n}. The reward structure we choose here is given by f (i) = 1 if i ∈ U and f (i) = 0 otherwise. We suppose that the initial state of the system is state (n, 1). The number of processors is fixed to 16, each with a failure rate λ = 0.00006 per hour and a repair rate µ = 0.1666 per hour. The number of buffer stages is fixed to 1024, each with a failure rate γ = 0.00131 per hour and a repair rate τ = 0.1666 per hour. The error tolerance is ε = 0.00001. In Figure 2, we plot the first five moments Y (t), as a function of t. Let us now consider higher values of t. Solving the linear system πA = 0 with π1 = 1, we obtain πD1 = 0.110475. The complexities of the two algorithms in Table 1 and Table 2 are mainly due to the products of matrix P by vector U(n -1, r). The number of such products is equal to NR for the algorithm of Table 1 and it is equal to KR, with K ≤ N, for the algorithm of Table 2. We show in Figure 3, the values of K and N obtained when R = 5, for different values of t = t J . As expected the value of N = N(t J ) increases with t J and the value of K constant after the first instant where K < N. This instant which we denote by T ε , does dependent on ε, and may be called the time to stationarity. It is defined as

1,1 0,1 λ τ τ τ τ γ γ γ b b b γ b µ µ µ µ µ µ µ µ λ λ 2 n,1 n-1,1 n,0 0,0 (n-1) n-1,0 1,0 λ n
T ε = inf{t ≥ 0|K < N(t)}.
In our example, it is between 50000 and 60000. A more detailed analysis shows that T ε = 55482. We start with the second term. Since L ≤ ≤ n -L, we have ≥ L and n -≥ L, we can write

H(n -)U( ) -HU = [H(n -) -H]U( ) + H[U( ) -U] ≤ H(n -) -H U( ) + H U( ) -U ≤ M ( H(n -) -H + U( ) -U ) ≤ ε 2 .
We thus have Concerning the first and the third terms, we use the fact that, for every 0 ≤ ≤ n, we have

H(n -)U( ) -HU = H(n -)U( ) + HU ≤ H(n -) U( ) + H U ≤ 2M 2 .
Let us define the sequence β(n) = sup 0≤ ≤n β(n, ). Since, by hypothesis, β(n) converges to 0 when n goes towards infinity and since L is fixed, we can determine N such that for any n ≥ N, we have β(n) ≤ ε/(8M 2 L).

We then have Putting together these two results, we obtain for every n ≥ max(2L, N), V (n) -HU ≤ ε, which completes the proof.

PI n˚1869

  A T O I R E S P U B L I C A T I O N I N T E R N E N o

Lemma 7 .

 7 For every t ≥ 0, we havet 0 (tu) u n du = t n+ +1

Figure 1 :

 1 Figure 1: State-transition diagram for a n-processor system.

t

  50000 60000 70000 80000 90000 100000 N 76621 91823 107015 122200 137379 152154 K 76621 84955 84955 84955 84955 84955

Figure 3 :

 3 Figure 3: Stationarity detection for different values of t when R = 5.

  , )[H(n -)U( ) -HU] ≤ n-L =L β(n, ) H(n -)U( ) -

L- 1 =0β≤ L- 1 =0β

 11 (n, )[H(n -)U( ) -HU] + n =n-L+1 β(n, )[H(n -)U( ) -HU] (n, ) H(n -)U( ) -HU] + n =n-L+1 β(n, ) H(n -)U( ) -HU]

Appendix

Let β(n, ), for n ≥ 0 and 0 ≤ ≤ n, be real numbers, H( ) be a sequence of matrices and U( ) be a sequence of column vectors with the same finite dimension We define the sequence of column vectors V (n) by

We have the following result.

Lemma 10. Assume that

Proof. The convergence of the sequences H( ) and U( ) implies that both sequences are uniformly bounded, i.e. there exists an integer M such that, for every ≥ 0,

Moreover, their convergence implies that for any ε > 0, there exists L such that for any , m ≥ L, we have

Irisa


Emilie Soret
email: emilie.soret@ircica.univ-lille1.fr

Laurent Clavier
email: laurent.clavier@telecom-lille.fr

Gareth W Peters

Ido Nevat

François Septier

Émilie Soret

SIMO communication with impulsive and dependent interference -the Copula receiver

come   L'archive ouverte pluridisciplinaire

Introduction

Impulsive interference are encountered in many situations, e.g. power line communications, ultra-wide band technology, or in dense networks. This significantly degrades the performance of the classical receivers [START_REF] Benmâad | Asymptotic performance of LDPC codes in impulsive non-gaussian channel[END_REF]. In this paper, we consider a simple detection problem in a block fading scenario. Each data symbol is transmitted over wireless channels and K = 2 versions of each symbol are received. We only consider the case K = 2 for clear analytical expressions and simple illustrations but this can be extended to higher dimensions. This transmission structure can be motivated by many different practical wireless communication systems like a rake receiver [START_REF] Niranjayan | The BER optimal linear rake receiver for signal detection in symmetric α-stable noise[END_REF], a Single-Input-Multiple-Output (SIMO) system [START_REF] Filippou | Optimal combining of instantaneous and statistical CSI in the SIMO interference channel[END_REF], a cooperative communication scheme involving multiple relays [START_REF] Chen | Alpha-stable multiple access interference modelling for amplify-and-forward multihop ad hoc networks[END_REF] or in impulse radio Ultra Wide Band systems with symbol repetitions [START_REF] Ghannudi | Alpha-stable interference modeling and cauchy receifor an ir-uwb ad hoc network[END_REF].

For a single transmitted symbol s n at time n, the received signal Y ∈ R K is Y = s n h n + I k + N n , where h n ∈ R K is the block fading channel coefficients, I n ∈ R K is the impulsive interference and N n i.i.d.

∼ N (0, σ 2 ) is the thermal noise. In this paper we make the assumption that the channel state is perfectly known and that interference is dominating. Besides we assume independence between different time instant n so that we will drop this index for clarity. The studied case can then be summarized by Y S + I, where S is a vector containing the repeated sample s and I the interference vector.

Many papers have considered the case where I is composed of independent and identically distributed samples. Depending on the impulsive interference distribution assumption, it is more or less complicated to derive the optimal receiver and sometimes suboptimal approaches are considered [START_REF] Beaulieu | Designing time-hopping ulwide bandwidth receivers for multiuser interference environments[END_REF].

In this paper we consider an α-stable interference distribution but we do not consider any longer that the components of I are independent. We take the example of a SIMO link : if a strong interference is received on one antenna, the probability of receiving a strong interference sample on another antenna is not negligible. This upper tail dependence can not be captured by traditional correlation function that, anyway, can not be used for α-stable random vectors. We propose to use the copula framework to model the dependence structure. It allows to separately model the marginal distributions and the dependence structure.

Copulae are a very useful way to model structures of dependence between random variables [7]. The fundamental result justifying this usefulness is the Sklar's Theorem : it ensures that under the condition that the cumulative distributions of the random variables are continuous, there exists a unique copula C such that ∀(x 1 , . . . , x d ), we have

H(x 1 , . . . , x d ) = C (F 1 (x 1 ), . . . , F d (x d ))) . ( 1
)
where H is the joint distribution of (X 1 , . . . , X d ). Hence, a copula is a function C : [0, 1] d → [0, 1] which couples the marginals F i between themselves. The name copula comes from this last remark. In Fig. 1 we represent the interference samples when I has independent components. The representation is done directly on the sample or after a transformation through the repartition function of the marginals (F i (.)) to have the representation of the copula. Received samples in the Copula space FIGURE 1 -In the left plot samples are independent and they form a cross. This can be explained saying that large values are rare and the occurrence of two large values on the same vector is very unlikely. In the right plot (X and Y axis are F i (y 1 ) and F i (y 2 )), the points are uniformly distributed which signifies the independent structure.

Archimedean copulae

In the following, we consider a particular class of bivariate Archimedean copulae. The interest of this class is, first of all, the easiness with which they can be constructed. The multivariate Archimedean copulae have the following form : for all

(u 1 , • • • , u d ) ∈ [0, 1] d , C(u 1 , • • • , u d ) = φ -1 (φ(u 1 ) + • • • + φ(u d )) . (2)
The function φ is called the generator of the copula and is a continuous and convex function such that φ(1) = 0. It appears that all Archimedean copula is symmetric in its variables. We will focus on two families of Archimedean copulae, both indexed by a single parameter. The Clayton and the Gumbel families of copulae model asymmetric dependence in tails. Definition 2.1. For all θ > 0,The Clayton copula of parameter

θ is defined on [0, 1] d by C(u 1 , • • • , u d ) = u -1/θ 1 + • • • + u -1θ d -(d -1) -θ .
In particular, it is obtained when φ -1 is the Laplace transform of a Gamma distribution.

In Fig. 2 we have a similar representation as in Fig. 1 but introducing the dependence structure of the Clayton copula. The

C(u 1 , • • • , u d ) = exp -( d i=1 (-log u i )) 1/θ .
In particular, it is obtained when φ -1 is the Laplace transform of a α-stable distribution.

In Fig. 3 we represent the dependence structure of the Gumbel copula on the received samples and after the transform through the marginals.

LLR for dependent variables

In the two-dimensional case and with a binary input, our system model in section 1can be written

y 1 = s + i 1 y 2 = s + i 2 , (3)
where s ∈ {-1, 1}. Two repetitions y 1 and y 2 of this transmitted bit are obtained and I = (i 1 , i 2 ) is a bivariate interference vector. The two coordinates i 1 and i 2 are not independent. The Log Likelihoood Ratio (LLR) for each Y ∈ R 2 is given by

Λ(y 1 , y 2 ) = log P(y 1 = s + i 1 , y 2 = s + i 2 | s = 1) P(y 1 = s + i 1 , y 2 = s + i 2 | s = -1) . (4)
Let f be the joint density of the couple (i 1 , i 2 ), (4) becomes

Λ(y 1 , y 2 ) = log f (y 1 -1, y 2 -1) f (y 1 + 1, y 2 + 1) . (5)

Independent interferences

In the left plot in Fig. 4, we illustrate the two decision regions in the case when interference is independent on the two dimensions and Cauchy distributed with location x 0 = 0 and scale δ = 1. The X and Y axis are the values of the components of the received vector Y. We consider two possible transmitted symbols {-1, 1} meaning that the transmitted vector corresponds to the points (1, 1) and (-1, -1). The white region corresponds to the decision 1, meaning that Λ ≥ 0 and the black one to -1, i.e., Λ < 0 A Gaussian noise would correspond to a linear boundary, corresponding to an Euclidean distance. Impulsiveness significantly modifies those boundaries and necessitate non linear operation to implement an optimal receiver.

Dependent interferences

If we now consider that i 1 and i 2 are dependent and that we can express this dependence through an Archimedean copula, the LLR will become

Λ(x, y) = log f i (x -1)f i (y -1)c(F i (x -1), F i (y -1)) f i (x + 1)f i (y + 1)c(F i (x + 1), F i (y + 1)) = Λ ⊥ (x, y) + Λ c (x, y), (6)
where c is the density of the copula and is defined by

c(u, v) = ∂ 2 C ∂u∂v (u, v); (7)
f i and F i are respectively the probability density function and the cumulative distribution of the interference. Λ ⊥ represents the independent part of the LLR. The second term

Λ c (x, y) = log c(F i (x -1), F i (y -1)) c(F i (x + 1), F i (y + 1)) (8)
is the part of the LLR depending on the copula and represents the dependence structure. It can however be tricky to derive because it also depends on the marginals.

In the case of the Clayton copula, the consequence on the decision region is shown in Fig. 4. We clearly see that the lower tail dependence significantly modifies the decision regions.

In the case of the Gumbel copula, the consequence on the decision region is shown in Fig. 5. We again clearly see that the lower tail dependence significantly modifies the decision regions.

Application to SIMO transmissions 4.1 Receiver design

The optimal receiver in terms of minimizing the Bit Error Rate (BER) is the Maximum Likelihood (ML) detector ŝ = FIGURE 5 -Decision region for independent Cauchy, Cauchy marginals and Gumbel copula and the difference between both (in white the areas where the dependence structure modifies the optimal decision).

arg max s∈Ω P(Y|s), where Ω is the set of possible transmitted bits. In the binary case, Ω = {-1, 1} and the problem is reduced to obtaining the sign of the LLR defined in [START_REF] Beaulieu | Designing time-hopping ulwide bandwidth receivers for multiuser interference environments[END_REF].

s = sign (Λ(x, y)) = sign (log Λ ⊥ (x, y) + Λ c (x, y)) , (9)
Fig. 6 compares the performance of the linear Gaussian receiver, a Cauchy receiver assuming independent Cauchy interference and a copula receiver that knows both the marginal and the dependence structures. In that case the dependence is captured by a Clayton copula. Obviously when the parameter gets close to zero, the dependence is low, the Cauchy receiver outperforms the Gaussian receiver and there is no need to introduce the dependence structure. However, when the dependence increases (θ gets larger), the performance of the Cauchy receiver quickly degrades when the copula receiver is able to maintain a better performance level. Fig. 7 shows similar results with the Gumbel copula. We finally apply our Copula receiver to the SIMO case. Interferers are uniformly distributed in a square around the recei- ver. The square is 10 by 10 and the mean number of interferers is 50. We consider normalized distances so the unity is not significant. The channels are Rice channels with a main path strength randomly chosen. The channel attenuation coefficient is 3. For the copula receiver, we chose a mixture of the Gumbel and Clayton copulas to ensure a symmetry in the upper and lower tail dependence. To have the symmetry, the parameters for the two copula are linked so that only one parameter has to be chosen to define the dependence structure. To observe the impact of including the dependence in the receiver, we vary this parameter. Fig. 8 shows the benefit of including the dependence structure to design the receiver. The gain in performance is limited but this is easily explained by the small dimension considered (only 2). Besides, we chose a Cauchy distribution for the marginals which is not the optimal choice. It has however proved to be close to the optimal in several situations [START_REF] Ghannudi | Alpha-stable interference modeling and cauchy receifor an ir-uwb ad hoc network[END_REF]. It is clear that taken the dependence structure into account allow a further gain compared to the independent receiver, which already gives better performance than the Gaussian receiver.

Conclusion

We proposed in this paper a way to model dependency in impulsive interference. Usual tools (based on correlation) do not allow to well capture the dependence structure of such an impulsive interference, especially when the α-stable model is used.

In the case of Cauchy marginals and copulae from the Archimedean family and with a binary input, we are able to derive analytical expressions of the decision rule based on the likelihood ratio. The results on the decision regions show that dependent interference has a significant impact on the optimal decision that we should make. Consequently, we compared receivers that takes this dependency into account to receivers that do not. We show that the latter can rapidly degrade if a dependence structure is present when the former manage to maintain good performances. We illustrate the possible benefit on a SIMO example. Many questions are still open, what is the best choice for the copula ? How can we estimate the interference parameters ? How can we implement such a receiver, especially if the dimension increases ?

The densification of networks and their heterogeneity make interference an important issue in wireless communication. The dependence structure is certainly a crucial point for an efficient implementation of such networks. Will copula play a role ?

Acknowledgement

This work is part of the research project PERSEPTEUR (financed by the French Agence Nationale de la Recherche ANR) and supported by IRCICA, USR CNRS 3380.

2 )

 2 Received samples in the copula plane

FIGURE 2 -

 2 FIGURE 2 -Interference samples for Cauchy marginals and Clayton copula. cross in the left plot tends to disappear and points, especially in the bottom left quadrant, are differently positioned. This results from the non zero asymmetric tail dependence introduced by the Clayton copula. Definition 2.2. For all θ ≥ 1,The Gumbel copula of parameter θ is defined on [0, 1] d by

FIGURE 3 -

 3 FIGURE 3 -Received samples for Cauchy marginals and Gumbel copula.

4 FIGURE 4 -

 44 FIGURE 4 -Decision region for independent Cauchy, Cauchy marginals and Clayton copula and the difference between both (in white the areas where the dependence structure modifies the optimal decision).

FIGURE 6 -

 6 FIGURE 6 -BER for Cauchy noise and Clayton copula as a function of the Clayton parameter.

FIGURE 7 -

 7 FIGURE 7 -BER for Cauchy noise and Gumbel copula as a function of the Gumbel parameter.

FIGURE 8 -

 8 FIGURE 8 -BER for SIMO model for mixture of Clayton and Gumbel copulae.


• • ` ` • Amaël Dupaix • Aurélie Guislain-Vitek • Aurélien Levret • Bérangère Foscarin • Caroline Tatard • Dominique Rieffel • Eve Afonso • Geoffroy Couval • Herminie Piernavieja • Javier Fernandez-de-Simon • Jean-Baptiste Pons • Jean-Francois Cosson • José Gonzalo Zapata Carbonell • Julie Montas • Lexa Villette • Marilyne Fontanier • Maxime Galan • Nathalie Charbonnel • Nicolas Bassin • Patrick Giraudoux • Thibault Moulin • Virgile Baudrot • Vincent Driget


Pierre-Loup Schaefer

Grégory Chagnon

Alexandre Moreau-Gaudry

Geometry-Based Model for U-Shaped Strain Gauges on Medical Needles

The knowledge of needle location during insertions is essential for the success of interventional radiology procedures. As the needle is susceptible to undergo deformations during its insertion into tissues, several methods have been proposed to monitor the needle deformed shape. Thus, instrumented needles with U-shaped strain gauges are currently being developed to reconstruct the shape of the needle from gauge acquisitions. These acquisitions are used in combination with gauge model to obtain estimate of the strain of the needle. The current modeling is limited as it does not consider the geometry of the gauge. This paper introduces a more complete model for U-shaped strain gauge which, on the difference of the current model, takes into account the width of the gauge. Thus, the impact of width modeling on the strain estimate can be measured and used to improve strain estimation accuracy. Results with real characteristics of instrumented needle devices show that the differences of strain estimate are around few percents. Finally, by taking into account the width and the length of the gauge our model includes the effets of the gauge size on the strain estimation and makes the miniaturization of the gauge less necessary.

I. INTRODUCTION

Needles insertions are common procedures carried out during interventional radiology surgeries. During these insertions, the needle can be subject to deformation due to its interaction with the tissues in the patient body. The deformations can result in a loss of accuracy of the needle location which might affect the success of the surgical procedure. Different technical solutions have been proposed in order to know the deformed shape of the needle. The physician can then use medical imaging, such as echography, x-ray or CT imaging to acquire images surrounding the insertion area and visualize the deformed needle shape. A recent technique, proposed to calculate the deformed shape of the needle, involves the use of strain measures acquired from strain sensors on the surface of the needle. Different kinds of strain sensors can be used such as Fiber Bragg gratings [START_REF] Abayazid | 3d flexible needle steering in soft-tissue phantoms using fiber bragg grating sensors[END_REF], [START_REF] Park | Realtime estimation of 3-d needle shape and deflection for mri-guided interventions[END_REF], [START_REF] Roesthuis | Three-dimensional needle shape reconstruction using an array of fiber bragg grating sensors[END_REF], [START_REF] Seifabadi | Real-time tracking of a bevel-tip needle with varying insertion depth: Toward teleoperated mri-guided needle steering[END_REF] or strain gauges as in the work of Bonvilain et al. [START_REF] Bonvilain | Characterization of strain microgauges for the monitoring of the deformations of a medical needle during its insertion in human tissues[END_REF] and Hammond et al. [START_REF] Hammond | Estimating surgical needle deflection with printed strain gauges[END_REF].

In the case of strain gauge technology, the U-shaped gauges are printed [START_REF] Hammond | Estimating surgical needle deflection with printed strain gauges[END_REF] or micro-etched [START_REF] Bonvilain | Characterization of strain microgauges for the monitoring of the deformations of a medical needle during its insertion in human tissues[END_REF] parallel to the needle in order to measure the axial strain on the surface, as illustrated in Fig. 1. As the needle undergoes deformations, the variation of the length of the gauge causes a change of its resistance. Thus, the strain ε gauge indicated by a gauge of initial resistance R gauge and resistance change ∆R gauge can be obtained using the relation:

ε gauge = ∆R gauge R gauge × GF (1)
1 Univ. Grenoble Alpes / CNRS / TIMC-IMAG UMR 5525, France.

Email: FirstName.Lastname@univ-grenoble-alpes.fr The term GF is the gauge factor and depends on the material characteristics of the gauge. According to the size of the gauge, hypothesis can be made on the dimensions to be considered. Thus, the length of the gauge can be taken into account [START_REF] Hammond | Estimating surgical needle deflection with printed strain gauges[END_REF] or not [START_REF] Bonvilain | Characterization of strain microgauges for the monitoring of the deformations of a medical needle during its insertion in human tissues[END_REF], depending on whether it is negligible compared to the length of the needle. The width of the gauge, however, is not taken into consideration in the work of Bonvilain et al. [START_REF] Bonvilain | Characterization of strain microgauges for the monitoring of the deformations of a medical needle during its insertion in human tissues[END_REF] and Hammond et al. [START_REF] Hammond | Estimating surgical needle deflection with printed strain gauges[END_REF]. This paper presents a model for U-shaped strain gauges which takes into account the incurved width of the gauge on the needle surface. Relations are established between the strain calculated from the gauge and the axial strain at the surface of the needle at the gauge location. Needle axial strain expressions calculated with and without the proposed model are compared and instrumented numerical characteristics from the works of Bonvilain et al. and Hammond et al. are used to quantify the impact of strain gauge width consideration on needle axial strain estimation.

A brief review of the literature is proposed in Section II. The model proposed in this paper is presented in Section III with the geometrical modeling of the strain gauge and the establishment of the strain calculated from the gauge in terms of axial strain of the needle. The results of this new model are compared with current results in Section IV and are discussed in Section V. Finally, conclusions on the presented model are summarized in Section VI.

II. PRIOR WORKS

The Infinitesimal Width gauge model (IW) denotes the modeling where the width of the gauge, characterized by the trace separation and the trace width, is ignored, such as in Bonvilain et al. and Hammond et al. In that case, the strain indicated by a gauge is assumed to be equal to the length variation of the gauge divided by its initial length. For a gauge located between endpoints s a and s b whose initial length is l 0 = s bs a , as shown in Fig. 1, this gives the following relationship [START_REF] Hammond | Estimating surgical needle deflection with printed strain gauges[END_REF]:

ε gauge = ∆l l 0 (2) = 1 l 0 s b s a ε(s)ds ( 3
)
where ε is the axial strain at the surface of the needle.

III. FINITE WIDTH GAUGE MODEL

This section proposes a geometrical modeling of a Ushaped strain gauge bonded on a needle parallel to its length taking into account the width of the gauge. The geometrical characteristics of the model are defined from the cross sections of the needle with the U-shaped strain gauge, as illustrated in Fig. 2. The outer radius of the needle is denoted Fig. 2: Section of the needle with the U-shaped strain gauge. The center of the needle cross section is denoted C. The outer radius of the needle is denoted r. The location angle of the gauge α is measured from the center of symmetry of the gauge. The variables associated with the gauge are the trace separation l s defined by the semi-angle ϕ and the trace width l w defined by the angle ψ. The gauge angle range is 2(ϕ + ψ).

r. The gauge is located at the angle α. The trace separation of the strain gauge is defined by its arc l s and its associated semi-angle ϕ. The trace width is defined by its arc l w and its associated angle ψ. The angles ϕ and ψ are determined by the following relationships:

ϕ = l s 2r , ψ = l w r (4)
We Under the hypothesis of linear bending deformation, the axial strain at the surface of the needle ε(s, α) located at length s and angle α can be expressed with the curvature κ(s) and the direction of bending θ (s):

ε(s, α) = rκ(s) cos(θ (s) -α) (5)
Equation [START_REF] Roesthuis | Three-dimensional needle shape reconstruction using an array of fiber bragg grating sensors[END_REF] shows that ε varies depending on the angle α. This means that the axial strain on the circumference of a cross-section is not constant where the needle is deformed. According to Fig. 2 the strain gauge covers an angle range of 2(ϕ + ψ). Thus, the axial strain of the surface on which the gauge is located is not constant during deformations.

The problem of the strain gauge response on a non-uniform strain field was addressed by Schajer [START_REF] Schajer | Use of displacement data to calculate strain gauge response in non-uniform strain fields[END_REF]. In his work, Schajer presents ε gauge the strain indicated by a strain gauge grid as the average of the axial strain field on each of the grid lines of the strain gauge. Thus, for a strain gauge with n grid lines of respective area A grid i we have:

ε gauge = 1 n n ∑ i=1 1 A grid i grid i εdS (6)
As a U-shaped gauge is a gauge grid composed of 2 grid lines (the junction between them is not taken into account), the relation ( 6) using the FW model becomes:

ε gauge = 1 2 × 1 l 0 s b s a 1 l w α-ϕ α-ϕ-ψ ε(s, φ )rdφ + α+ϕ+ψ α+ϕ ε(s, φ )rdφ ds (7)
The complete calculation of ( 7) is presented in Fig. 3 for convenience. Finally, the expression of the strain indicated by the gauge ε gauge with the FW model is:

ε gauge = F(ϕ, ψ) × 1 l 0 s b s a ε(s)ds (17)
where

F(ϕ, ψ) = cos ϕ + ψ 2 sinc ψ 2 (18)

IV. STRAIN ESTIMATE COMPARISON

The expressions of ε gauge in ( 3) and ( 17) are directly proportional, with the function F(ϕ, ψ) defined in (18) as the coefficient of proportionality. 1 The estimates of the axial strain ε with the IW and FW models are denoted respectively [ε] IW and [ε] FW . From relations (3) and (17) it results that:

[ε] IW = F(ϕ, ψ)[ε] FW (19)
Equation ( 19) shows that the axial strain estimate with the IW model is equal to the axial strain estimate with the

ε gauge = 1 2 × 1 l 0 s b s a 1 l w α-ϕ α-ϕ-ψ ε(s, φ )rdφ + α+ϕ+ψ α+ϕ ε(s, φ )rdφ ds (8) = 1 l 0 s b s a 1 2ψ α-ϕ α-ϕ-ψ ε(s, φ )dφ + α+ϕ+ψ α+ϕ ε(s, φ )dφ ds (9) = 1 l 0 s b s a rκ(s) 2ψ α-ϕ α-ϕ-ψ cos(θ (s) -φ )dφ + α+ϕ+ψ α+ϕ cos(θ (s) -φ )dφ ds (10) = 1 l 0 s b s a rκ(s) 2ψ [sin(φ -θ (s))] α-ϕ α-ϕ-ψ + [sin(φ -θ (s))] α+ϕ+ψ α+ϕ ds (11) = 1 l 0 s b s a rκ(s) 2ψ (sin(α -ϕ -θ (s)) -sin(α -ϕ -ψ -θ (s)) + sin(α + ϕ + ψ -θ (s)) -sin(α + ϕ -θ (s)))ds (12) = 1 l 0 s b s a rκ(s) 2ψ (-2 cos(θ (s) -α) sin(ϕ) + 2 cos(θ (s) -α) sin(ϕ + ψ)) ds (13) = sin(ϕ + ψ) -sin(ϕ) ψ 1 l 0 s b s a rκ(s) cos(θ (s) -α)ds (14) = 2 cos ϕ + ψ 2 sin ψ 2 ψ 1 l 0 s b s a ε(s)ds (15) = cos ϕ + ψ 2 sinc ψ 2 × 1 l 0 s b s a ε(s)ds (16)
Fig. 3: Expression of strain response ε gauge in terms of axial strain ε for a strain gauge located between positions s a and s b on the needle at angle α.

FW model multiplied by the coefficient F(ϕ, ψ). Therefore, the difference between the axial strain estimates for IW and FW models depends exclusively of the gauge width parameters and is independant of the deformation. As the strain gauges modeled have an angle range smaller than half of the circumference of the needle, the following inequalities for F can be established:

0 < F(ϕ, ψ) < 1 ( 20
)
which associated with (19) gives:

|[ε] IW | < |[ε] FW | (21)
Thus, the amplitude of the strain estimate is inferior when the width of the U-shaped gauge is not considered. The relative difference of strain estimation between the FW and IW models expressed in percent, denoted η, is given as: Results show that when the gauge width is not taken into account the larger the strain gauge is, the larger the amplitude of the axial strain is underestimated. Thus, for the needle instrumented in Bonvilain et al. where the rangle angle of the strain gauge is 14 • the relative difference of strain is only of 0.4 % whereas it reaches 3.0% in the case of Hammond et al. where the range angle is 40 • .

η(ϕ, ψ) = [ε] FW -[ε] IW [ε] FW × 100% (22) = (1 -F(ϕ, ψ)) × 100% (
Therefore, the proposed model allows to estimate the impact of width modeling by quantifying the difference in strain estimation. It can then be used to correct the strain estimate value and then improve the accuracy. Moreover, as our model offers the possibility to take into account the length and the width of the gauge, the miniaturization of the gauge becomes less necessary as it is possible to correct the bias induced by both of these dimensions with our model.

VI. CONCLUSION

This paper introduced a new model for U-shaped gauge taking the width of the gauge into account. It is assumed that improving the strain gauge model is improving as well the accuracy of the strain estimation. Thus, the difference of strain estimation with the new model has been quantified and strain correction can be applied on instrumented needle from their characteristics.

Fig. 1 :

 1 Fig.1: Needle with a U-shaped strain gauge attached to its surface parallel to the needle length. The gauge is located between endpoints s a and s b and has an initial length of l 0 = s bs a . The trace separation of the gauge is l s and the trace width is l w .

  23) The values of the function F and the relative difference of strain estimation η for the needle instrumentations used by Bonvilain et al. and Hammond et al. are presented in Table I and Table II. These results show that the relative difference of the strain between the IW and FW model is 0.4 % for Bonvilain et al. and 3.0 % for Hammond et al.

Figure 4

 4 presents the plot of η(ϕ, ψ) according to the previously defined trace separation semi-angle ϕ and trace width angle ψ.

Fig. 4 :

 4 Fig. 4: Relative difference of strain estimation η(ϕ, ψ) between FW and IW gauge models according to angles ϕ and ψ. The angles ϕ and ψ are expressed in degrees. The values of η for Bonvilain et al. and Hammond et al. are plotted on the function graph of η.

  called this model Finite Width gauge model (FW) as the width of the gauge is taken into account. The values of the Finite Width gauge model parameters for the needles and strain gauges used by Bonvilain et al. and Hammond et al. are presented in Table I.

TABLE I :

 I Needle and U-shaped strain gauge characteristics for instrumented needles devices of the literature.

	Author	Needle radius r	Trace separation l s	Trace width l w	Trace separation semi-angle ϕ	Trace width angle ψ	F(ϕ, ψ)
	Bonvilain et al.	600 µm	50 µm	50 µm	2.39 •	4.77 •	0.996
	Hammond et al.	635 µm	150 µm	150 µm	6.77 •	13.53 •	0.970

TABLE II :

 II Effect of the gauge width consideration on strain estimation. Numerical values of the relative strain difference of strain estimation η are presented for instrumented needles devices of the literature.

	Article	Gauge angle range Relative difference of strain estimation η(ϕ, ψ)
	Bonvilain et al.	14 •	0.4 %
	Hammond et al.	40 •	3.0 %
	indicated by the gauge. It is then believed that this expression
	coming from a more complete model increases the accuracy
	of the needle axial strain estimate.

If ϕ → 0 and ψ → 0 then the width of the gauge becomes infinitesimal and F(ϕ, ψ) → 1. As a consequence, (3) and (17) are identical and the FW model is equivalent to the IW model in this case.

VII. ACKNOWLEDGMENTS This work is part of the project GAME-D, financed by the French National Agency for Research (ref: ANR-12-TECS-0019) and supported by Laboratory of Excellence CAMI (ref: ANR-11-LABX-0004-01).

The authors would like to thank Cecilia Hughes for her English reviews and corrections.


P Villette

Eve Afonso

G Couval

A Levret

M Galan

C Tatard

J.-F Cosson

P Giraudoux

J.-F Giraudoux

Metabarcoding of Bacterial Pathogens in a Rodent Pest: Which Organ?

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

If we pool all organs within an animal, we find an average of 5.4 bacterial species per animal* *5.4 ± 0.7 (mean ± standard error), which is 2.9 ± 0.9 species (bootstrapped mean and 95% confidence interval) greater than the mean bacterial richness of the heart, the organ with the highest bacterial richness.

We found different bacteria in different organs

Spleen

  throughput DNA sequencing of the V4 region of the 16s rRNA bacterial gene using the Illumina MiSeq platform. Sequences were grouped into Operation Taxonomic Units (OTUs) using a 97% similarity threshold, and taxonomy was assigned using the Silva SSU Reference database.

  All five organs had similar bacterial species richness * *Mean bacterial species richness ranged from 1.5-2.5 species/organ, and did not differ significantly between organs (generalized linear mixed-effects model for average OTU richness within organs, with animal as a random effect. The model was fit using maximum likelihood, a Poisson distribution, and log-link function). Error bars are standard error. are necessary for detecting all the bacteria present in the animal* *Bootstrapped mean bacterial richness as organs are included in the richness calculation for each animal. Error bars are bootstrapped 95% confidence intervals (999 iterations for means and confidence intervals). Increasing mean richness with increasing nmber of organs indicates that organs differ in the bacterial they host, and that bacteria are not concentrated in one or two organs Species RichnessIn addition, all 5 organs are necessary to detect host-population-level differences in bacterial assemblages. The animals sampled here are from two different populations; we found that organ bacterial richness does not vary with location (generalized linear mixed-effects model for average OTU richness within organs, with animal and location as random effects, fit using maximum likelihood, a Poisson distribution, and log-link function), but differences in pooled, liver, and lung bacterial assemblages map to host populations (ordination and PERMANOVA analysis of Jaccard dissimilarities). Livers and lungs are prone to having no detectable bacteria, making them poor candidates for studying host-population level bacterial assemblages, thus pooling organs within animals is the "best" choice.

We collected 13 water voles from Franche-Comte, and extracted and sequenced bacterial DNA from their hearts, lungs, livers, spleens, and kidneys* We found 25 potentially pathogenic bacteria in our voles Bartonella Bacterial parasites that are transmitted by ticks and fleas. The causative agent of cat-scratch disease is a member of this genus Acinetobacter Common soil bacteria that can cause infections in already-sick animals Mycoplasma A diverse group, bacteria in this genus can infect mammals, birds, reptiles, fish, insects, and plants Other bacteria we found:

	These
	Avibacterium
	Pasteurellaceae
	Helicobacter
	Aerococcus
	Chryseobacterium
	Cornyebacterium
	Leptospira
	Peptococcus
	Sphingomonas
	Streptococcus
	Ureaplasma

are bacteria that we know can cause disease in animals, and may be causing disease in our voles. The 5 most prevalent bacteria we found: Treponema The causative agents of yaws and human and rabbit syphilis are members of this genus