H Hamdane
email: hedi.hamdane@ifsttar.fr

R Anderson

C Masson

M Llari

T Serre

Assessment methodology of Active Pedestrian Safety Systems: an estimation of safety impact

Keywords: Pedestrian accident, Active Safety System, accident reconstruction, multi-body simulation

Introduction

Devoid of any protection, pedestrians are highly vulnerable to road accidents against a vehicle. To enhance their protection, new safety-based technologies have been introduced in the vehicle market. These on-board systems are developed to prevent crashes from occurring or reduce their severity by reducing the impact speed. Several methods assessing these systems have been presented (e.g. [START_REF] Lindman | Benefit estimation model for pedestrian auto brake functionality[END_REF]). This research is focusing on assessing the benefit of Active Pedestrian Safety Systems (APSS's) for pedestrian injury mitigation. Researchers have established a relationship between impact severity and variations in speed impact [START_REF] Rosén | Pedestrian injury mitigation by autonomous braking[END_REF][START_REF] Anderson | Potential benefits of forward collision avoidance technology[END_REF]. This project is examining the effect of speed reduction on variations in impact conditions. Outlines of the assessment method are presented here illustrated with one example.

Methods

The first step consists of gathering a sample of real vehicle/pedestrian crashes provided by in-depth crash investigation. A considerable level of details is required to reconstruct numerically the pre-crash sequence including trajectories of the vehicle and pedestrian prior to the collision and the eventual obstacles. Each crash is modelled by representing the vehicle and pedestrian involved and the road environment. An APSS is then virtually represented by the parameters of the sensor and actuator. Once modelling has been set up, all the required components of each sub-model (crash environment, vehicle, pedestrian, and sensor and actuator technology) are implemented through a computational simulation and so interacting in a virtual environment identical to the real world crash scenario. This batch simulation provides a set of data displaying a new impact speed distribution. This distribution is estimated according to the actuation of the emergency braking manoeuvre. The last step is to estimate change in injury outcome using the HIC. These changes are calculated through the use of multi-body system (MBS) software, MADYMO ® . The real accident is firstly simulated in order to obtain the actual risk. Finally, by adding the effect of the emergency braking manoeuvre, the simulation enables to find out if the risk is reduced.

Results and Discussion

To show the possibilities of this assessment method, an accident case has been selected. In the original configuration of the crash, the HIC was 1645.6 for an impact speed of 37 km/h. With an APSS fitted in the vehicle, the impact speed is reduced to 8.9 km/h and the pedestrian head doesn't hit any part of the vehicle and hit the ground resulting in a very low HIC value of 28.2.

Conclusions

In this research, a formalized assessment methodology has been presented and illustrated with one case to forecast the safety benefits of APSS. This method is based on confronting these systems to real accident configurations through computational simulations. The safety impacts of these systems is then estimated by comparing injury outcome with and without the system enabled for the crash set.

Figure 1

 1 Figure 1 Scheme illustrating the crash scenario

Figure 2

 2 Figure 2 Crash simulation with and without APSS


Xiao Yi Zhou

Franziska Schmidt

Bernard Jacob

Extrapolation of traffic data for development of traffic load models: assessment of methods used during background works of Eurocode

INTRODUCTION

Transport demand is growing worldwide. Indeed the road freight transport has increased in Europe by 45.7% between 1995 and 2008 and this strong growth trend seems likely to keep continuing at the same rate of about 2.7% per year (European Commission 2010). It is important to ensure that the European transport network, which has been mainly designed using the national codes, can sustain this continuing growth in demand.

The currently used normal load model, LM1, in Eurocode 1, Part 2 was first calibrated with a two weeks heavy traffic dataset from Auxerre (A6 motorway, Paris to Lyon, France) in the late 1980s; it was then re-calibrated by O'Connor et al (2001) with several representative European traffic datasets recorded at France weigh-in-motion (WIM) sites. The increasing traffic load effect has been found as the augmentation of traffic flow and probability of multiple trucks on bridge simultaneous. Therefore, the load model needs to be periodically re-assessed by current traffic, because of the wide changes in traffic volume, composition of traffic, vehicle weights and sizes, to ensure a satisfactory safety level for the design of new bridges; and also the quality of WIM data has increased greatly in the last decade due to improved technologies and the development of specifications regulating accuracy levels. Accurate prediction of extreme load effects expected during the proposed or remaining lifetime of a structure is a key issue for the design or assessment of highway bridges.

This paper reviews several methods to obtain these extreme effects: the extrapolation methods (fitting distribution to the upper tail data and Rice's formula) implemented in the background study of the normal load model, LM1, in Eurocode 1 Part 2, and recently extensively used methods (block maxima method and peaks over threshold method). The methods are applied on WIM recorded traffic to extrapolate characteristic values of long return period. Comparisons are made on the extrapolated extreme gross vehicle weight, total load on lane, and induced traffic load effect. We compare the difference between various prediction methods and between past and recent recorded data in France.

This paper is organized as follows: Section 2 reviews the methods. Section 3 describes the WIM data used in this study. Section 4 compares the extrapolation of GVWs, TLLs and TLEs with different methods based on the three sets of WIM data and gives some comments and section 5 gives some conclusions.

Extrapolation of traffic data for development of traffic load models: assessment of methods used during background works of the Eurocode X. Y. Zhou, F. Schmidt & B. Jacob

IFSTTAR, France

ABSTRACT: This paper reviews several extrapolation methods (for example fitting a distribution to the upper tail of the data and Rice's formula) implemented in the background study of the normal load model, LM1, in Eurocode 1, Part 2. Recent extensively used extreme value theory (EVT) based methods (block maxima method and peaks over threshold method) in other domains such as meteorology (for flood levels, quantity of rain, etc) are also reviewed. The methods are applied on weigh-in-motion (WIM) data to extrapolate characteristic values of long return period. Two different types of comparisons are made on extrapolating extreme gross vehicle weight (GVW), the total load on lane (TLL) and the induced traffic load effect (TLE). The first focuses on differences between the various prediction methods, and the other on changes between past and recently recorded data in France. The results show that the EVT based method is the better one for extrapolation.

KEYWORDS: extrapolation (prediction) method, gross vehicle weight, total load on lane, traffic load effect, characteristic value.

In bridge engineering, the GVW of a heavy truck is particularly important to assess load effects on structures, as a single-heavy-truck combined with common trucks governs the traffic loading scenarios for short-to medium-span bridges; while multiple trucks on bridge simultaneously governs the traffic loading scenarios for medium-span to large span bridges. The extreme traffic load effects induced by these loading events should be obtained by extrapolating. Let X be a random variable with distribution F , and 1 ,, n xx , an identically and independent distributed sample from X . The maximum value over the "n-observation" period is 1 max , , Figure 2. This distribution has a well known limit, given in Figure 1 (Weibull, Gumbel, Fréchet). (1996) and [START_REF] Flint | Extreme traffic loads on road bridges and target values of their effects for code calibration[END_REF]

Basic of theory

The methods of statistical inference used in the literature to predict extreme traffic load effects are numerous. The following critique of the literature distinguishes four types of methods: those fitting distribution to the upper tail data, those extrapolating traffic load effects by Rice formula, those dealing with the generalized extreme value distribution based block maximum method and those fitting generalized Pareto distribution to peaks over threshold.

Fitting a Gaussian distribution to the upper tail

The distribution of n M can be derived exactly according to the parent distribution function F of the sample [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF]. However, raising an initial distribution function to a power leads to various types of extreme value distribution (see Figure 1), because of the tail of the initial distribution function. Only the upper tail of parent distribution function seriously contributes to the maximum distribution function [START_REF] Bailey | Basic principles and load models for the structural safety evaluation of existing road bridges[END_REF] which is obtained by raising F to a power (see Figure 2). [START_REF] Jacob | Methods for the prediction of extreme vehicular loads and load effects on bridges[END_REF] fits a half normal distribution and Gumbel distribution to the high tail in order to predict the characteristic load and load effects of longer return periods. To calibrate AASHTO and develop a new load model, [START_REF] Nowak | Effect of truck loads on bridges[END_REF] fits a normal distribution to the upper tail of the ratio of load effect to HS20 load effect to extrapolate the average 75-year maximum load effect. However, neither the plotting position method nor the rules of threshold selection are described. [START_REF] Fu | Truck loads and bridge capacity evaluation in China[END_REF] fit a Gumbel distribution to extrapolate load effects for long return periods.

There is a common agreement on using multimodal normal distribution to fit the curve of the GVW: O'Brien et al. ( 2010) fit normal distribution to the high upper tail, O'Connor & O'Brien (2005) use trimodal, while [START_REF] Mei | Bimodal renewal processes models of highway vehicle loads[END_REF] utilize bimodal to fit GVW data. Thus a normal distribution is adopted to fit the upper tail of GVW for prediction. Its cumulative distribution function is:

2 2 1 ; exp 2 2 x Fx (1)
where the parameter vector is , , the mean and standard deviation values, respectively. The Ryear return level, z, corresponding to exceeded probability of is:

1 1 z (2)
where -1 is the generalized inverse of the normal distribution function. [START_REF] Flint | Extreme traffic loads on road bridges and target values of their effects for code calibration[END_REF] considered Rice's formula to predict longer term return level of the load effect for the case of free traffic and congested traffic. To apply Rice's formula, the load effect should be assumed to Gaussian stationary process. [START_REF] Ditlevsen | Traffic loads on large bridges modeled as white-noise fields[END_REF] proved that a load effect from a long traffic load (several days) on a large bridge with slowly and smoothly varying influence lines could be modeled as a Gaussian random process.

Rice formula

Then, under the Gaussian hypothesis, the mean rate () vx of excedences of a level 0 x (number of times the level x is crossed with increasing values) during a reference period (3) where is the vector of parameters: m is the mean of the process, is its standard deviation and is the standard deviation of the differentiated process. Then the characteristic value for a long return period R can be extrapolated based on the definition of return period which is the mean period between two occurrences. Therefore the corresponding value is:

1 , zR (4)
The critical issue in using the Rice formula is to estimate the parameters. [START_REF] Cremona | Optimal extrapolation of traffic load effects[END_REF] provided an approach to obtain the optimal extrapolation of traffic load effect. By taking the logarithm of the Rice formula, the parameter estimation problem is simplified to the identification of the parameters of a second order polynomial function. The same as in fitting a distribution to the upper tail, only the upper tail of the level crossing histogram is important to extrapolate the R-year return level. [START_REF] Cremona | Optimal extrapolation of traffic load effects[END_REF] suggested determining the optimal number of class intervals of the level-crossing histogram to use for fitting the Rice formula via Kolmogorov test. Getachew (2005)

Block maxima method

As shown in Figure 2, the higher power to which we raise the initial distribution function, the less data that contributes to the maximum distribution function. In other words, very small discrepancies in the estimate of F can lead to substantial discrepancies for n F . Fu & You (2009) states that 500 million data entries are needed to accurately predict a 75 year maximum with 4000 average daily truck traffic (ADTT): it is obviously impossible to obtain accurate extreme values in that way.

However, statisticians have found that n F as- ymptotically tends to three classical extreme value distributions: Gumbel, Frechet and Weibull, which correspond to three different tail behaviors of the initial distribution function F, exponential, polynomial and bounded, see Figure 1 [START_REF] Kotz | Extreme value distributions: Theory and applications[END_REF][START_REF] Castillo | Extreme value and related models with applications in engineering and science[END_REF]. A unification of the three families of extreme value distribution into a single family known as generalized extreme value (GEV) distribution has been widely used in recent years to avoid choosing which of the three families is most appropriate for the data at hand. Moreover, reducing the power will increase the stability of prediction. So the block maximum method may be the best method that fits the GEV distribution to periodic, thus the GEV distribution is fitted to the daily maxima in this paper. Its cumulative distribution function is:

1/ ; exp 1 z Gz (5)
where the parameter vector is ,, , the shape, local and scale parameters, respectively. The R-year return level, z, corresponding to a probability of exceedance is:

1 log 1 z (6)

Peaks-over-threshold method

In hydrology, finance and wind engineering disciplines, it is extensively agreed that the way of using information is rather uneconomical and unreasonable in the block maxima method. The discarding of some of the largest observation in the block, given that only the block maximum is considered, represents a loss of information. Peaks-over-threshold (POT) method is an approach that can avoid deciding distribution type and efficiently using upper tail data. [START_REF] Crespo-Minguillon | A comprehensive traffic load model for bridge safety checking[END_REF] use POT approach to study weekly maximal traffic load effects. [START_REF] James | Analysis of traffic load effects on railway bridges[END_REF] analyzes traffic load effects on railway bridges by POT approach. Asymptotic theory reveals that the excesses over sufficiently high threshold have an approximate distribution following generalized Pareto (GP) family which is:

1/ ( ; ) 1 1 y Hy (7) where y x u is the excess of x over threshold u , and the parameter vector is , , the shape and scale parameters, respectively. The R-year return level, z, corresponding to m observations over R-year is:

1 u z u m ( 8
)
where u is the proportion of data over the threshold u . For this study two sets of WIM data are studied to extrapolate the extreme value of GVW and traffic load effect on a given bridge. One was recorded in the 1986, and the other consists of modern data that was recorded in 2010. In order to eliminate possible measurement errors, these data were scrubbed first before processing with cleaning rules (Sivakumar et al., 2010;Getachew 2005;[START_REF] Enright | Simulation of traffic loading on highway bridges[END_REF]), such as GVW greater than 3.5t, vehicle speed in the range of 35 km/h to 160 km/h etc..

5-axle

Proportion of vehicle types

The first set was recorded by the same WIM station on A6 highway at Auxerre in France, and this site has 4 lanes of traffic (2 in each direction). However, two lanes were recorded for one week from 26 th May to 2 nd June 1986. It contains 46049 trucks after filtration. The second set of data was from a piezo-ceramic weigh-in-motion system on the A9 motorway near Saint Jean de Védas, South-East of France in 2010. Only the upstream traffic lanes are recorded. In total 835468 trucks from January 2010 to May 2010 (GVW greater than 3.5 t) were recorded in the three lanes, with an average daily truck flow of 6217 trucks.

Tables 1 and2 show the traffic remarkably evolved from 1986 to 2010. The average daily truck traffic (ADTT) has more than doubled in these 24 years; the average yearly growth rate is about 4.5% that is larger than the average growth rate of 2.7%. The composition of traffic has also hugely changed as shown in Figure 3, the proportion of 5-axle truck increases by about 30%, while the proportion of 4axle truck decreases by about 20%, it is caused by an increase in the gross vehicle weight limit from 36 to 40t in 1992. 4-axle tandem trucks were replaced by 5-axle trimdem trucks.

EXTRAPOLATIONS AND COMMENTS

To re-assess the load model of Eurocode, three types of extrapolation have been done in this paper, which are the prediction of extreme GVWs, TLLs and TLEs. The prediction of GVW is based on single truck weight, the prediction of TLL is based on total loads on a single lane with length 50m, and the TLE is the induced mid-span moment of a simply supported structure. The influence lines are illustrated in Figure 4, and the time history of TLL and TLE are illustrated in Figure 5 (It is a double y-axis figure, on which the left y-axis is mid-span moment and the right is total traffic load on lane). As proved in the authors' paper [START_REF] Zhou | Assessing confidence intervals on extreme traffic loads[END_REF], stable prediction can be achieved based on weekday's traffic, all extrapolations are based on weekday WIM data.

In order to demonstrate the efficiency of extrapolation, three different R-year return levels have been extrapolated: 1000-year return for each data set, 24year return level for the data of 1986 and 1-year return level for the data of 2010. Because of the very short measured period, the number of daily maxima were not enough to ensure to obtain reasonable ex-trapolation, the hourly maxima has been used to apply block maxima method. The predictions are compared with those in previousbackground studies. The differences between characteristic values for the past and this study are given.

Gross vehicle weight

The GVWs were analyzed for the purpose of extrapolating the extreme traffic load and calibrating the traffic load model. Figure 6 shows the GVW histogram for these three sites' data. As shown, the modal has moved from right to left as time goes on, more details are given in Table 3, which shows summary statistics for the data. The mean GVW of all trucks has decreased from 31.9 t to 26.6 t; and the same phenomena has been found on events like daily maxima and top 5% trucks. While this phenomenon is in contradiction with the continuous trend of truck load increase, some explanations are: (i) a better precision of WIM sensors and stations (more efficient calibration), (ii) more effective enforcement above all since 2008 in France with the WIM based scanning [START_REF] Marchadour | Development and implementation of a WIM network for enforcement in France[END_REF], and (iii) two different sites. Using the previously mentioned prediction methods on GVW as the basis for comparing the impact of the models and the evolution of traffic, the R-year return levels of GVWs are given in Table 4. From these sources some general observations can be made:

Firstly, as expected, the extrapolated values increase from 1986 to 2010 as the traffic flow and the proportion of maxi-code 5-axle trucks increase.

Secondly, the simply fitting of a normal distribution to upper tail data almost provides much higher estimates than the other two EVT based methods. This is surely due to the selected normal distribution: the shape parameters of the two EVT based methods are negative, which means the underlying parent distribution has a bounded tail, while the normal shape has exponential infinite bound. This shows once again the efficiency of the EVT based methods. Finally, the comparison of characteristic values for 24-year return period of the data of 1986 and 1-year return period of the data of 2010 show that these values are quite similar for the same methods. This means that the extrapolation can predict well the future events, i.e. the return levels can be The GEV method brings about the Weibull distribution which is bounded (the bounds for Table 4 are 73.17 t for 1986 and 93.12 t for 2010). So, there is a discussion that must take place about whether the GVW is bounded or not. This discussion was made in the Joint Committee on Structural Safety (JCSS) in the 80s. Ditlevsen claimed that because the whole truck population is finite, the GVW distribution must be bounded. Jacob and Rackwitz inclined for a non upper bounded distribution, because the GVW possible values are infinite, and no upper limit may be chosen without truncating some abnormal loads, or imposing an arbitrary asymptotic value much too high. The question is not fully decided, and the choice should carefully take into account the available information and the consequences on the extrapolation.

Total load on a lane length

The objective of this section is to predict the maximum of the total load on a lane length for a typical medium span bridge of 50 m. As shown in Figure 5, process of total load on lane is a square-wave; therefore the Rice formula and fitting a Gaussian to upper tail are not a good choice. The block maxima method is chosen to perform the extrapolation. Table 5 summarizes the distribution parameters and characteristic values for different cases. The extrapolated values are compared with those of past studies, see Table 5. Indeed, in the last four columns of Table 5, several extrapolations of past studies are given: column "A6-1" gives the difference between our extrapolated values and those originally extrapolated by Calgaro in the background studies for Eurocode 1 (traffic data is from the A6 highway in 1986); column "A6-2" is the difference between our extrapolated values and those obtained by Alan O'Connor [START_REF] O'connor | Report of studies performed on the normal load model, LM1, ENV 1991-3 Traffic loads on bridges[END_REF] in his benchmark on extrapolation methods (the traffic is that of highway A6 in 1998). The last column, called "LM1 Design effects", gives the difference between our extrapolated values and the design effects of the load model 1 of Eurocode 1. Some comments can be made:

Firstly, stable predictions have been obtained due to the large size of the samples. Secondly, we note that the extrapolated values are lower than the effects of load model 1 of Eurocode 1, in any cases. If we consider the recent WIM data (of 2010), we see that the 1000year return level calculated with this data is between 40% and 50% inferior to that obtained with the load models of Eurocode 1, which is on the safe side.

Traffic load effect

The load effect considered in this part is the bending moment of a simple supported bridge of 50 m length. The partial conclusions that can be made based upon the results given in Table 6 are: First, the Rice's formula provides larger estimates than other methods, and it states that the total load on a lane length is bounded, which is reasonable as the traffic flow, distribution of GVW is stationary. Secondly, generally the characteristic load effects provided by modern data are similar to those for the 1980s' data.

Figure 1

 1 Figure 1 Comparison of extreme value distributions.

Figure 2

 2 Figure 2 Extreme value distribution PDF for varying n 2.1 Prediction methods used in the background studies of Eurocode In the background studies for the development of the Eurocode for bridge loading (EC1.2 2003), Bruls et

  adopted Cremona's approach to analyze traffic load effects on bridge induced by both surveyed and Monte Carlo simulated traffic data. & O'Brien (2005) compared the predicted extremes of simply supported moment for span lengths of 5, 10, 20, 50, 100 and 200 m by the Rice formula, classic Gumbel and Weibull extreme value distribution. There is approximately a 10% difference between Rice and others. It demonstrates the importance of appropriate selection of an extreme value prediction method.

  Figure 3 Proportion of vehicle types

Figure 4

 4 Figure 4 Influence lines used in calculation

Figure 5

 5 Figure 5 Time history of total load and induced load effect

  Figure 6 GVW histogram

  with the data of the 1980s as with the data of 2010.

Table 1

 1

	The basic information
	Name of road	A6	A9
	Time of recording	May, 1986	Jan. -May, 2010
	Site location	Auxerre	St Jean-de-Vedas
	No. of measured lane	2	3
	Record period (days)	7	138
	No. of trucks recorded	46049	835468
	Av. hourly flow per lane	107	259
	Av. daily flow per lane	2558	6217

Table 2

 2 Proportion of vehicle types (classified by no. of axles)

	No. of axles	No. of trucks	A6, 1986 Percent (%)	No. of trucks	A9, 2010 Percent (%)
	4	14306	31.07	82926	9.93
	5	19835	43.07	628709	75.25

Table 3

 3

	Gross vehicle weight statistics (in tons)
		A6,1986	A9,2010
	Max of all trucks	76.1	74
	Mean of all trucks	31.865	26.603
	COV of all trucks (%)	0.436	0.407
	Mean of daily maxima	62.120	57.543
	COV of daily maxima (%)	0.084	0.127
	Mean of second mode	41.387	36.505
	COV of second mode (%)	0.148	0.093

Table 4

 4 Extrapolated characteristic value of GVW (in tons)

						Return
			Distribution parameters	period
	Site	Method	Shape	Scale or standard deviation	Location or mean	(years) 1000 Short*
	A6 1986	Normal GPD GEV	-0.09 -0.23	8.71 2.57 4.54	28.28 53.44	79.99 74.37 78.42 75.56 72.51 71.79
	A9 2010	Normal GPD GEV	-0.20 -0.10	8.64 5.57 4.85	47.68 -44.62	100.3 94.69 81.61 75.76 82.12 72.41
	*24-year, 1-year return levels for A6-1986 and A9-2010 data
	respectively.

Table 5

 5 Extrapolated characteristic value of load on lane (in kNs) 1 : value extrapolated by Calgaro (quoting from O'Connor et al, 1998) * 2 : value extrapolated by O'Connor (quoting from O'Connor et al, 1998)Table 6 Extrapolated characteristic value of single lane's traffic load effect (in kN.m)

			Distribution parameter		Return period	Difference (%)
	Site	Method	Shape or standard deviation	Scale or stand-ard deviation of differentiate process	Location or mean	Short	1000 (years)	A6-1* 1	A6-2* 2	LM1 Design effects
	A6 1986	GEV (see 2.2.3)	-0.56	149.2	808	1072	1072	-35.58	-42.74	-45.03
	A9 2010	GEV	-0.31	126.1	661.6	1033	1055	-36.61	-43.64	-45.90
			Distribution parameter		Return period	Difference (%)
	Site	Method	shape	Scale	Location	Short	1000 (years)	A6-1	A6-2	LM1 Design effects
	A6 1986	Rice GPD GEV	-31973 0.20 -0.13	96892456 292.5 517.7	4470372.9 -5629	16262 14826 8829	19873 26385 9178	99.84 123.17 2.87	41.67 58.20 -27.08	26.12 40.84 -35.08
	A9 2010	Rice GPD GEV	10229.3 -0.13 -0.10	3821.8 591.7 599.2	2689502 -4764.5	8563 9128 8226	10229 9631 9455	-7.70 31.46 -4.92	-34.57 -6.81 -32.60	-41.75 -17.04 -39.99

*

Thirdly, among the three methods, the GP distribution based POT approach provides better extrapolations. Finally, the extrapolated characteristic values assessed based on short measurements are more prone to variations. This is the case for the return levels assessed by means of the traffic of 1986.

CONCLUSIONS

In this study, various prediction methods are applied to extrapolate the characteristic GVWs, TTLs and TLEs for long return periods. The results show that it is very important to choose the adapted prediction method for the considered effect: for example, POT is the optimal for bending at mid-span of a singlespan structure of 50m. In generally, the EVT based method is the better one in the studied cases. Moreover, a bigger size of sample gives more stable extrapolations, while Monte-Carlo simulation or the Bootstrapping method could be a useful method to reduce the variation of a sample by expanding the size of sample. The comparisons between extrapolations and LM1 design effects show that the normal load model LM1 is sufficient for modern traffic.

ACKNOWLEDGEMENT

The authors greatly acknowledge the financial support provided by the Marie Curie Initial Training Network TEAM (Training in European Asset Management) project


Antoine Ferlin
email: antoine.ferlin@ifsttar.fr

Philippe Bon
email: philippe.bon@ifsttar.fr

Simon Collart-Dutilleul
email: simon.collart-dutilleul@ifsttar.fr

Virginie Wiels
email: virginie.wiels@onera.fr

Parallel verification of temporal properties using dynamic analysis

come

I. INTRODUCTION

The development of critical software is constrained by certification standards. The standard depends on the application domain. For instance, DO-178 is the certification standard for avionics software. In railway transportation, IEC 50128 is the certification standard. Even if standards define goals for each step of the software development, companies have responsibility for choosing the methods to achieve these goals. Known methods are proposed by standards but are not mandatory.

The verification phase is one of the steps of the software life-cycle. Several kinds of verification methods exist and can be classified using two criteria: formal or not, static or dynamic. Formal verification means the use of computable mathematical rules based on a language with unambiguous grammar and a defined semantics. These rules allow a mathematical verification of a given property. Static verification means that the program to be verified is not executed. On the contrary, dynamic verification requires the execution of the program.

For instance, classic verification means are typically not formal: the test is dynamic whereas reviews are static. We can notice that a lot of formal techniques are static: B methods [START_REF] Wiels | The B-book: Assigning Programs to Meanings[END_REF], abstract interpretation [START_REF] Cousot | Abstract interpretation: a unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF], model check-ing [START_REF] Clarke | Automatic verification of finite-state concurrent systems using temporal logic specifications[END_REF]. . . RuntimeVerification1 , which is the purpose of this paper, can be classified as formal and dynamic. This work follows a PhD thesis done at AIRBUS and at ONÉRA, the French aerospace lab, about verification of temporal properties [START_REF] Ferlin | Combination of static and dynamic analyses for the certification of avionics software[END_REF], [START_REF] Ferlin | Vérification de propriétés temporelles sur des logiciels avioniques par analyse dynamique formelle[END_REF]. Several proprietary embedded programs have been studied in order to determine the properties which are difficult to be verified. But for reasons of industrial confidentiality, the detailed results cannot be published. According to this study, it appears that temporal properties are complex to be verified with the current industrial practices. No automatic verification method is currently used. Indeed, only code reviews are processed. The goal of the work was to propose a new equipped method to verify these temporal properties, in order to reduce time and cost of verification.

The resulting method is based on runtime verification. A specific language has been defined to formalize encountered temporal properties. This one is grounded on an aggregation of Linear Temporal Logic (LTL) and regular expressions, which are well adapted for sequence properties. The defined approach consists of two major steps: transforming the temporal property into a non-deterministic Büchi automaton, using Ltl2ba [START_REF] Gastin | Fast LTL to Büchi automata translation[END_REF], and then executing the Büchi automaton on an execution trace to verify the property. Because LTL has a semantics on infinite traces, the finite execution trace is transformed into an infinite one by looping over the last state of the trace. In order to counter the side effects, statistical information is computed to help to the interpretation of the results, for unclear cases. The topic of this paper is not the language definition, the verification phase or the finite trace problem. Readers interested in this topic will find this information in [START_REF] Ferlin | Vérification de propriétés temporelles sur des logiciels avioniques par analyse dynamique formelle[END_REF]. This approach has been reused in the railway context. In order to save time during the verification phase, this paper presents a parallelized version of this method. This new method is based on the divide-and-conquer strategy using parallel execution and the result should be predictable. But, this strategy has to be validated, because we cut a partial execution trace, and because the analysis of a sub-trace naturally depends on the previous sub-traces. In addition, the fusion operation of each sub-trace analysis requires the computation of additional features in comparison with the classic sequential method. We aim to test this approach on traces coming form an European Rail Traffic Management (ERTMS) proprietary simulator [START_REF]European rail software applications[END_REF]. However, because this platform is currently evolving with the PERFECT project2 , we first propose to do experimentation on a generated trace.

After this introduction, section II discusses this approach with regard to the related works and the industrial context. Section III summarizes the approach defined in [START_REF] Ferlin | Vérification de propriétés temporelles sur des logiciels avioniques par analyse dynamique formelle[END_REF]. Section IV defines some notations in this paper, then section V formalizes the classic approach which consists in executing a Büchi automaton on an execution trace and section VI formalizes the parallelized approach. Then, section VII presents experiments and comparison between the two methods. Finally, section VIII concludes this article discussing efficiency of the parallelized approach versus the sequential approach, and identifying the current limitations of the parallel method.

II. CONTEXT

A. Related Work

The existing works on runtime verification can be classified in two categories: online methods which perform the verification during execution of the program and a-posteriori methods, which perform verification on execution traces. There is a lot of work on online verification, especially for Java programs [START_REF] Drusinsky | The temporal rover and the atg rover[END_REF], [START_REF] Patrick | An overview of the mop runtime verification framework[END_REF], [START_REF] Havelund | A rewriting-based approach to trace analysis[END_REF], [START_REF] Havelund | Monitoring programs using rewriting[END_REF] and in the aspect-oriented programming community [START_REF] Stolz | Temporal assertions using aspectj[END_REF]. Existing works are based on rewriting techniques [START_REF] Havelund | Monitoring programs using rewriting[END_REF] or specific techniques, such as translation of LTL formula into state machines [START_REF] Amorim | Efficient monitoring of omegalanguages[END_REF], [START_REF] Giannakopoulou | Automata-based verification of temporal properties on running programs[END_REF].

Reduction of verification time is a crucial issue in allowing us to deal with industrial software. In online verification, many efforts are made in this way. [START_REF] Bodden | Efficient hybrid typestate analysis by determining continuation-equivalent states[END_REF] proposes to reduce the number of monitors, using static analysis. [START_REF] Berkovich | Gpu-based runtime verification[END_REF] proposes a parallel verification of several three-valued-LTL properties encoded in a monitoring system using the GPU3 . [START_REF] Pellizzoni | Hardware runtime monitoring for dependable cots-based real-time embedded systems[END_REF], [START_REF] Haitao Zhu | Predictable runtime monitoring[END_REF], [START_REF] Zilles | A programmable co-processor for profiling[END_REF] separates monitors from the program to be verified.

The approach of this paper is an offline approach; execution traces come from a railway simulation framework. Verifications are done on execution traces rather than online, for several reasons. One of the most important is that the verification phase must not disturb the execution of the program to be verified. Actually, for real-time programs, a reduction of the execution speed can change the veracity of a property. In this context, there are three factors depending on the trace and affecting the verification time: the size of the trace, the number of variables inside the trace, and the format of the trace. The last factor is not optimizable because the simulator is proprietary.

In several existing works, execution traces are obtained by listening to all variables of the program to be verified, even if variables are not necessary to verify the specified properties [START_REF] Barringer | Eagle does space efficient ltl monitoring[END_REF], [START_REF] Pnueli | Psl model checking and run-time verification via testers[END_REF]. These approaches lead to generation of big traces, because some variables and states are unnecessary.

Static analysis may be used to minimize the size of execution traces [START_REF] Ferlin | Combination of static and dynamic analyses for the certification of avionics software[END_REF]. Static analysis detects all program points where variables relevent for the verified property are modified. An observation point is defined by collecting instruction at each needed location. But reduction of trace is not sufficient if the simulation goes on for several hours because some variables perpetually evolve and hence the size of the trace dramatically grows.

Consequently, the verification must be improved. [START_REF] Geimer | Scalable parallel trace-based performance analysis[END_REF] proposes an approach to verify several traces coming from a parallel application. Properties are verified using parallel processes (one process for each trace) instead of merging all traces into a unique trace to be analysed. This approach allows the verification of traces coming from parallel programs. If splitting the generation trace is not possible, another solution consists in splitting the trace after its generation and verifying the property on each piece of trace. This is our proposition.

B. Railway Transport context

Nowadays, each new line built inside the European Union must comply with both national rules and ERTMS/ETCS (European Railway Transport Management System/European Train Control System). ERTMS/ETCS specification is a European proposition about embedded systems and communication between trains and infrastructures. This new technological and legal context must be implemented in different European states. Indeed, line operation rules must comply with the ERTMS/ETCS norm.

In this context, the considered software are critical and their verification is essential for avoiding collisions or nearcollisions [START_REF] Aguirre | Application of evidential networks in quantitative analysis of railway accidents[END_REF]. To do this, a proprietary simulator allows the simulation of the train behaviours. This proprietary simulator virtualizes control centres, railway with stations, markers, a train with control command and the communication with the control center. A train can be driven by an operator as a real train or by a defined scenario. We propose to adapt and to parallelize the approach of [START_REF] Ferlin | Vérification de propriétés temporelles sur des logiciels avioniques par analyse dynamique formelle[END_REF], in order to analyse temporal properties on execution trace of the simulator.

III. VERIFICATION APPROACH

In this section, the approach defined in [START_REF] Ferlin | Vérification de propriétés temporelles sur des logiciels avioniques par analyse dynamique formelle[END_REF] is summarized. The goal of this work has been to verify temporal properties on avionics embedded software.

A. Language

Instead of using a complex language, we chose to define a language which is well adapted to our needs. Hence, the first step of our work was to study the documentation of a lot of avionics software to extract temporal properties. For reasons of confidentiality, this study has never been published.

After extracting temporal properties, they were classified according to an extension of Dwyer's classification [START_REF] Dwyer | Property specification patterns for finite-state verification[END_REF]. This expansion was necessary because frequency properties, properties with time interval, were not present in the original classification. The goal of this classification was to determine what kind of logic and operators were necessary.

Finally, the defined language is a combination of LTL (Linear Temporal Logic), regular expressions, and some additional operators for properties on numerical variables (integer, float): comparison between two numerical numbers, addition, substraction, multiplication, division. Operators on variables have been defined to change the variable scope: if x is a variable

• x?T refers to the last time when x was modified,

• x?C refers to the number of changes the variable suffered,

• x?n refers to the value of the variable n changes ago,

• x?n?T refers to the time when x was modified n changes ago,

• x?n?C refers to the number of modifications n changes ago.

B. Trace Generation

In order to generate the trace, a dynamic analysis platform was used. This virtualizes the hardware of the software to be verified. For efficiency reasons, it was not possible to listen to a variable of the program. Observation points are used instead. As a consequence, in order to correctly verify a property using a set V of variables, it was necessary to set an observation point each time a variable of V changed.

Static analysis was used to find all necessary observation points. A plugin for Frama-C [START_REF]The frama-c platform for static analysis of c programs[END_REF], called Breakpointer, was implemented to do this task. It was based on the Frama-C plugin Value Analysis which performs a semantic analysis. Hence, pointers on variables are taken into account.

Instead of using this avionics simulator, the ERTMS simulator will be used to generate execution traces. Each trace is stored inside a database. Each trace maps to one table with a timestamp, a type of message, and the associated values. The type of message describes which variable is modified.

C. Verification

The verification phase consists in:

• Transforming the temporal property into a Büchi automaton using Ltl2ba

• Verifying the property by executing the generated Büchi automaton on the trace execution.

If the temporal property follows a pattern, a statistical automaton is generated and executed concurrently to the Büchi automaton. This automaton is a state-transition machine, each transition contains a propositional logic formula to activate the transition and a list of counter assignments to change the value of defined global counters representing interesting information defined for each property pattern.

IV. RECALL

A. Notations

In the following sections, some notations are used to help formalization of both approaches. The meaning of these notations is described hereafter:

• P(Q) is the set of subsets of a set Q • the set E is a Cartesian product such as E = E 1 × E 2 × ... × E n , then E |i is projection of E on E i . In other words, E |i = E i . If E i is a Cartesian product itself such as E i = E i,1 × ... × E i,m , hence E |i,j
refers to E i,j . This notation is applicable on function ; f :

X 1 × ... × X n → Y 1 × ... × Y m is a function. f |i refers to projection of f on Y i . • if E is a set, E ω is the infinite Cartesian product of E.
• a; b such as (a, b) ∈ N 2 and a < b is the set of integers {a, a + 1, ..., b}.

• a trace σ is a sequence of states. A state σ i ∈ Σ, where Σ is the alphabet, at index i ∈ 0; |σ| -1 of the trace σ is a total function which returns the value of a given variable. |σ| is the length of the trace σ. In this paper, a trace is considered as a word based on Σ alphabet in the meaning of formal language theory [START_REF]Handbook of Formal Languages[END_REF].

• if K is a set of values, then K ¿ = K ∪ {¿},
where ¿ is the unknown value;

In addition, item hereafter, to avoid confusions, we will speak about element to define a trace state and about state to define a automaton state.

B. Büchi automaton definition [8]

In this section, formal definitions of an automaton [START_REF] S C Kleene | Representation of events in nerve nets and finite automata[END_REF], of a Büchi automaton and of a statistical Büchi automaton are recalled. After, we schematized an execution of a Büchi automaton.

1) Büchi automaton:

Definition 1: An automaton is classically defined as a fiveuplet A = (Q, Σ, →, q 0 , F ) where:

• Q is a set of states • Σ an alphabet • →⊆ Q × Σ × Q is a transition relation • q 0 ∈ Q is the initial state; • F ⊆ Q is a set of accepting states.
The accepting condition of a finite word by an automaton is: Definition 2: A finite word w ∈ Σ * is recognized by an automaton A = (Q, Σ, →, q 0 , F ), if and only if there is a sequence (q) i∈ 0;|w| , which starts with the initial state q 0 , such that for all i ∈ 0; n -1 , (q i , w i , q i+1 ) ∈→, such that q n ∈ F .

A word w which is recognized by A is written w ∈ L(A), where L(A) is the set of words recognized by A, A Büchi automaton is a classic automaton with a special accepting condition which allows the handling of infinite traces. The accepting condition of a trace is: Definition 3: An infinite word w ∈ Σ ω is a word inside L(B), where B = (Q, Σ, →, q 0 , F ) is a Büchi automaton, if and only if:

• there is a sequence (q) i∈N such that ∀i ∈ N, (q i , w i , q i+1 ) ∈→

• ∀j ∈ N ∃k ∈ N such that k > j and q k ∈ F .

In this paper, a property has to be verified on an execution trace. Hence, the corresponding Büchi automaton is executed on this trace. Consequently, the alphabet used is built from the trace states.

In order to compute statistical information on a given infinite word, a statistical Büchi automaton is defined as an extension of a Büchi automaton. When a formula of a given transition is true, then primary operations, defined below, are done on some given counters. At the end of the execution of the statistical automaton, counters quantify properties which are orthogonal to the temporal property. The statistical Büchi automaton can be deterministic or not. Actually, a statistical automaton is not an automaton used to analyse a program whose the behaviour is not determined as in [START_REF] Bulychev | Monitorbased statistical model checking for weighted metric temporal logic[END_REF]. The word statistical only refers to the statistical information which are computed.

Example 1: For instance, the property to be verified is (♦ e), which means e infinitely often occurs. The number of occurrences of e inside a given trace is statistical information.

Contrary to the work in [START_REF] Ferlin | Vérification de propriétés temporelles sur des logiciels avioniques par analyse dynamique formelle[END_REF], in this article, the statistical Büchi automaton replaces the classic Büchi automaton when a temporal property maps to the pattern of the statistical automaton. Before defining a statistical Büchi automaton, we define statistical operations.

Definition 4: C is a set of integer variables which will be called counters.

Λ C : (C → Z ¿ ) → (C → Z ¿
) is a set of action on C, depending on current value of all variables. Operations can be:

• doing nothing,

• assigning a constant or another counter value, or an expression which can be:

• addition, subtraction, multiplication, division of counter/constant • minimum, maximum of counter/constant In the case of parallel execution of a statistical deterministic Büchi automaton, the unknown value (¿) is used and allows the symbolic computation of the value of counters.

As a definition, the statistical operation λ of a given transition is a list of actions.

Finally, the statistical Büchi automaton is defined as follows:

Definition 5: A statistical automaton is a five-uplet A = (Q, Σ, →, q 0 , F, C 0 ) where: • Q is a set of states • Σ an alphabet • →⊆ Q × Σ × Λ C × Q a transition relation • q 0 ∈ Q is the initial state; • F ⊆ Q a set of accepting states; • C 0 : C → Z ¿ ,

is a function which returns initial value of each element of C

Hereafter, the term word is replaced by execution trace.

2) Execution of a Büchi automaton: We use three cases of Büchi automaton: the deterministic one, the statistical deterministic one and the non-deterministic one. Before giving a formal definition of an execution of a Büchi automaton for each case, we illustrate these three cases with examples.

We speak about a step of computation when the Büchi automaton consumes an element in order to enable accessible transitions. When the Büchi automaton is deterministic (figure 1), at each step of computation, there is only one state before the consumption of the mapping element (current element) and at most one state after the consumption of the element. Hence, the execution of the Büchi automaton can be seen as a path with several nodes. Each node is a computation step and there is only one link between two nodes. The execution of a statistical deterministic Büchi automaton corresponds to the execution of a deterministic Büchi automaton with a set of counters with their value for each current state. In figure 2, counters with their value are modelled by C i .

A B B A C σ 0 σ 1 σ 2 σ 3 σ 4 σ 5 computation step
A B B A C σ 0 σ 1 σ 2 σ 3 σ 4 σ 5 computation step C 0 C 1 C 2 C 3 C 4
Fig. 2: Execution of a statistical deterministic Büchi automaton

If the Büchi automaton is not deterministic, then there will be a set of current states instead of one state. The path becomes a lattice instead of a sequence of automaton states. In figure 3, each rectangle is a set of current states. Rectangle of level 1 is the initial state, the set {A}. One of level 2 is the set {A, B} after consumption of the first element. And so on... Schematically, a computation step consists in going from a rectangle of level i to the nearby one (i + 1). Hence, the formal definition of an execution of a program has to contain two sets of state(s). The following sections take this fact into account.

A B A B B A A C A σ 0 σ 1 σ 2 σ 3 σ 4 σ 5 level 1 {A} level 2 {A, B} level 3 {A, B} level
The execution of a statistical non-deterministic Büchi automaton is not described. Consequently, in the rest of the article, we will speak about statistical automaton for a Büchi statistical deterministic automaton.

V. CLASSICAL EXECUTION OF BÜCHI AUTOMATON ON

TRACES

Now, we formalize a sequential execution of a nondeterministic Büchi automaton and a statistical Büchi automaton.

A. Execution of a non-deterministic automaton on trace

The formalization of the execution of a non-deterministic Büchi automaton is a core material for the other cases.

In the approach in [START_REF] Ferlin | Vérification de propriétés temporelles sur des logiciels avioniques par analyse dynamique formelle[END_REF], regular expression are translated into deterministic Büchi automaton by an AIRBUSproprietary tool, whereas LTL formulas are translated into nondeterministic Büchi automaton by Ltl2ba.

As a reminder, the formal definition 3 requires a sequence of states, hence the definition is sufficient for describing the first case. On the contrary, generated Büchi automata are non-deterministic in the second case. To deal with this case, definition 3 is sufficient. However, the complete formalization requires the definition of a sequence of set of states instead of a sequence of states.

Formalization of an execution of a non-deterministic Büchi automaton on a trace is based on figure 3. An automaton execution is a sequence of computation steps. Each computation step is a pair of state sets: one set before the trace state consumption and one set after. Formally, the definition of an execution of a Büchi automaton on a trace is: Definition 6: An infinite word w ∈ Σ ω is a word inside L(B nd ), where B nd = (Q, Σ, →, q 0 , F ) is a non-deterministic Büchi automaton, if and only if there is a sequence R i∈N ∈ P(Q) 2 such that:

• R 0 = ({q 0 }, {q 0 }) (1)
• ∀i ∈ N * , R i ∈ P(Q) 2 , such that:

• R i|1 = R i-1|2 (2) • ∀r i ∈ R i|2 , ∃r i-1 ∈ R i|1 , such that (r i-1 , w i , r i ) ∈→ (3)
• ∀j ∈ N, ∃k ∈ N, k > j and ∃R k|1 ∩ F = ∅. (4)
Property 2 ensures the consistency of the automaton execution. Property 3 is the part of the definition of an accepted word, which handles the existing path, for a given Büchi automaton. Property 4 is the other part of the definition of an accepted word, which requires an accepting state which is infinitely often reached.

If an infinite word w is not recognized by a nondeterministic Büchi automaton, then ∃k such that R k|2 = ∅ and for all k > k, R k = (∅, ∅), or else the sequence of the set of the current automaton states contains a finite number of accepting states.

In figure 3, we obtain the following sequence:

R 0 = ({A}, {A}) R 1 = ({A}, {A, B}) R 2 = ({A, B}, {A, B}) R 3 = ({A, B}, {A, B}) R 4 = ({A, B}, {A, C})
In our industrial railway context, handled finite traces are transformed into infinite ones by looping over the last state of the trace. This choice allows the use of existing tools such as Ltl2ba without modification. The finite trace problem is not the purpose of this article and was broached in [START_REF] Ferlin | Vérification de propriétés temporelles sur des logiciels avioniques par analyse dynamique formelle[END_REF].

The limit of the classic verification of a trace with a Büchi automaton is that the sequence R i∈ 0;|σ|-1 is sequentially computed. Hence, it slows the verification process.

B. Execution of a deterministic statistical Büchi automaton on trace

In this section, an execution of a statistical deterministic Büchi automaton on a trace is formalized. Definition 7: An infinite word w ∈ Σ ω is a word of L(B S ), where B S = (Q, Σ, →, q 0 , F, C 0 ), if and only if there is a sequence

R S i∈N ∈ (Q) 2 × (C -→ Z ¿ ) 2 such that : • R S 0 = (q 0 , q 0 , C 0 , C 0 ) (5)
• ∀n ∈ N * , R S n ∈ (Q) 2 × (C -→ Z ¿ ) 2 such that: • R S n|1 = R S n-1|2 et R S n|3 = R S n-1|4 (6)
• (R S n|1 , w i , λ n , R S n|2 ) ∈→ (7) • R S n|4 = λ n (R S n|3 ) (8)
• ∀j ∈ N, ∃k ∈ N, k > j and R S k|1 ∈ F . (9) Property 6 ensures consistency of automaton execution. Property 7 is a part of the definition of an accepted word. Property 9 is the other pat of the definition of an accepted word, which requires an infinitely-often-reached accepting states. In property 7, w n is the current element and λ n is the current operation applied to the counters (property 8).

If we deal with figure 2, then the sequence will be:

R S 0 = (A, A, C 0 , C 0 ) R S 1 = (A, B, C 0 , C 1 ) R S 2 = (B, B, C 1 , C 2 ) R S 3 = (B, A, C 2 , C 3 ) R S 4 = (A, C, C 3 , C 4 )
Despite the fact that trace analysis naturally looks like a single sequential process, we will show that parallel execution of the Büchi automaton improves the computing efficiency.

VI. PARALLEL EXECUTION OF BÜCHI AUTOMATON ON FINITE TRACES

A. Principle

Until now, when a trace has been verified by execution of a Büchi automaton, it has been necessary to compute R 0 , ..., R i-1 , before computing R i , because of the equality R i-1|2 = R i|1 . We propose the following new approach:

1) splitting the execution trace in several pieces, 2) executing a Büchi automaton B on each piece of the trace 3) merging the results of each execution of B to determine if the trace belongs to L(B).

Figure 4 synthesises the approach on a trace which is split into two pieces, at trace state 3. Two identical Büchi automata are executed on traces σ 0 ..σ 2 and σ 3 ...σ n . The execution of the Büchi automaton on the first trace follows the rules defined in the previous section. The execution of the Büchi automaton on the second trace begins with Q ({A, B, C}) as the set of initial states.

The verification of the property on the entire trace is performed using a function which acts as a short-cut between the beginning and the end of each part of the trace.

For instance, if the execution trace stops after trace state 5, we can see in Figure 4 that:

• state A before element 3 leads to nothing after element 5.

• states B and C before element 3 lead to A and C after element 5.

Because at level 3 the set of current states is {A, B}, the set of states after trace state 5 is {A, C}.

After presenting the principle, we will formalize it on nondeterministic and statistical deterministic Büchi automata. For reasons of clarity, the formalization will be done on a trace which is split into two sub-traces, but this one is generalizable to several divisions.

A B A B A A B A B C A C level 1 {A} level 2 {A, B} level 3 {A, B} level 3 Q level 4 {A, B} level 5 {A, C} σ 0 σ 1 σ 2 σ 3 σ 3 σ 4 σ 5 σ 6 (a) non-deterministic A B B A B A C B A C level 1 {A} C 0 level 2 {A, B} C 1 level 3 {A, B} C 2 level 3 Q I 0 level 4 {A, B} I 1 level 5 {A, C} I 2 σ 0 σ 1 σ 2 σ 3 σ 3 σ 4 σ 5 σ 6 (b) statistical
Fig. 4: Principle of parallel approach

B. Execution of a non deterministic Büchi automaton

The trace σ is split into two sub-traces at element index c. The Büchi automaton will be executed on the both sub-traces. To do so, two sequences of computation steps will be defined: R n∈ 0;c and R n∈ 0;|σ|-c+1 . R n∈ 0;c is the sequence of computation steps which begins at element σ 0 as defined in section V. Hence, for all n ∈ 0; c ,

R n = R n .
The sequence of computation steps on the second trace is defined as:

Definition 8: Sequence R n∈ 0;|σ|-c+1 ∈ (P(Q)) 2 × (Q -→ P(Q)) 2 is such that :
• Definition of the first two components of R n as following:

• R 0|1 = Q (10) • ∀n ∈ 0; |σ| + 1 -c , R n|1 = R n-1|2 . ( 11
) • ∀n ∈ 0; |σ|+1-c , ∀r n ∈ R n|2 , ∃r n ∈ R n|1 , such that (r n , σ c+n+1 , r n ) ∈→ (12)
• The definition of the last two components of R n corresponds to the α function:

• ∀n ∈ 0; |σ| + 1 -c , R n|3 = α n , (13)
• ∀n ∈ 1; |σ| + 1 -c , R n|4 = R n+1|3 , ( 14
) ∀q ∈ Q, α 0 (q) = {q}. (15) α n (q) = {q |∃q ∈ α n-1 (q), (q σ c+n+1 , q ) ∈→}. ( 16
)
Concretely, R is the sequence of the execution of the Büchi automaton which begins at element index c and ends at element index |σ|-1. The initial set of states is Q. Function α n records the link between the set of states at element c + 1 and the set of states at element index |σ| -1.

Property 10 means that the sequence of computation steps is initialized with Q as the initial set of states. Properties 11 and 12 respectively ensure the consistency of the execution and th acceptance of a trace by the automaton. α is initialized with the identity function (property 15). For each state, α is recursively computed depending on activated transitions (property 16).

In figure 4a, the sequence will be:

R 0 = ({A}, {A}) R 1 = ({A}, {A, B}) R 2 = ({A, B}, {A, B}) R 0 = ({A, B, C}, {A, B}, α = id, {α(A) = ∅, α(B) = {A, B}, α(C) = {A, B}}) R 1 = ({A, B}, {A, C}, {α(A) = ∅, α(B) = {A, B}, α(C) = {A, B}}, {α(A) = ∅, α(B) = {A, C}, α(C) = {A, C}})
Execution of the Büchi automaton on the first sub-trace can be performed independently of the execution on the second one. The results of the analysis of each sub-trace have to be merged using the α function.

Using the same approach, it is possible to divide the trace into more than two sub-traces. Actually, the prefix of the trace is classically analysed and we use the definition of R for each other sub-trace.

C. Execution of a statistical deterministic Büchi automaton

Let us recall some notations: σ is an execution trace, B S = (Q, Σ, →, q 0 , F, C 0 ) is a statistical deterministic Büchi automaton.

The execution sequences look like R defined in section VI-B. An additional function is computed as α function, in order to built the operation applied on counters between the beginning and the end of each piece of trace.

R S

n∈ 0;c is the execution sequence which begins at trace state σ 0 as defined in section V-B. Hence for all n ∈ 0; c , R S n = R S n . The other sequence of computation steps is built as the previous-section one, taking into account the computation of statistical information. 2 , where, for all n ∈ 0; |σ| + 1 -c :

Definition 9: R S n∈ 0;|σ|+1-c ∈ (P(Q)) 2 × (Q → P(Q)) 2 × (Q → Λ C )
• R S n|1 = R n|1 , R S n|2 = R n|2 • R S n|3 (α function) is close to R n|3 .
The last property about α defined at section VI-B is the only one which is modified:

• ∀q ∈ Q, α n (q) = q such that ∃q ∈ α n-1 , (q , σ c+n+1 , λ n,q , q ) ∈→ • property R S n|4 =R S n+1|3 is preserved • ∀q ∈ Q, R S 0|5 (q) = U nknown. • ∀q ∈ Q, R S n|6 (q) = λ (n,q ) (R S n|5 (q)). • ∀q ∈ Q, R S n|6 (q) = R S n+1|5 (q).
In figure 4b, if we take into account the statistical information, the sequence will be:

R 0 = ({A}, {A}, C 0 , C 0 ) R 1 = ({A}, {A, B}, C 0 , C 1 ) R 2 = ({A, B}, {A, B}, C 1 , C 2 ) R 0 = ({A, B, C}, {A, B}, α = id, {α(A) = ∅, α(B) = {A}, α(C) = {B}}, I 0 , I 1 ) R 1 = ({A, B}, {A, C}, {α(A) = ∅, α(B) = {A}, α(C) = {B}}, {α(A) = ∅, α(B) = {C}, α(C) = {A}}, I 1 , I 2 )
If transition of X at Y implies the statistical computation λ X,Y , then:

I 2 (A) = U nknown I 2 (B) = λ A,C (I 1 (B)) = λ A,C (λ B,A (U nknown)) I 2 (C) = λ B,A (I 1 (C)) = λ B,A (λ C,B (U nknown))}

D. The merging operation

After executing the Büchi automaton on each part of the trace, the result has to be generated using α function.

Let F, the set of states at element |σ|-1. Hence F = {q ∈ R |σ|-1|2 , such that ∃q ∈ R c|2 et q ∈ α n (q )}.

Theorem 1: The Verification of a property by sequential execution of a Büchi automaton is equivalent to verification of a property by parallel execution of a Büchi automaton.

If the trace is split into more than two sub-traces, then the results will have to be sequentially merged.

Proof: Function α is constructed such that it allows the selection of states of R |σ||1 coming from the set R 0 . The fusion operation consists in determining states of R 0 which have to be kept. The states to be kept are those of the set R c . As a reminder, we have R c ⊂ R 0 = Q. Hence, the selection of states of R |σ||1 coming from the set R 0 is equivalent to determining the generated states, for each state of R c . By definition of F, therefore, we have F = R |σ||1 .

Hence, the verification result can be computed using the end of trace algorithm defined in [START_REF] Ferlin | Vérification de propriétés temporelles sur des logiciels avioniques par analyse dynamique formelle[END_REF] and F.

E. Limitations

In [START_REF] Ferlin | Vérification de propriétés temporelles sur des logiciels avioniques par analyse dynamique formelle[END_REF], specific operators on variables have been defined in order to use a counter of variable modifications or the value of a variable at a given number of its modification before the current trace state. See the examples given in section III-A.

Because these elements are recursively computed, parallel verification cannot be used. A possibility to handle these operators could consist in adding new variables inside the trace which map to these specific asked values.

VII. EXPERIMENTS

In this section, the experiments results about the parallelverification implemented approach are presented. Tests are performed on generic traces, for the reasons given in section I. The goal of this section is to test the efficiency of the approach, and the to compare it, in term of execution time, with the classic approach.

A. Implementation

Due to industrial confidentiality, it is not possible to provide the source code of the program which implements the approach of this paper. Nevertheless, some technical elements are given below. The prototype is written in C++ and is 32, 900 length (number of lines of code), with 13, 000 length of header file.

1) A distributed architecture: The implementation of both approaches (sequential and parallel) follows a distributed architecture, in order to efficiently share the available hardware resources. Indeed, the reading action of the trace is separated from the Büchi automaton execution. Figure 5 illustrates the relation between the processes of the sequential approach. The execution process and the reading process communicate through a TCP/IP interface. The execution process can load information from one or several readers. Actually, trace reading can be split and shared out among several reader processes. Figure 6 illustrates the communication relations between the processes of the parallel approach. The fusion operation is not presented in this figure. Each execution process can interact with one or several reader processes. In figure 6 only the case where each execution process maps with a reader process is presented. 2) The trace format: For this prototype, the trace format is based on Xml. The trace is split into several blocks which contain the same number of states |b|. Hence, when the execution process needs a state σ i , it calls the reader process which returns the full block containing σ i . The first state σ i of this block contains the sub-trace σ i ...σ i+|b|-1 . The use of blocks has two explanations.

Firstly, it limits the network latency. Test have shown that verification time increases dramatically if states are sent one by one. Secondly, Büchi automaton split follows the block segmentation. Actually, traces are partial 4 . But elements σ i and σ i+|b|-1 are complete 5 for each block b beginning at σ i .

B. Load experiments

Two experiments are done: without statistical information, and with statistical information. The tests are done on the same generic execution trace for the two experiments. This trace has 10 million states and contains only one variable called x. For each trace state, the value of the variable changes. The domain of variation of the variable is -10; 10 . In order to do verification on different-sized traces, the verification can be stopped after a given number of states. Hence, we can verify a prefix of the given trace to simulate smaller traces.

Figure 8 gathers the results of the two experiments. The trace has been split into 2, 4, 8 or 10 sub-traces. The verification has been done for each sub-trace and the maximal time for each class of division has been placed in the graph. The memory process was launched in a 24Go-RAM computer, whereas the execution processes have been run on a 4Go-RAM computer with a Core-5i processor (2Ghz, 3Mo cache, 64-bits).

1) Without statistical information:

The targeted property to verify on this trace is (x 10 ∧ x -10) which means that constraint x 10 ∧ x -10 has to be true for each state of the trace. The goal of this test is to compare a sequential verification with a parallel verification of this property on the trace.

The Büchi automaton of this property, which is obtained by Ltl2ba, is defined by figure 7a.

According to figure 8a, comparison with the sequential verification, verification time is reduced by 40% when the trace is split into two sub-traces. If the trace is split into ten subtraces, verification time is reduced by 90%. The fusion time is weak: less than 0.01 seconds when the trace is divided into ten pieces.

2) With statistical information: The targeted property to verify is (♦ (x = 5)), which means that the constraint x = 5 occurs infinitely often. During the analysis of the trace, the number of states where the property x = 5 occurs is computed using the count counter whose initial value is 0.

Figure 8b shows that the verification of the property is more efficient when the trace is split than when the trace is sequentially browsed. Time verification is divided by two when the trace is split into two pieces,and divided by four when the trace is split into four pieces of trace... The fusion time is weak: less than 0.01 of a second when the trace is divided into ten pieces. This approach seems to be more efficient than the sequential one. The fusion operation which allows the merging of intermediate results requires less than 1 second. The next step of this work will consist in verifying traces which come from real industrial cases, with more complex properties, and traces with more variables. Actually, Büchi automata of complex A start

x 10 ∧ x -10

VIII. CONCLUSION

In this paper, we formalize the execution on a trace of a non-deterministic Büchi automaton, of a deterministic one and of a statistical one. This formalizations are based on the classic formal definition of a Büchi automaton, which is not complete to entirely describe each computation step. This formalization is the basis of the main contribution of this paper which consists in parallelizing the execution of a Büchi automaton on a trace.

The analysed trace is split into several sub-traces. Then, the Büchi automaton is executed on each sub-trace. The execution on the first sub-trace is classic, whereas the execution on each other sub-trace is different. Indeed, the Büchi automaton is initialized with the set of automaton states instead of the initial state. During the execution of the Büchi automaton on a sub-trace, a short-cut is computed between an automaton state at the beginning of the execution and its generated set of automaton states at the end of the execution. After that, this short-cut is used to perform the fusion of the results of the first sub-trace with the results of the second one, in order to deduce the global result.

This new approach has been used to carry out experiments on generic traces of different sizes. The verification times of the parallel approach are compared to the ones of sequential approach. The verification time is approximatively divided by the number of sub-traces analysed. As a result, this new approach allows the splitting of the verification of a property on an execution trace to limit time and memory requirements. The next step consists in verifying a property on a real industrial case with the parallel approach and comparing time requirements with the sequential one.

A theoretical improvement moderation of this parallel approach is that the specific operators on variables developed in [START_REF] Ferlin | Vérification de propriétés temporelles sur des logiciels avioniques par analyse dynamique formelle[END_REF] cannot be used because they imply an on-the-fly computation which is not parallelizable. However, we show that the verification time is better with the parallel approach than with the classic one. In addition, this approach includes computation of statistical information on traces. The relevance of the divide-and-conquer strategy using parallel execution had to be proved in this context. The computing performance is a validation.

A perspective could be to enrich the trace with necessary data to allow the verification of properties with the specific operators. Another direction of research could be to determine after how many splits we no longer save verification time. Finally, this new parallel method has to be confronted with a real industrial case.

Fig. 1 :

 1 Fig. 1: Execution of a deterministic Büchi automaton

Fig. 3 :

 3 Fig. 3: Execution of a non-deterministic Büchi automaton

Fig. 6 :

 6 Fig. 5: Sequential implementation Büchi automaton execution 1

Fig. 7 :

 7 Fig. 7: Büchi automata of the experiments

Fig. 8 :

 8 Fig. 8: Maximal time requirement for verification of a trace

Runtime Verification,

2001-2014, www.runtime-verification.org

Performing Enhanced Rail Formal Engineering Constraints Traceability, http://perfect.ifsttar.fr/Site/

Graphic Processing UnitActually, verification time increases linearly with the size of traces, and with the average number of variables which are modified at each state.

An element only contains modified variables

All variables with their current values are inside these elements.

ACKNOWLEDGEMENT

The work is funded by the French national research agency (ANR) in the context of the PERFECT Project. This project is also supported by the French I-TRANS competitive pole.


A Van Hamme

A El Habachi

W Samson

R Dumas

L Chèze

B Dohin

Gait parameters database for young children: The influences of age and walking speed

Keywords: healthy children, joint moments, joint power, regression analysis

Background. Reference databases are mandatory in orthopaedics because they enable the detection of gait abnormalities in patients. Such databases rarely include data on children under seven years of age. In young children, gait is principally influenced by age and walking speed. The influence of the age-speed interaction has not been well established. Therefore, the objective of the present study is to propose normative values for biomechanical gait parameters in children taking into account age, walking speed, and the age-speed interaction.

Methods. Gait analyses were performed on 106 healthy children over a large age range (between one and seven years of age) during gait trials at a self-selected speed. From these gait cycles, biomechanical parameters, such as the joint angles and joint power of the lower limbs, were computed. Specific peak values and the times of occurrence of each biomechanical gait parameter were identified. Linear regressions are proposed for studying the influence of age, walking speed and the age-speed interaction.

Findings. Most of the regressions achieved good accuracy in fitting the curve peaks and times of occurrence, and the normal reference targets of biomechanical parameters could be deduced from these regressions. The biomechanical gait parameters of a pathological case were plotted against the normal reference targets to illustrate the relevance of the proposed targeting method.

Interpretation. The normal reference targets for biomechanical gait parameters based on agespeed regressions in a large database might help clinicians detect gait abnormalities in children from one to seven years of age.

Introduction

During the first years of independent walking, considerable changes occur in joint kinematics and dynamics [START_REF] Chester | The identification of age-related differences in kinetic gait parameters using principal component analysis[END_REF][START_REF] Chester | A comparison of kinetic gait parameters for 3 -13 year olds[END_REF][START_REF] Dominici | Control of foot trajectory in walking toddlers: adaptation to load changes[END_REF][START_REF] Grimshaw | The 3-dimensional kinematics of the walking gait cycle of children aged between 10 and 24 months: cross sectional and repeated measures[END_REF][START_REF] Ivanenko | Kinematics in newly walking toddlers does not depend upon postural stability[END_REF]. Gait modifications have been studied to better understand gait maturation during the growth of children [START_REF] Samson | Biomechanical maturation of joint dynamics during early childhood : Updated conclusions[END_REF][START_REF] Sutherland | The development of mature gait[END_REF]. One of the difficulties in understanding gait maturation is the availability of an agematched reference databases for children, as suggested by Chester et al. [START_REF] Chester | An extended index to quantify normality of gait in children[END_REF].

The following reference databases for gait in children have been published: the temporal distance, kinematic and dynamic gait parameters of 10 toddlers aged 13.5 to 18.5 months old [START_REF] Hallemans | 3D joint dynamics of walking in toddlers A cross-sectional study spanning the first rapid development phase of walking[END_REF]; the ground reaction force patterns of more than 7000 children aged 1 to 13 years old [START_REF] Müller | Static and dynamic foot characteristics in children aged 1-13 years: a cross-sectional study[END_REF]; and the kinematic and dynamic parameters of 20 Chinese children aged 7 to 12 years old (Bacon-Shone and [START_REF] Bacon-Shone | Gait of normal Hong Kong Chinese Children: the bootstrap approach[END_REF]. These studies demonstrated the influence of age on biomechanical gait parameters. In mid childhood, "sagittal joint kinematics, moments and powers are predominantly characterized by speed of progression, not age", as reported by [START_REF] Stansfield | Sagittal Joint Kinematics, Moments, and Powers are predominantly characterized by speed of progression, not age, in normal children[END_REF]. The major relevance of the latter study is that dimensionless walking speed should be preferentially considered rather than age to compare healthy and pathological gaits in children. These conclusions are based on children aged 7 to 12 years old and could be different for younger children. Moreover, Schwartz et al. (Schwartz et al., 2008) described the gait of 83 typically developing children walking at a wide range of speeds and displayed spatio-temporal, kinematic, kinetic and electromyographic data for children between 4 and 17 years old. In this study, the influence of speed variation on this population was obvious from the graphs presented, but the age influence was ignored. [START_REF] Stansfield | Regression analysis of gait parameters with speed in normal children walking at self-selected speeds[END_REF] proposed a regression analysis of biomechanical gait parameters as a function of walking speed; however, the results were obtained with relatively low determination coefficients (R²)

A C C E P T E D M A N U S C R I P T (i.e., less than 0.3 except for the temporal distance parameters and components of the ground reaction force (GRF)). The study included 16 children aged 7 to 12 years old.

Clinical indices based on kinematic data (the Gillette Gait Index [START_REF] Schutte | An index for quantifying deviations from normal gait[END_REF], the Gait Deviation Index (Schwartz and Rozumalski, 2008) and the Gait Profile Score [START_REF] Baker | The gait profile score and movement analysis profile[END_REF]) and on dynamic data [START_REF] Rozumalski | The GDI-Kinetic: a new index for quantifying kinetic deviations from normal gait[END_REF] have been proposed. Age and speed variations were not considered in their calculations, and there are difficulties in comparing healthy and pathological gaits in children because of walking speed differences.

The objective of the present study was to propose normative values for biomechanical gait parameters taking into account age, walking speed, and the age-speed interaction. This is not easily achievable with group corridors (i.e., mean +/-standard deviation) unless a very large number of age-speed groups are considered, but can be more simply achieved by establishing normal reference targets based on regression models. Therefore, this study establishes a large database of more than 100 young children (from one to seven years old). Because the objective of the paper is to measure the influence of age and walking speed, the database was collected to provide a wide range of ages and speeds (the children walked at a self-selected speed). The influence of these factors was analysed using regression models that link age, walking speed, and their interaction on biomechanical gait parameters (kinematic and dynamic data). These regressions are important for studying the influence of the tested factors; however, they are somewhat impractical for clinical application. Normal reference targets were constructed based on regression models that allowed the pathological biomechanical gait parameters of children to be plotted against the normative values, taking into account age, walking speed, and the age-speed interaction. The relevance of the method is illustrated with one pathological case.

A C C E P T E D M A N U S C R I P T

Materials and Methods

Population and experimental set-up

Gait analysis was performed on 106 healthy children (from one to seven years old). The participant characteristics are presented in Table 1. One child could be measured several times during its growth. All of the children were independent walkers from the first examination, and clinical examination did not reveal any orthopaedic or neurological disorders. The local ethics committee approved the study. The children were included in the study after clinical examination and when their parents consented to involvement after having been informed about the protocol.

Twenty-four skin markers were fixed on anatomical landmarks of the pelvis (the right and left anterior and posterior superior iliac spines) and the lower limbs (the great trochanter, medial and lateral epicondyles, anterior tibial tuberosity, medial and lateral malleoli, calcaneus, first and fifth metatarsal heads and hallux) [START_REF] Samson | Biomechanical maturation of joint dynamics during early childhood : Updated conclusions[END_REF][START_REF] Samson | 3D joint dynamics analysis of healthy children ' s gait[END_REF].

The children walked barefoot at a self-selected speed. Fifteen to twenty gait trials were measured for each subject using a Motion Analysis  system with eight Eagle  cameras (Motion Analysis Corporation, Santa Rosa, California, USA) at 100Hz and three Bertec  force platforms (Bertec Corporation, Columbus, Ohio, USA) at 1000Hz. Only trials with valid dynamic data were selected (i.e., one foot and only one foot on one forceplate), providing between one to six gait trials per gait analysis. In total, 1253 gait cycles were retained.

Data processing

After filtering (low-pass zero-lag, 4th-order, Butterworth filter, with a 6-Hz cut-off frequency), the marker trajectories were obtained in an Inertial Coordinate System (ICS) [START_REF] Wu | ISB recommandations for standardization in the reporting of kinematic data[END_REF]. The hip joint centre localisation was determined using the regression models established by Harrington et al., selecting only the data from healthy children [START_REF] Harrington | Prediction of the hip joint centre in adults , children , and patients with cerebral palsy based on magnetic resonance imaging[END_REF]. The inertial parameters were determined using the regressions established by Jensen [START_REF] Jensen | Changes in segment inertia proportions between 4 and 20 years[END_REF]. The three orthogonal axes (X, Y, and Z) corresponding to each segment coordinate system (SCS) were built following the International Society of Biomechanics recommendations [START_REF] Wu | ISB recommendation on definitions of joint coordinate system of various joints for the reporting of human joint motion--part I: ankle, hip, and spine[END_REF]. The quaternion was extracted from the attitude of these axes in the ICS. The angular velocities of the proximal and distal segments were obtained in the ICS using quaternion algebra and were subtracted to compute the (relative) joint angular velocity, . The net 3D joint moments, M, were computed in the ICS by bottom-up inverse dynamics [START_REF] Dumas | A 3D Generic Inverse Dynamic Method using Wrench Notation and Quaternion Algebra[END_REF], with the force platform's data resampled at 100Hz. The power, P, was computed in 3D by the dot product between M and .

A C C E P T E D M A N U S C R I P T (
The joint moments, M, were expressed in the joint coordinate systems [START_REF] Desroches | Expression of joint moment in the joint coordinate system[END_REF], and M and P were re-sampled on a percentage of the gait cycle and were expressed using the dimensionless scaling strategy [START_REF] Hof | Scaling gait data to body size[END_REF], with the leg length (the distance from the ground to the great trochanter) used as a metric value. The walking speed was defined from the initial contact of one foot to the next initial contact of the same foot (one gait cycle)

and was expressed with a dimensionless parameter [START_REF] Hof | Scaling gait data to body size[END_REF]. The moments are in units of N.m/ , the powers are in units of , the GRF is in units of and the walking speed is in units of (with m 0 indicating the body mass, l 0, the leg length and g, the acceleration of gravity).

The gait trials were not averaged per subject. Data from both right and left strides were included, taking into account the sign conventions. The peak values and the corresponding times of occurrence were identified on the curves displaying kinematic and dynamic parameters (Table 2). The times of occurrence were expressed as a percentage of the gait cycle. The coefficients of the linear regression models were estimated considering age, walking speed, and the age-speed interaction as the model inputs as well as the peak values

A C C E P T E D M A N U S C R I P T
and times of occurrence as the model outputs. We computed the confidence intervals of the regression models to define the normal reference targets.

Data analysis

Regression models

To establish a link between the biomechanical gait parameters (the peak values and times of occurrence) and age, walking speed and their interaction, the regression model used was as follows:

where Y was the estimated output variable, age and speed were the input variables and a, b, c, and d were the regression coefficients. These coefficients were estimated using the least squares method. The determination coefficients (R²) and t-test p-values were calculated to evaluate the goodness and relevance of the fit, respectively. The regressions of body mass, m 0 , and leg length, l 0 , with age are also provided to allow for comparison with previous studies that did not use a dimensionless parameter (Table 3).

Normal reference targets

The aim of the study was to propose normal reference targets for clinical use. These reference targets were achieved on the regression residual (i.e., using the difference between the measured output value and its estimation by the regression). The calculation of the standard deviation of the residuals allows for the estimation of the confidence interval of the output variable. The confidence interval was calculated for the peak value and time of occurrence, defining an ellipse of confidence regarding the estimated output variable. For each estimated output variable, the confidence interval at 95% was computed using the following formula: for the superimposing of the normal reference targets (i.e., the ellipses of confidence centred on each biomechanical parameter peak) on the patient gait curves. As an illustration, a pathological gait was evaluated using the normal reference targets. The patient had cerebral palsy with right hemiplegia. He was seven years old and walked at 0.39 .

Results

Regression models

Most of the regressions were significant (the p-values of the t-test were less than 0.05). R² was greater for the peak values than for the times of occurrence. R² values higher than 0.4 were obtained on some peaks for the values and times of occurrence, especially for the knee and hip dynamic parameters. The standard deviation values, allowing for the computation of the normal reference targets, are provided in Table 3. All of the results of the regression analysis are available in the supplementary materials. The results of a global sensitivity analysis [START_REF] Plischke | An effective algorithm for computing global sensitivity indices (EASI)[END_REF] are also presented. This analysis was performed on age, walking speed and their interaction and provides information about the contribution of each parameter in biomechanical gait parameters.

Normal reference targets, an application example

Figure 1 shows an example of the application of the normal reference targets to a pathological gait analysis (cerebral palsy with right hemiplegia). The peak values and times of occurrence were calculated for the ankle, knee and hip power on the lower limbs using the regression model (inputs: age 7 and walking speed 0.39 ). Normal reference targets were built for each estimated value, taking into account the confidence interval of the peak values and the times of occurrence. For the ankle, the second peak target was not represented because the regression analysis was not significant (see supplementary materials). The gap between the normal and pathological data was gauged with the proposed targeting method.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

10 During the stance phase, the left ankle power is close to zero, whereas the right ankle absorbs negative power more than in healthy children. For the knee and hip, the right hemiplegic limb peak powers are lower than the means of the normal data, and the peak powers on the left lower limb appear to compensate for this phenomenon (even if most of the peaks still fall within normal limits). The child with a pathological gait developed a larger negative or positive power than that of the healthy children, depending on the side, for the knee and hip, especially during the swing phase.

Discussion

This study presents a large biomechanical gait output database of healthy young children (younger than seven years old) walking at a self-selected speed. A regression analysis was performed to estimate the normal reference targets for each peak of the biomechanical gait parameters for the healthy children. An application of the normal reference targets in an illustrative pathological case was proposed. More than showing the influence of age and speed, the objective of this study is to provide targets based on the age-speed regressions to compare pathological cases (of given ages and walking at typically lower speeds) with a reference.

Only walking at a self-selected speed could be analysed with this population (i.e., speed could not be imposed on very young children). Yet, a large range of speeds was achieved ([0.1-0.7 ]), which was a similar result to that of a previous study on older children (Schwartz et al., 2008), where very slow to very fast conditions were imposed.

Although the proposed linear regression model was still relatively simple (taking into account age and speed), it provided a better approximation of peak values than only considering walking speed [START_REF] Stansfield | Regression analysis of gait parameters with speed in normal children walking at self-selected speeds[END_REF]. More complex models could be explored in future studies to fit the experimental data better (e.g., including the exponential or logarithmic

A C C E P T E D M A N U S C R I P T functions)

. This simple linear model provided significant regressions for most of the studied biomechanical gait parameters. The difference between R² for the peak values and for the times of occurrence could be explained by the fact that the variability of the times of occurrence is higher than the variability of the peak values. The best correlations were found for the dynamic values, perhaps because of the better repeatability of the dynamic rather than the kinematic parameters [START_REF] Steinwender | Intrasubject repeatability of gait analysis data in normal and spastic children[END_REF]. The best R² values were obtained for the power values. The calculation of power, taking into account the joint angular velocity, was obviously linked to the walking speed and could explain the greater R² values. Stansfield et al. [START_REF] Stansfield | Regression analysis of gait parameters with speed in normal children walking at self-selected speeds[END_REF] did not find a significant regression for the second peak of the vertical GRF, assuming that this peak could be linked with the body's control of stability instead of with the maintenance of speed. Including age in the model, the regression on the second peak of the vertical GRF was significant and had an acceptable R² value in comparison with that of Stansfield's results (i.e., R² = 0.18 vs. no significance, respectively). These differences could be explained by our younger population.

The standard deviations were large because of the high variability of gait in young children.

The calculation of the confidence intervals provided normal reference targets, allowing for the evaluation of any gait cycle and defining the normative values for comparison with pathological cases. These normative values take into account age, walking speed, and the agespeed interaction, which was not possible with simple group corridors (i.e., mean +/-standard deviation). By including all of the gait cycles, the complications linked with group definitions (based on the age and/or speed) and, especially, a definition of the boundaries for groups were avoided. The visual comparison with an illustrative pathological case illustrates the clinical potential of the targeting method. The patient was seven years old, which corresponded to the upper limit of our database of healthy children. Application to younger patients and to a larger number of patients with different pathologies is necessary to confirm the suitability of the

A C C E P T E D M A N U S C R I P T

normal reference targets. These targets, based on the age-speed regressions established for a large database, could be a complement to existing clinical indices (e.g., the Gillette Gait Index [START_REF] Schutte | An index for quantifying deviations from normal gait[END_REF] and the Gait Deviation Index (Schwartz and Rozumalski, 2008)).

Conclusion

This study presented a large biomechanical gait parameters database of young healthy children (including more than 100 children) and proposed an original regression of these parameters with age, walking speed, and the age-speed interaction. The regressions were calculated for the peak values of the biomechanical gait parameters and their times of occurrence. A method was proposed to define normal reference targets that might help clinicians detect gait abnormalities in children from one to seven years of age. [2-3[ [3-4[ [4-5[ [5-6[ [6-7]

A C C E P T E D M A N U S C R I P T

  Figure: 1

  [Y-1.96*SD(Y-Y mes ); Y+1.96*SD(Y-Y mes )], where Y mes is the measured output variable (the peak value and time of occurrence). The knowledge of these standard deviation values allows A C C E P T E D M A N U S C R I P T

Figure 1 :

 1 Figure 1: Ankle, knee and hip powers (in ) of a hemiplegic child for the left

Table captions Table 1 :

 captions1 Participant characteristics according to age. For each group, age boundaries (e.g., [1-2[ for the 1st group) mean that all children between their 1 st birthday plus one day and 2 birthday were included. An exception was accepted for the last group (i.e.,[6][7]) in which one child who had already had his 7 th birthday was included.

Table 2 :

 2 Description of peak identification and corresponding abbreviation

Table 3 :

 3 Regression models for mass, m 0 ; leg length, l 0 ; and the biomechanical gait parameters with R²> 0.1 for the peak value and the time of occurrence: a, b, c, and d: Regression coefficients, R²: determination coefficient, p: p-value (****: p<10 -5 ), SD: standard deviation of errors. N/A: Not Applicable. Abbreviations (e.g., A_A2) are detailed in Table1.

	Table 1:
	Age Group (years old)	[1-2[

Table 2 :

 2

	Total

Table 3 :

 3

Research highlights

 This study included more than 100 healthy children between one to seven years old.

 Our study considered both age and speed of progression influence on children gait.

 Regression models were processed for the peak values and their times of occurrence.

 An original method of normal reference targets was proposed.

 The normal reference targets might help clinicians in detecting gait abnormalities