Keywords: Pantograph-catenary interaction, Variability, Dynamic simulation

Currently, pantograph-catenary dynamic simulations are mainly based on deterministic approaches. However, the contact force between catenary and pantograph depends on many key parameters that are not always quantified precisely and can vary in time and space. To get a better chance of addressing extreme or combined critical conditions, methodologies to consider variability are thus necessary. Aerodynamic forces and geometrical irregularities of catenaries are thought to be significant sources of variability in measurement and this paper proposes methods to take them into account. Results are compared with measurements to correlate the effect of the considered parameters with experimentally observed variability. Finally, a virtual certification example is shown, with a study of the influence of speed on the impact of variability.

Introduction of variability into pantograph-catenary dynamic simulations

Olivier Vo Van * ‡ † , Jean-Pierre Massat ‡ , Christophe Laurent π and Etienne Balmès † † PIMM, Arts & Metiers ParisTech, Paris, France; π Vibrateam (Vibratec Group), Lyon, France; ‡ SNCF research department, Paris, France

1 Introduction
Today, the minimization of rolling stock impact on infrastructure is a leading objective in the railway sector. In this context, requirements on current collection quality given by the Technical Specifications for Interoperability (TSI) are likely to become much restrictive in the next few years. However, performance of existing pantographs has almost reached their physical limits, especially in the case of train sets with multiple pantographs.

Since inline test give limited knowledge of the system response and are limited due to conflicts with commercial traffic, virtual tools [START_REF] Bruni | The pantograph-catenary interaction benchmark[END_REF] are certain to represent a growing part of the effort on optimisation and sensitivity studies. Examples of optimisation studies are [START_REF] Gomez | A model to obtain optimal designs of railway overhead knuckle junction using simulation[END_REF], which aims to limit breakages and to improve current collection quality under knuckle junctions, or [START_REF] Pombo | Influence of pantograph suspension characteristics on the contact quality with the catenary for high speed trains[END_REF] in which the parameters studied are pantograph suspension characteristics. For sensitivity analyses, Ref. [START_REF] Kim | State sensitivity analysis of the pantograph system for a highspeed rail vehicle considering span length and static uplift force[END_REF] analyses span length, static uplift force and pantograph system parameters, and Ref. [START_REF] Aboshi | Installation Guidelines for Shinkansen High Speed Overhead Contact Lines[END_REF] analyses catenary geometry through installation error indices.

Based on SDTools libraries, SNCF has developed an efficient and powerful finite element library for the simulation of pantograph/catenary interaction named OSCAR [START_REF] Balmès | OSCAR 1.0 user's guide[END_REF]. This software is certified against EN50318 and validated against measurements. It is routinely used for design studies. Although OSCAR has been developed on deterministic approach, it is quite clear that a basic statistical strategy would give a higher confidence in results and be useful for industrial use. This paper presents the first developments of this new methodology, introducing three kinds of irregularities, namely catenary geometry irregularity, catenary wear and transient aerodynamic forces.

The case study is based on the French pantograph CX made by Faiveley running under the east high-speed line between Paris and Strasbourg allowing a maximum speed of 320 km/h. Some homologation measurements done by SNCF's testing agency (Agence d'Essais Ferroviaires) and others resulting from maintenance measurement campaigns are used to establish the source of variable parameters impacting performance.

Geometrical irregularities of catenaries, contact wire wear and aerodynamic forces are first introduced as physical sources of variability and characterized by measurements. Section 3 then proposes methodologies to take these sources into account. Results of simulations are then compared to each other and with measurements to identify the main physical phenomena. In section 4, the proper selection of observation criteria is shown to be a critical aspect of the proposed analysis. Finally, section 5 details an example typical for virtual certification highlighting the importance of vehicle speed.

2

Sources of variability analysis

Three physical perturbations and associated measurements are studied in the following section:

 geometric irregularities, which correspond to the actual height of the contact wire, given in [cm],  contact wire wear which is defined by the contact wire section, given in [mm 2 ],  aerodynamic irregularities, defined as a perturbation around the mean aerodynamic applied force, given in [N].

Geometric irregularity

Geometric irregularity characterizes differences between nominal and real geometry of the catenary particularly at mast and dropper positions. The nominal geometry is defined by perfectly regular sag of 1/2000 and a fixed height identical for every mast. This study focuses on parameters like dropper length and mast height, which affect the contact wire static height. Figure 1 shows that the measurement is much less regular than the theoretical design. The spatial and spectral representations are given for measurements in solid line and for the nominal height of the catenary in dashed line.

Measurements are low-pass filtered using a 20-points moving average method. This corresponds to a filter at 0.25 m -1 or 4m, which is smaller than the smallest interdropper distance. Higher frequency irregularities are not taken into account since they cannot be caused by geometric irregularities. In-line measurements are used to determine a mean wear over the contact wire. These global values permit to define coefficients of a wear law defined experimentally [START_REF] Bruni | Pantograph-catenary dynamic interaction in the medium-high frequency range[END_REF] as:

= 1 × 1 + 0 -× 0 × 0 × × 2 × ( ) × 2 • (2.1) with ( ) = 1 2 (2.2)
where A is the surface removed by one pantograph passage, F c the load, i the current and V the train speed. V 0 , i 0 and F 0 are used for normalization and taken equal to nominal values, i.e. V 0 =3m/s, i 0 =426A, F 0 =190N. k 1 , k 2 ,  and  are the adimensional parameters defined to fit with measurements. H is the hardness of the softer of the two materials, H=110HR. Finally, k r is a coefficient used to interpolate the relationship between R and F c , k r =4.97 10 -2 Ω.N 1/2 .

(a) (b)

Aerodynamic irregularities

The aerodynamic irregularity is defined as a time varying vertical force fluctuating around the mean aerodynamic force applied on the pantograph. This variation is due to changes in different parameters such as train speed or wind speed and direction. Available measurements are used to compute the mean aerodynamic force. They are filtered at 1 Hz and give a reasonable indication of the variation range for the aerodynamic force. For the considered case at 320km/h shown on Figure 2, one has a mean of 96N and a range of variation between 83N and 109N with a 95% confidence interval. The use of the same data for irregularity studies will be discussed in the next section.

Figure 2. Measured aerodynamic forces (normal lines) and 95% probability interval (bold lines)

Integrating uncertainty in simulation tools

The method to account for each uncertainty depends on how each parameter influences the dynamic behaviour and will be detailed in subsections. All considered cases however share a random signal generation process. The classical method used here is to affect a random phase  n , in the [-] interval, to a spectrum of known amplitude A n and to compute the response by inverse Fourier transform.

where zero mean is achieved by setting A 0 =0 and the spectrum is constrained to be symmetric (that is A N-n =A n and N-n=-n) to obtain a real realization by inverse Fourier transform.

Geometric irregularities

Between droppers and spans the geometry is only linked to tension, beam section and gravity which are well known. The proposed strategy is thus to use random values for contact wire height at dropper and span locations. Wear induces height modifications whose dynamic impact is very low for long wavelengths. One can thus neglect geometry modifications induced by wear.

( ∆ ) = 2 + =1 (3.1)
One thus first generates a random signal with equation (3.1) using the amplitude spectrum shown in dashed line in Figure 1(b), which corresponds to the filtered measurement. Given this signal, an inverse static computation is performed by modifying dropper length of the nominal catenary in order to converge to the target position of contact wire at dropper and span locations. The resulting catenary mesh is then used in dynamic simulations, taking into account geometric irregularities. Figure 3 shows an example of sag obtained from the nominal catenary model and the target height after a static computation. Figure 3b illustrates the difference in sag and demonstrates that the target geometry is reached with a gap of less than 0.1mm along the contact wire. The highest variations come from the approximation of the masts height, which is negligible.

Contact wire wear

The spatial distribution of wear is not really random. Contact force is indeed usually higher in the same area of the contact wire. One thus prefers to generate a wear profile by using equation (2.1). One applies the contact force computed with the nominal geometry as F c . The wear profile A(x), obtained for one cycle, is then multiplied by a target number of cycles before renewing the contact force computation. The full lifespan of the contact wire is thus separated into n linearly spaced simulations, n representing the number of cycles varying between 0 and 10 8 .

As shown in Figure 4, wear induces a perturbation dz to the vertical position of the contact point. In the time integration, solved with Newmark implicit method, one thus considers

K • q + C • q + M • q = k c z cat + dz -z panto (3.2)
where K is the stiffness matrix, C the damping matrix, M is the mass matrix, z cat , height

(b) (a)
of the contact wire corrected by the section change dz and z panto , height of the pantograph. q is the displacement vector and k c the contact stiffness. The second effect is a section reduction of the contact wire which leads to a decrease of the local stiffness. As OSCAR uses FEM, stiffness variations cannot have higher frequency than the one induced by the element size. A constant stiffness is applied for each beam element of the contact wire, by averaging the random section wear signal over the element length. Both phenomena are taken into account at the same time in this study, which does not aim to analyse their impact separately.

Aerodynamic irregularities

As seen previously, measured aerodynamic forces are filtered at very low frequency. While defining a probability distribution for the aerodynamic force would be a classical approach, it does not contain any spectral information so that a different approach was preferred here. As shown in Figure 5, one defines a design force spectrum by extending the measurements to higher frequencies and uses this spectrum to generate a vertical load on the upper mass of the pantograph model that is random in time. The extension was done by roughly taking the mean value of the spectrum amplitude computed at 1Hz as fixed amplitude between 1Hz and 100Hz. The frequency under 1Hz corresponds to a spatial frequency under 0.012m -1 for a nominal speed of 83m/s, which is under the span frequency, equal to 0.0185m -1 (=1/54m). The value 100Hz has been taken arbitrary. understanding of the impact of irregularities in particular due to the presence of phase differences with the measurements. To address this issue different comparison tools have been considered. Of the four detailed here, two are used to compare signals together, and two are used to compare signals to measurements.

Criteria selection

Statistical values

The main goal of pantograph-catenary simulation is to evaluate the current collection quality. The usual approach is to divide the contact-force standard deviation by mean contact force, which has to be as low as possible and under 0.3 to meet standards. This is a statistical criterion which, assuming the load to follow a normal distribution, corresponds to a probability to lose contact of less than 1%. This criterion is called coefficient of variation in statistics. For n simulations, the first statistical moments can be considered as discrete variables of i, i∈ [1, n] on which standard deviation can be computed: where and E denote the standard deviation and expectation operators computed over the n simulations. S3 and M3 are the main criteria to define

 S 1 = (F m,i ), with F m,i the mean contact force of the i th simulation,  S 2 = (F c,i )), with F c,i the time signal of contact force of the i th simulation,  S 3 = ((F c,i )/F m,i ),  Min 3 =Min((F c,i )/F m,i ),  Max 3 =Max((F c,i )/F m,i ),  M 3 = E((F c,i )/F m,i ),  M nom = (F c, nominal )/F m,

Spectrum of differences

In order to compare simulations together more precisely than using a simple scalar criterion, the computation of the difference between two contact force spectra gives information about which frequencies are more affected by modifications. All simulations are compared with the nominal model. For each type of irregularity, results are averaged over a hundred simulations.

Frequency band analysis

Frequency Band Analysis (FBA), see more details in [START_REF] Massat | Simulation tools for virtual homologation of pantographs[END_REF], is a comparison tool used to get a sharper discrimination of unsynchronized signals than a simple comparison of Power Spectral Density (PSD) without observing shift of frequency but only mean amplitude in frequency bands. FBA gives a graphical representation of variance against frequency. Defining a range of variation over a set of simulations gives information about the amplitude of variation of energy in a frequency band, which can easily be compared with measurements. In our case, frequency bands will be centred on the harmonics of median span frequency, which appeared to be the main contributors in the PSD.

Frequency coherence function:

The aim of the coherence function [START_REF] Massat | Simulation tools for virtual homologation of pantographs[END_REF],

, 2 ( ) = ( ) 2 ( ) • ( ) (4.1)
is to assess and quantify if a linear relationship exists between two sets of frequency spectrums. As the numerator is the magnitude squared of a cross power spectrum and the denominator a multiplication of two power spectrum (real valued), the coherence function is real and has values between 0 (no linear relationship between the two signals) and 1 (perfect linear relationship between the two signals). The result is usually given in percentage.

In practice, each signal is decomposed in I equally spaced samples centred on a post. Inter and auto spectra are then computed as the mean over I samples

= 1 • ; • ; * =1 (4.2)
with

; = ; • -• 2 • • -1 =0 (4.3)
where  X is the discrete Fourier transform of the time signal x,  F E is the sample frequency,  f k is the frequency of evaluation of the coherence function,   x,y 2 (f k ) is the coherence function (without units) between x and y at frequency f k ,  S xy (f k ) is the Cross Power Spectrum between the signals x and y at frequency f k ,  S xx (f k ) Power Spectrum of the signal at frequency f k , Since the model and test are random, one has l simulation and m measurements. The final coherence at a frequency is the average l x m coherences. The curve will be shown in black. The process is repeated for each measurement with respect to all others and the average curve is shown in red. It is assumed that simulation coherence with measurements cannot be higher than measurements with each other. To be able to compare with nominal results (i.e. with no irregularity), its coherence with measurements is shown as a green line.

To ease the comparison of various cases a scalar measure of coherence is introduced by computing the surface under the curve defined by  x,y 2 (f k ) for frequencies between 0 and 20Hz. This scalare represents the percentage of coherence between the group of simulations and the group of measurements for this range of frequencies.

Sample results

Simulations are performed with OSCAR for a nominal model containing 3600 elements and 3200 nodes, for a section containing 14 spans. Train speed is fixed at 320km/h with a simulation over 700m (10 spans). Time consumption is about 2 minutes, which gives the possibility to make 3000 simulations for all irregularity cases in around 100 hours. The work is distributed over several computers.

Comparison of simulations

The most widely used criterion is /F m . Its mean M 3 must be as small as possible for a good collection. M 3 -M nom gives the mean impact of irregularities. The range of variation defined by Min 3 and Max 3 gives additional information of the tendency. Associated values are given in Table 1. As a complement Figure 6 gives a histogram of /F m . Geometric irregularity has a strong impact on current collection quality. This impact can be strongly positive or negative as the values Max 3 and Min 3 show. The distribution is almost centered on the nominal value and the standard variation induced by this irregularity on /F m is ±3.7% of the mean value (S 3 /M 3 ).

Wear irregularity is the second most influent irregularity at 320km/h according to Figure 6. But the way it is defined here only induces positive impact on current collection. The contact force is globally smoothed where wear is high, which is located on the area of high force. Consequently, standard deviation of contact force (S 2 ) is reduced.

Aerodynamic irregularity has a very small impact on standard deviation of the contact force (S 2 ). The main modification is induced on mean force (S 1 ), which can be explained by the way the simulation is performed. The aerodynamic signal is indeed averaged to zero over the whole simulation, which is different from the period on which statistical values have been computed. The impact of aerodynamic irregularity can therefore be taken into account simply by varying the mean force with ±13N as seen previously, which leads to a standard variation on /F m of ±9% of M 3 . Figure 6. Distribution of the geometry (black), aerodynamic (grey) and wear (white) irregularities for 1000 simulations each

The second result, shown in Figure 7, is the spectrum of contact force differences for each irregularity. Span frequency and harmonics are given with vertical dotted lines. It appears that the impact of geometric irregularity is focused on span frequencies, which is consistent with the fact that the energy of contact force is higher at these frequencies.

The same behaviour is observed with wear irregularity in a more ambiguous way. It is although not the case of aerodynamic irregularities. This observation confirms the previous conclusion about the lower impact of aerodynamic irregularity on contact force variation around its average. Figure 7. Mean of spectra of contact force differences for geometry (solid line), aerodynamic (dashed line) and wear (dotted line) irregularities, with span harmonic frequencies (vertical lines)

Measurement-simulation comparison

In order to compare model and test variations, the FBA and coherence criteria give visual support.

On the FBA graph in the left column of Figure 8, surfaces with transparency represent the range between maximum and minimum observed values for corresponding set of simulations (in dark grey) and measurements (in light grey). Nominal behaviour is shown in black dotted line. One can see the energy variation of the contact force at a specific frequency, namely, the range of variance at this frequency. Geometric irregularity represents far better the amplitude of variations observed in measurements than the two other which are nearly unnoticeable. For frequencies higher than 0.1m -1 , wear irregularity seems to be more representative than aerodynamic irregularity. The low level of FBA variability in the aerodynamic sensitivity analysis is coherent with the low level of spectrum differences in Figure 7.

On frequency coherence charts in the right column of Figure 8, frequency dependent spectrum linearity is assessed between the set of measurements and the set of simulations. Coherence between simulations and measurement can at best be equal to coherence between measurements themselves. Wear and aerodynamic irregularities clearly improve global coherence, which is represented by the grey area. This means that phase and frequency variations induced by these irregularities are representative of measured variations unlike what happens with geometry irregularity.

Values of the scalar coherence criterion, given in Table 2, give a similar trend than the qualitative conclusion on frequency dependence. In particular, the mean coherence between geometry-irregularity simulations is quite low. One thus expects that geometry irregularity induces significant non-linearities into the pantograph-catenary system. Understanding why FBA and coherence criteria lead to different conclusions on the representativity of variability models seems an interesting direction for future work.

Example of use in virtual certification

Virtual certification is a promising approach to use numerical simulations. The European project PantoTRAIN (see [START_REF] Unife | PantoTRAIN WP3 deliverable 3[END_REF]) led to a full procedure to extend pantograph certification to different catenary systems. The recent certification of a new pantograph in France and Spain was an opportunity to lead a blind test in parallel to inline tests. This pantograph has the particularity to be insensitive to the airflow and has a twice lower bow mass. Wear is very low on the catenary used for the inline tests. Geometry irregularity is thus the main known source of variability.

The analysis is based on three sets of dynamic contact forces:

 Inline tests, which are dynamic forces obtained directly from an equipped pantograph on the train at four different speeds,  Simulations with nominal geometry, which correspond to simulations made at the same speeds and mean forces as the inline tests with the nominal catenary defined in section 2.1  Simulations with measured geometry, which correspond to simulations made at the same speeds and mean forces as the inline tests with the geometry measured under the catenary modelled Figure 9 displays the simulation value of (F m -3 with black triangles in comparison to the measured values in light grey triangles for the same values of mean forces as shown by the circles. The tests show a dispersion of mean force which generates dispersion on (F m -3This criterion is almost the same as /F m , since it follows the same goal, but has to be over zero and consequently as high as possible. The curves F m,min and F m,max correspond to the minimal and maximal values of the mean contact forces authorized by the recent TSI.

Table 3. Deviation on  depending on the speed and presence of geometry irregularity shows three comparisons between dynamic inline tests, simulations with nominal geometry and simulations with measured geometry over the section used for the model.

The difference on the prediction of /F m between tests and simulations is computed for the four speeds. At each speed, the value of /F m is averaged. The left column is obtained by comparing inline tests with the nominal catenary and the central one with the measured height of the catenary. The right one compares simulations obtained with measured and nominal geometries. The difference between both nominal and measured geometries is only of 2% at 320km/h. At every other speed, the difference between simulations is between 14% and 17%, which means that either the speed or the mean force (or both) has an influence on the geometric irregularity impact on contact force.

Inline test are based on different sections and each of them has a geometry that differs from that used in simulations. In the central column, the value of /F m of inline tests however seems to be always lower than simulations. This implies that the representation of geometric irregularity used here includes cases that have significantly different impact from what is actually found in test. One thus includes too many geometry configurations within the way the geometry irregularity is defined.

Introduction of speed and mean force variations

Since the impact of variable geometries was shown to depend on speed, one seeks to further characterize this dependence. In order to compare simulation results to those obtained for the virtual certification, one computes two sets of simulations at 300km/h and 320km/h with the mean force F m,300 obtained at 300km/h and the thousand geometries previously generated. Figure 10(a) shows the distribution of mean force for variable geometries at three different design points. The mean force distribution show little impact from speed, or the mean force. For the average of /F m and , as show Figure 10(b) (c), speed has clearly a higher influence than mean force. It is difficult to draw more conclusions from this figure.

Figure 11 shows the histograms of /F m for two speeds and a correlation plot where each point is associated with a geometry and the resulting /F m at the two speeds.

The correlation plot is elongated along the diagonal, which implies that there is a strong correlation for both variation of mean force and speed. One can observe a higher dependence between two different mean forces F m than between two different speeds. This qualitative observation can be quantified by the correlation defined as

1,2 = 1 , 2 1 × 2 (6.1)
which gives r Fm =96,4% and r Speed =88.9% where r Fm is the correlation between /F m of the sets (320km/h, F m,320 ) and (320km/h, F m,300 ) and r Speed is the correlation between /F m of the sets (300km/h, F m,300 ) and (320km/h, F m,300 ). Speed thus has more influence than mean force when considering the impact of geometry on contact force. One now seeks to evaluate how representative the random geometry is of the measurement by comparing the results of virtual certification with those of this method. By getting in the same conditions of speed and force as those which led to the results of Table 3, one can define a range of variation of /F m at (300km/h, F m,300 ), which can be deduced knowing the variation of /F m at (320km/h, F m,320 ) from the nominal. Figure 12. Dependence on (speed, mean force) of the criterion /F m between two design points, (300 km/h, F m,300 ) against (320 km/h, F m,320 ). shows the combined influence of speed and mean force. The nominal values of /F m for both design points are indicated in thin solid line. The intersection is thus the position of the nominal geometry. In order to locate the position of the geometry irregularity used for the Table 3, one draws a vertical thick dotted line 2% higher than the nominal value for the design point (320km/h, F m,320 ). One then qualitatively define the limits of values of /F m at (300km/h, F m,300 ) taken around this line. These limits are indicated in horizontal thin dotted lines. The range of variation thus defined is around [-8% 8%] around the nominal value. The result of 15% obtained in virtual certification is clearly out of range. One concludes that the random generation of irregularity is not very representative of the measurements used. One thus seeks to explain the origin of this gap. The last difference between simulations computed for virtual certification and these ones is the length of simulation; three catenary sections were used for the certification. Figure 13 shows the contact force at (300km/h, F m,300 ) for measured and nominal geometries. Three stretches have been defined to compute the relative deviation between the /F m of both signals. The first stretch is the same as the one computed for the variability study with a length of 400m. The relative deviation of 7% is included in the interval defined previously in Figure 12. The random geometry generation method thus includes the measured case.

The second stretch has the same length as the first one, namely 400m. It is located on and after an overlap section. The variation of /F m is 27%. The last stretch is the whole simulation, on which the variation of /F m is 17%, approximately the value given in Table 3 (15%). One observes that the criterion /F m is not precise enough to compare several geometries and the location of variations can be discriminating, particularly for the overlap section.

To confirm this conclusion, that the location has a strong influence on current collection quality, and justify further developments, Figure 14 shows the variation of /F m computed over 400m along the catenary. Higher values of local /F m are clearly located around overlap sections. This observation encourages developing more accurate criteria which take into account the location on a catenary section and avoid averaging critical areas.

Conclusion

The paper introduced variability in pantograph/catenary simulations due to three different types of irregularities and the following conclusions can be drawn. Wear irregularity has a low impact on current collection quality, even in the entirely determinist case defined here, which amplifies this impact. This kind of irregularity improves significantly the frequency coherence. Simulations show a positive impact on current collection. Taking this kind of irregularity into account is thus not a priority.

Aerodynamic irregularity only impacts the mean contact force, without inducing any significant effect over time. As the energy variation of the load is low, this kind of irregularity nearly does not affect current collection quality if the average is removed. Taking a range for the averaged load is thus representative of the impact of aerodynamic irregularity.

Geometry irregularity induces variations of amplitude that are representative of measurements for low frequencies, but does not improve coherence with measurements. The variation of the normative criterion σ/F m indicates that this kind of irregularity has a high impact on current collection, which can be good or bad. As this irregularity has the most impact on current collection, further studies on appropriate representation are clearly needed.

The last section first showed that speed has more influence than mean force on the impact of the geometric irregularity. It then appeared that the proposed random geometries might include non-representative cases, so that further work on geometry characterization is needed. Finally, current criteria used to define the current collection quality were shown to vary depending on location.

The relatively simple examples of this paper demonstrated the importance of introducing variability in pantograph catenary models. Quantification of uncertainty, introduction of other variable parameters, further sensitivity studies, more careful statistical interpretation of tests and redefinition of criteria are clear perspectives for future work.

Figure 1 . 2 . 2

 122 Figure 1. Filtered measurements (solid line) and nominal (dashed line) contact wire height (a) and spatial-spectrum of contact wire height (b)

Figure 3 .

 3 Figure 3. Comparison between contact wire height measured and computed by inverse static method: superposed heights (a) and differences of heights under droppers (b)

Figure 4 .

 4 Figure 4. Contact wire section of a new wire (a) and a worn wire (b), A being the area and G the centre of gravity

Figure 5 .

 5 Figure 5. Aerodynamic force spectrum measured (dotted lines) and spectrum shape defined (solid line)

Figure 8 .

 8 Figure 8. Frequency band analysis (left) and frequency coherence (right) of geometry (a), aerodynamic (b) and wear (c) irregularities.

Figure 9 .

 9 Figure 9. Comparison between inline tests and simulations with nominal geometry (a) and comparison between inline tests and simulations with measured geometry.

Figure 10 .

 10 Figure 10. Mean force (a), standard deviation (b), /F m (c) distributions for one thousand geometric irregularities at three different design points (speed, mean forces): 300km/h and F m,300 (black), 320km/h and F m,300 (grey), 320km/h and F m,320 (white).

Figure 11 .

 11 Figure 11. Dependence on (speed, mean force) of the criterion /F m for two groups of design points: (320 km/h, F m,300 ) against (320 km/h, F m,320 )(a) and (320 km/h ,F m,300 ) against (300 km/h, F m300 )(b). The nominal values are indicated by the solid lines.

Figure 12 .

 12 Figure 12. Dependence on (speed, mean force) of the criterion /F m between two design points, (300 km/h, F m,300 ) against (320 km/h, F m,320 ).

Figure 13 .

 13 Figure 13. Contact forces along three catenary sections at (300km/h, F m,300 ) for nominal (black) and measured (grey) geometries and relative deviations between /F m of measured and nominal geometries for three different zones.

Figure 14 .

 14 Figure 14. Measured contact forces at 320km/h (a) and local /F m defined from the contact force 200 m before and after the associated PK. Overlap sections are shown with vertical dotted lines.

Table 1 .

 1 Results of statistical values of sets of contact forces simulated

		M nom	M 3	Min 3	Max 3	S 3	S 1	S 2
	Geometric irregularity		0.2391	0.211	0.270	0.009	0.183	1.20
	Wear irregularity	0.2427	0.2377	0.233	0.242	0.003	0.01	0.34
	Aerodynamic irregularity		0.2441	0.240	0.250	0.002	0.628	0.181

Table 2 .

 2 Percentage of coherence for each type of irregularity

	Mean Coherence between	Simulations/	Nominal/	Measurements/	Simulations/
	[0-20]Hz	Measurements	Measurements	Measurements	Simulations
	Geometric irregularity	38.91%			64.95%
	Wear irregularity	52.37%	32.13%	58.11%	95.37%
	Aerodynamic irregularity	50.60%			95.10%

Table 3 .

 3 Deviation on  depending on the speed and presence of geometry irregularity

		Relative deviations on /F m
	Speed [km/h]	Nominal geometry /Inline	Measured geometry	Measured geometry /
		tests	/Inline tests	Nominal geometry
	260	-15%	-0.7%	14%
	280	-24%	-9%	17%
	300	-26%	-14%	15%
	320	-16%	-14%	2%

Acknowledgements

This paper describes work undertaken by SNCF Research and Innovation Department with the precious support of SDTools and Vibratec companies. The Authors would like to thank particularly Jerôme Antoni for his contribution on the signal-processing field.


Marc Hoffmann

Mauricio Labadie

Charles-Albert Lehalle

Gilles Pagès

Huyên Pham

Mathieu Rosenbaum

OPTIMIZATION AND STATISTICAL METHODS FOR HIGH FREQUENCY FINANCE

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

High Frequency Finance requires increasingly sophisticated algorithms in order to manage a series of orders from the point of view of an agent (a trader) that wishes to execute his or her orders in a limited amount of time, under various trading objectives and market constraints, possibly at very fine temporal scales. Section 1 (contributed by M. Labadie) gives a general framework for operating such a task, by means of a stochastic control program (using HJB equations) in order to deal with inventory constrainsts. Section 2 (contributed by Gilles Pagès) proposes a new on-line optimisation algorithm for posting passive orders in a limit order book. Section 3 (contributed by H. Pham) proposes a tick-by-tick price model that incoporates microstructure effects at fine scales and that diffuses at coarse scales and for which a stochastic control algorithm for executing orders is developed. Finally, Section 4 (contributed by C.A. Lehalle) offers a wider view for order book dynamics reviewing some recent approaches and suggesting an opening in the direction of the recent and attractive mean-field game theory.

HF trading under inventory and directional risk constraints

1.1. High-frequency trading and market-makers.

The high-frequency traders (HFTs) are trading algorithms, with the fastest response or lowest latency, that buys and sells assets in electronic and listed markets, trying to capitalise ephemeral arbitrage opportunities.

One special kind of HFTs are market-makers. HFTs offer liquidity to the market, i.e. they place both a buying and a selling limit order on the Limit Order Book (LOB). Since they are offering liquidity, their buying price (bid price) is smaller than their selling price (ask price); therefore, each time they buy and sell they earn the spread (difference between ask and bid prices). However, market-makers suffer execution risks since they cannot control when and how their orders will be executed.

Stochastic model and Hamilton-Jacobi-Bellman equation (HJB)

The optimisation programme for a market-maker is thus to maximise the profit by trading the spread over and over throughout the day whilst minimising both the execution and inventory risks. We extend the model of Avellaneda and Stoikov [START_REF] Avellaneda | High-frequency trading in a limit-order book[END_REF], where the only state variable is the mid-price S(t), to a general Markovian setting with several continuous state variables Y (t) = (S(t), Σ(t), Z(t), . . . ). In our case Z(t) is the half-spread of the market and Σ(t) the volatility of the mid-price. We use the convention of small letters for the values of the state variables at time t, e.g. y = Y (t). The market-maker places two limit orders of one share on each side of the LOB: a selling limit order at S(t) + δ + and buying limit order at S(t) -δ -. The arrival of buying and selling market orders to the LOB is assumed to be Poisson and independent; the distribution of how many levels in the LOB are comsumed by a single market order is exponential (see Avellaneda and Stoikov [START_REF] Avellaneda | High-frequency trading in a limit-order book[END_REF] for details).

Assume the utility function is the PNL (profit and loss) X(t) + S(t)Q(t). Since the only variables the market-maker can control is where to place her quotes, i.e. (δ + , δ -), then her value function is (see Fodra and Labadie [START_REF] Fodra | High-frequency market-making with inventory constraints and directional bets[END_REF][START_REF] Fodra | High-frequency market-making for multi-dimensional markov processes[END_REF])

u(t, y, q, x) = max (δ + ,δ -) E t,y,q,x [X(T ) + Q(T )S(T )] , (1)
where

E t,y [ • ] := E[ • | Y (t) = y].
The associated HJB equation of ( 1) is

(∂ t + L) u + max δ + ∈A Ae -k[z+δ + ] u(t, y, q -1, x + (s + δ + )) -u(t, y, q, x) (2)
+ max

δ -∈A Ae -k[z+δ -] u(t, y, q + 1, x -(s -δ -)) -u(t, y, q, x) = 0 , u(T, y, q, x) = x + qs ,
where L is the infinitesimal operator associated to the continuous state variables Y (t). The jump part of the equation ( i.e. the max terms) corresponds to the infinitesimal operator for the discrete state variable Q(t). In the presence of inventory-risk aversion ε and transaction costs α, the HJB equation takes the form (see Fodra and Labadie [START_REF] Fodra | High-frequency market-making with inventory constraints and directional bets[END_REF][START_REF] Fodra | High-frequency market-making for multi-dimensional markov processes[END_REF])

(∂ t + L) u + max δ + ∈A Ae -k[z+δ + ] u(t, y, q -1, x + (s + δ + ) -α) -u(t, y, q, x) (3)
+ max

δ -∈A Ae -k[z+δ -] u(t, y, q + 1, x -(s -δ -) -α) -u(t, y, q, x) = ενσ 2 q 2 , u(T, y, q, x) = x + sq -εηzq 2 .
There are similar models in the literature, where the controls are formally explicit in terms of a certain ODE system, but they can only be computed numerically. For more details see e.g. Cartea and Jaimungal [START_REF] Cartea | Risk metrics and fine tuning of high frequency trading strategies[END_REF] and Géant et al [START_REF] Guéant | Dealing with the inventory risk. a solution to the market making problem[END_REF].

Exact solution without inventory penalty

When there is no inventory penalty, the exact solution of the HJB equation ( 2) is (see Fodra and Labadie [START_REF] Fodra | High-frequency market-making with inventory constraints and directional bets[END_REF][START_REF] Fodra | High-frequency market-making for multi-dimensional markov processes[END_REF]) u(t, y, q, x) = u hold (t, y, q, x) + u mm (t, y) ,

u hold := x + q(s + ∆) , ∆ := E t,y [S(T )] -s u mm := 2 e E t,y T t A K e -KZ cosh(K∆) dξ , ψ * := δ + * + δ - * = 2 k (market-maker's spread) , r * := p + + p - 2 = s + ∆ (centre of spread) . (4)
From (4) we get the following observations. (i) The solution u has two strategies: one buy-and-hold u hold , which is the expected value of the portfolio if the trader keeps the position until the end, and one market-maker u mm , which is the gain from trading the spread very fast and frequently (hence the integral). (ii) Since u mm > 0, the market-maker gets a better PNL by trading the spread than a buy-and-hold investor. (iii) The market-marker's spread ψ * is constant and inveresly proportional to k: when k increases, the market orders are less likely to consume several levels of the LOB, hence the market-maker needs to offer better prices to get executed, which means a smaller spread. (iv) The centre r * of the spread is not the current mid-price s but E t,y [S(T )], i.e. the conditional expectation of the value of the mid-price. This means that the market-maker plays a price arbitrage: she tilts her quotes, accepting an inventory risk, in order to capitalise the current mispricing in the long run.

Asymptotic solution with inventory penalty

In the case of inventory penalty ε and transaction costs α, the solution and controls are not explicit. However, using perturbation methods on ε we can find the explicit expresion of the controls up to first-order (see Fodra and Labadie [START_REF] Fodra | High-frequency market-making for multi-dimensional markov processes[END_REF]):

ψ α * = 2 k + 2α + 2επ + O ε 2 (market-maker's spread), (5)
r α * = s + ∆ -ε E t,y T t H(α)πdξ + 2qπ + O ε 2 (centre of spread) , ∆ := E t,y [S(T )] -s , H(α) := 4 e Ae -K(Z+α) sinh(K∆) , π := ηE t,y [Z(T )] + νE t,y T t Σ 2 dξ .
From (5) we get the following observations. (i) The spread ψ α * has a linear correction via the unitary transaction cost α and the unitary inventory risk π. (ii) The current inventory q affects the centre r α * of the spread, compensating the directional bet via ∆. (iii) The transaction cost α makes the spread wider, but the total gain for trading the spread is constant and independent of α. (iv) If α > 0 then the market-maker will be less executed due to a wider spread. If everyone in the market have the same behaviour then the market spread widens, i.e. the available liquidity diminishes and the trading costs for market orders increases. (v) If α < 0 (as in the case of rebates in certain alternative stock markets) then the market-maker reduces her spread. If the rebate is big enough then the market-maker can have a zero-width spread, i.e. she trades one leg with a market order. With this strategy, called rebate arbitrage, all the profit comes from the rebate -α.

2. Optimal posting price of limit orders : learning by trading

Introduction

Our aim is to propose an on-line optimization posting method of passive orders in a limit order book, without specifying the limit order book dynamics. We consider an agent who wants to buy during a short period [0, T ] a quantity Q T of traded assets and we look for the optimal distance to the "fair" price (usually mid-, possibly best bid) where he has to post his limit order to minimize the execution cost (the case of a sell limit order may be obtained by symmetry). We assume the agent is a learning trader: he searches for the optimal posting distance (the price that balances adverse selection and non execution risk) by successive trials, errors and corrections. The derived recursive optimization procedure yields the best price adjustment to apply to an order for a given stopping time (reassessment dates) given past market observations.

We model the execution process of orders by a Poisson process (N

(δ)
t ) 0≤t≤T with intensity Λ T (δ, S) depending on the fair price (S t ) t≥0 and the distance of order submission δ. The execution cost results from the sum of the price of the executed quantity and a penalty function depending on the remaining quantity to be executed at the end of the period [0, T ]. This penalty κ Φ(Q) models the excess cost induced by crossing the spread and the resulting market impact of this execution. An optimal distance δ * is defined as a minimizer of the execution cost. This leads to an optimization problem under constraints which we solve by using a recursive stochastic procedure with projection (see [START_REF] Kushner | Stochastic approximation methods for constrained and unconstrained systems[END_REF] and [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF]). We prove the a.s. convergence of the constrained algorithm under additional assumptions on the execution cost function (and light assumptions on the fair price dynamics (S t ) t∈[0,T ] ). From a practical viewpoint, we give criteria on the model parameters that ensure the viability of the algorithm and the existence of a cost minimizer) relying on a functional co-monotony principle presumably satisfied by the fair price process.

Modeling and design of the algorithm

We focus our work on the problem of optimal trading with limit buy orders on one security. We only model the execution flow which reaches the price where the limit order is posted on a short period T with a Poisson process, namely

N (δ) T with intensity Λ T (δ, S) := T 0 λ(S t -(S 0 -δ))dt, (6)
where 0 ≤ δ ≤ δ max (δ max ∈ (0, S 0 ) is the depth of the order book), λ : [-S 0 , +∞) → R + is a finite nonincreasing convex function and (S t ) t≥0 is a càdlàg stochastic process modeling the dynamics of the the best opposite price of a security stock. Over the period [0, T ], we aim to execute a portfolio of size

Q T ∈ N invested in the asset S. The execution cost for a distance δ is E (S 0 -δ) Q T ∧ N (δ) T
. We add to this execution cost a market impact penalty function Φ : R → R + , non-decreasing and convex, with Φ(0) = 0, to model the impact effect (additional cost )induced by the execution of the remaining quantity

Q T -N (δ) T +
. The resulting cost function reads

C(δ) := E (S 0 -δ) Q T ∧ N (δ) T + κ S T Φ Q T -N (δ) T + , δ ∈ [0, δ max ] (7)
(where κ > 0 is a free tuning parameter). Prior to solving min 0≤δ≤δmax C(δ) by implementing a stochastic gradient descent numerically, we show that under natural assumptions on Q T and κ, the function C is twice differentiable and strictly convex on [0, δ max ] with C (0) < 0 and that the derivatives C (k) , k = 0, 1, 2, have a representation as an expectation. Consequently,

argmin δ∈[0,δmax] C(δ) = {δ * }, δ * ∈ (0, δ max ] and δ * = δ max iff C is non-increasing on [0, δ max ].
The representation of C as an expectation means that there exists a Borel functional H : [0,

δ max ]×D ([0, T ], R) → R such that C (δ) = E H δ, (S t ) t∈[0,T ] , δ ∈ [0, δ max ]
where D ([0, T ], R) → R denotes the set of càdlàg functions from [0, T ] to R (see [START_REF] Laruelle | Optimal posting price of limit orders: learning by trading[END_REF] for an explicit form for H). In practice to implement numerically the recursive procedure, we have to replace the "copies" S (n) by copies S(n) of a time discretization S = ( Sti ) 0≤i≤m , typically an Euler scheme of a Brownian diffusion process. Then, with an obvious abuse of notation for the function H, we can write the implementable procedure as follows:

δ n+1 = Proj [0,δmax] δ n -γ n+1 H δ n , S(n+1) ti 0≤i≤m , n ≥ 0, δ 0 ∈ [0, δ max ], (8)
where S(n) ti 0≤i≤m are copies of ( Sti 0≤i≤m either independent or sharing "ergodic" properties, namely some averaging properties in the sense of [START_REF] Laruelle | Stochastic approximation with averaging innovation applied to finance[END_REF]. In the first case, one will think about simulated data after a calibration process and in the second case to a direct implementation using a historical high frequency database of best opposite prices of the asset S.

Several a.s. convergence results are established in [START_REF] Laruelle | Optimal posting price of limit orders: learning by trading[END_REF] under various types of assumptions on the sequences ( S(n) ti ) 0≤i≤m : a first setting is to assume that the ( S(n) ) are i.i.d. which correspond in the real world to the Monte Carlo simulation of a formerly calibrated model for the dynamics of S; in a second setting we assume that we work with a historical high frequency dataset of prices of the asset assumed to share rather light averaging properties. In both cases we show that δ n a.s.

-→ δ * . Several numerical experiments on a true dataset ar reproduced and analyzed in [START_REF] Laruelle | Optimal posting price of limit orders: learning by trading[END_REF].

Criteria for the convexity and monotony at the origin

In this section, we look for simple criteria involving the parameter κ, that imply the requested assumption on the execution cost function C (or C) without never needing to really specify the process S. The original form of the criteria implying that C (0) < 0 and C ≥ 0 cannot really be used in practice since they involve ratios of expectations of functionals combining both the dynamics of the asset S and the execution parameters in a highly nonlinear way. The key establishing the criteria is the functional co-monotony principle established in [START_REF] Laruelle | Optimal posting price of limit orders: learning by trading[END_REF][START_REF] Pagès | A functional co-monotony principle with an application to peacoks and barrier options[END_REF] for a wide class of diffusions and their associated time discretization schemes.

Definition 2.1 (Co-monotony). A Borel functional

F : D([0, T ], R) → R is non-decreasing (resp. non- increasing) if ∀α, β ∈ D([0, T ], R), (∀t ∈ [0, T ], α(t) ≤ β(t)) =⇒ F (α) ≤ F (β) (resp. F (α) ≥ F (β)).
Two functionals F and G on D([0, T ], R) are co-monotonic if they have the same monotony. A functional F has polynomial growth if

∃ r > 0 s.t. ∀α ∈ D([0, T ], R), |F (α)| ≤ K (1 + α r ∞
) . Definition 2.2 (Functional co-monotony principle). A càdlàg (resp. continuous) process (S t ) t∈[0,T ] satisfies a functional co-monotony principle if for every pair F, G : D([0, T ], R) → R of co-monotonic P S -a.s. continuous Borel functionals (for the sup-norm over [0, T ]) with polynomial growth such that F (S), G(S) and F (S)G(S) ∈ L 1 , one has

E F (S t ) t∈[0,T ] G (S t ) t∈[0,T ] ≥ E F (S t ) t∈[0,T ] E G (S t ) t∈[0,T ] . (9)
Among processes satisfying such a functional monotony principle, the (fractional) Brownian motion, Brownian one-dimensional diffusions, Liouville processes Lévy processes, etc (see [START_REF] Pagès | A functional co-monotony principle with an application to peacoks and barrier options[END_REF]).

Theorem 2.3. Assume that (S t ) t∈[0,T ] satisfies a functional co-monotony principle. Assume that λ(x) = Ae -kx , x ∈ R, A, k ∈ (0, +∞). Then the following monotony and convexity criteria hold true.

(a) Monotony at the origin:

C (0) < 0 as soon as Q T ≥ 2T λ(-S 0 ) and κ ≤ 1 + kS 0 kE [S T ] (Φ(Q T ) -Φ(Q T -1)) . (b) Convexity. Let ρ Q ∈ 0, 1 -P(N µ =Q T -1) P(N µ ≤Q T -1) |µ=T λ(-S0) . If Φ = id, assume that Φ satisfies ∀x ∈ [1, Q T -1], Φ(x) -Φ(x -1) ≤ ρ Q (Φ(x + 1) -Φ(x)). If Q T ≥ 2T λ(-S 0 ) and κ ≤ 2 kE [S T ] Φ (Q T ) , then C (δ) ≥ 0, δ ∈ [0, δ max ], so that C is convex on [0, δ max ] or
for the Euler scheme of (S t ) t∈[0,T ] with step T m , for m ≥ m b,σ of a diffusion with constant diffusion coefficient. (c) The same inequalities hold for tany R m+1 -time discretization sequence (S ti ) 0≤i≤m which satisfies a comonotony principle.

Semi Markov model for market microstructure and HF trading

Markov renewal model for stock price

The dynamics of the mid-price are given by the càd-làg piecewise constant process:

P t = P 0 + 2δ Nt n=1 J n ,
where P 0 is the opening mid-price, 2δ > 0 is the tick size, N t is the point process associated to the tick times (T n ) n (i.e. the price jump times), and (J n ) n is the marks sequence valued in {+1, -1} indicating whether the price jumped upwards (J n = +1) or downwards (J n = -1) at time T n , called price direction. We use a Markov renewal approach to model the marked point process (T n , J n ). (J n ) is an irreducible Markov chain with symmetric transition matrix of diagonal terms 1+α 2 , and antidiagonal terms 1-α 2 , where α ∈ [-1, 1) represents the correlation between two consecutive price returns. Estimation on real data leads to a negative value of α, meaning that the stock price exhibits a short-term mean-reversion, which is a well-known stylized fact about high-frequency data, also called microstructure noise. Denoting by S n = T n -T n-1 the sequence of inter-arrival times associated to (N t ), we assume that conditionally on the sign J n J n+1 = ±1 of two consecutive price directions, (S n ) n is an independent sequence with distribution: F ± (s) = P[S n+1 ≤ s|J n J n+1 = ±1], and density f ± (s). The distinction between F + and F -models the clustering of trading activity, that is the irregular spacing of tick times according to the mean-reverting or trend period of price jumps. We denote by F = 1+α 2 F + + 1-α 2 F -the unconditional distribution of (S n ). The intensity function of price jump in the same (resp. opposite) direction is then given by:

h + = 1+α 2 f+ 1-F (resp. h -= 1-α 2 f- 1-F ).
Markov embedding. Let us define the last price jump direction: I t = J Nt , and the elapsed time since the last jump: S t = t -sup{T n : T n ≤ t}. Then, the price process (P t ) is embedded into a Markov system (P t , I t , S t ) with three observable state variables. Diffusive behavior at macroscopic scale. The scaled price process P (T ) t = P tT √ T , t ∈ [0, 1], exhibits a diffusive behavior given by the following functional central limit theorem: lim

T →∞ P (T ) (d) = σ ∞ W,
where W is a Brownian motion, and

σ 2 ∞ = 1 ∞ 0 tdF (t) 1+α
1-α is the macroscopic variance. Estimation, simulation and more properties of this Markov renewal model are studied in [START_REF] Fodra | Semi-markov model for market microstructure[END_REF].

Market order flow modeling

We model the market order flow of the Limit Order Book (LOB) by a marked point process (θ k , Z k ) k , where (θ k ) k represents the arrival times of the market order, and (Z k ) k , valued in {±1} is the side of the exchange in the LOB: Z k = +1 (resp. -1) means that the trade is exchanged at the best bid (resp. ask) price, thus associated to a market sell (resp. buy) order. The existing dependence between the market order flow and the stock price is modeled as follows. The counting process (M t ) t associated to (θ k ) k is a Cox process with conditional intensity λ(S t ), where λ(s) is a decreasing function, representing the fact that trade arrivals should occur less often in a period of price stability, i.e. when S t is large. The dependency between market order trade and stock price is modeled through the relation:

Z k = Γ k I θ k -, where (Γ k ) k is an i.i.d. sequence on {±1} with mean ρ ∈ [-1, 1]. This means that ρ = Corr(Z k , I θ - k
), with the interpretation that when ρ > 0 (resp ρ < 0) most of the trades arrive on the strong (resp. weak) side +I θ k -(resp. -I θ k -) of the LOB. Estimation on real data usually leads to negative value of ρ, which formalizes the observed fact that most of the trades arrive in the unfavorable side of the LOB for the market maker, inducing the so-called adverse selection risk.

The market making application

The agent strategy consists in placing continuously limit orders of small size L at the best available price, while respecting an inventory constraint: her control is described by a process ( ± t ) t valued in {0, 1} with the convention that when ± t = 1, a limit order of size L is posted at time t on the side ±I t-. When a market order arrives, she is executed on the side ±I t-of the LOB, of ± t times a random variable distributed as ϑ ± L (dk), with support on {0, . . . , L}. The objective of the market maker is to maximize over her limit orders strategies the expected gain:

E P N L T -CLOSE(Y T ) -η • RISK 0,T , (10)
where P N L T = X T + Y T P T is the portfolio value at terminal date T , evaluated at the mid price P , with X the cash holdings, and Y the inventory constrained to lie in a given bounded set Y, CLOSE(Y T ) = -(δ + ε)|Y T |, is the liquidated value of the inventory at terminal date with ε ≥ 0 a fixed fee paid at each transaction, and RISK 0,T = T 0 Y 2 t • d[P ] t represents the market risk by holding the inventory during the whole trading period, weighted by the penalization parameter η ≥ 0. The value function v to this Markovian stochastic control on finite horizon [START_REF] Gatheral | No-dynamic-arbitrage and market impact[END_REF], involving the state variables (t, P t , I t , S t , X t , Y t ), can be written in the reduced-form: v(t, p, i, s, x, y) = x + yp + ω(t, s, yi), where ω(t, s, q) is the unique bounded viscosity solution to:

∂ t + ∂ s ω + 2δ(h + -h -)q -4δ 2 η(h + + h -)q 2 + max q-L∈Y N + trd [ ] + N + jmp [ ] ω + max q+ L∈Y N - trd [ ] + N - jmp [ ] ω = 0, ω(T, s, q) = -|q|(δ + ε),
and the operators N 's are given by:

N ± trd [ ] ω := 1 ± ρ 2 λ(s) [ω(t, s, q ∓ k ) -ω(t, s, q) + k (+δ -ε)] ϑ ± L (dk) N ± jmp [ ] ω := h ± (s) [ω(t, 0, ±q -L ) -ω(t, s, q) + L (-δ -ε)] .
The solution (value function and optimal control) of the market making problem is then completely characterized in terms of an integro-ordinary differential equation involving three observable state variables; current time t, elapsed time s from last jump, and strong inventory q = yi. We provide a convergent numerical scheme for solving this IODE, and illustrate our results with some computational tests displaying the shape of the optimal policies, see [START_REF] Fodra | High frequency trading in a markov renewal model[END_REF].

Dynamical models for orderbook dynamics: which ones and what for

This work is mainly based on the paper Efficiency of the Price Formation Process in Presence of High Frequency Participants: a Mean Field Game analysis, a joint work with PL Lions, JM Lasry and A Lachapelle [START_REF] Lachapelle | Efficiency of the Price Formation Process in Presence of High Frequency Participants: a Mean Field Game analysis[END_REF].

The emergence of a need in orderbook dynamics modelling

To manage the increase of complexity following the multiplication of venues and dispersion of liquidity across all trading venues available, market participants more intensively used electronic trading and quantitative (i.e. optimal) slicing of large orders [START_REF] Almgren | Optimal execution of portfolio transactions[END_REF]. Besides, the need for cross-venue arbitrage of latency and liquidity increased, giving birth to a new kind of middle-men: the High Frequency traders [START_REF] Menkveld | High Frequency Trading and The New-Market Makers[END_REF]. A better understanding of the orderbook dynamics improves: the understanding by regulators and policy makers of the contribution of each market participant to the efficiency of the PFP (Price Formation Process); the existing market impact models ( [START_REF] Gatheral | No-dynamic-arbitrage and market impact[END_REF] Recently, some papers focused on this topic: a Dynamic Model of the Limit Order Book [START_REF] Roşu | A Dynamic Model of the Limit Order Book[END_REF]; a Fokker-Planck description for the queue dynamics of large tick stocks [START_REF] Gareche | A Fokker-Planck description for the queue dynamics of large tick stocks[END_REF]; Price Dynamics in a Markovian Limit Order Book Market [START_REF] Cont | Price Dynamics in a Markovian Limit Order Book Market[END_REF]; Measuring the resiliency of an electronic limit order book - [START_REF] Large | Measuring the resiliency of an electronic limit order book[END_REF]; State-Dependent Model of Order Book Dynamics (forthcoming, by Huang, Rosenbaum, Lehalle).

The Mean Field Game approach as a synthesis

MFG (Mean Field Games) are designed to model this type of dynamics [START_REF] Lasry | Mean field games[END_REF]. An MFG model of endogenous volatility has been proposed in "Large investor trading impacts on volatility" [START_REF] Lions | Large investor trading impacts on volatility[END_REF]. An MFG-based model has already been designed at a larger time scale in "High Frequency Simulations of an Order Book: a Two-Scales Approach" [START_REF] Lehalle | High Frequency Simulations of an Order Book: a Two-Scales Approach[END_REF]. We wanted to provide a framework at the smallest one. In this paper we model participant behaviours in the orderbook; obtain "master equations" describing the emerging dynamics; use them to describe the outcome of the mixing of HFT and Institutional Investors in the same orderbook.

One queue modelling

Sellers only, they arrive in the "game" according to a Poisson process of intensity λ, they compare their value function u(x) (where x is the size of the queue) to zero to choose to wait in the queue (when u(x) > 0) or not, the queue is consumed by a Poisson process of intensity µ(x), in case of transaction, a "pro-rata" scheme is used ("equivalent" to infinitesimal possibility to modify orders): q/x of the order is part of it; relaxed for FIFO in few slides. Dynamics. The value function evolves following these four main possible events (we have price impact P (x) and waiting cost proportional to c q): u(x) = (1 -λ1 u(x)>0 dt -µ(x)dt) • u(x)

← nothing happens + λ1 {u(x)>0} dt • u(x + q) ← new entrance + µ(x)dt • q x P (x) + (1 -q x )u(x -q) ← transaction -c q dt ← waiting cost

To solve it at the kth order, we will perform a Taylor expansion of u(x + q) and u(x -q) for small q at the kth order). At the second order, we obtain: 0 = µ(x) x (P (x) -u) -c + (λ1 {u>0} -µ(x))u + q 1 2 (λ1 {u>0} -µ(x))u + µ(x) x u .

This corresponds to a (trivial) shared risk Mean Field Game monotone system with N = 1. The mean field aspect does not come from the continuum of agents (for every instant, the number of players is finite), but rather from the stochastic continuous structure of entries and exits of players. Solution for a specific form of µ(x). At queue sizes x * such that x * = µ(x * )P (x * )/c, u sign changes. Moreover, for the specific case µ(x) = µ 1 1 x<S + µ 2 1 x≥S : There is a point strictly before S where u switches from negative to positive. It means that participants anticipate service improvement.

Going further

In the paper [START_REF] Lachapelle | Efficiency of the Price Formation Process in Presence of High Frequency Participants: a Mean Field Game analysis[END_REF], a model with two coupled queues is developed and analysis in a first step, and then extended to more than one class of trading agents. We then use this model to understand the mechanism of liquidity dynamics in a limit orderbook when traditional investors are mixed with high frequency traders. Such results can help regulators and policy makers to understand the stakes and the potential consequences of their decisions to adjust the market microstructure (see [START_REF] Lehalle | Market Microstructure in Practice[END_REF] for more details about practical aspects).

  [START_REF] Almgren | Direct Estimation of Equity Market Impact[END_REF] [21][START_REF] Lillo | Econophysics -Master Curve for Price -Impact Function[END_REF]), a key component of the PFP dynamics; and the capability to use more sophisticated optimal execution schemes ( [4][START_REF] Pagès | Optimal split of orders across liquidity pools: a stochatic algorithm approach[END_REF] [13][START_REF] Guéant | Dealing with the inventory risk: a solution to the market making problem[END_REF]).

Figure 1 .

 1 Figure 1. All the components of the one queue model.


Lech Pawłowski

Review of solution precursor spray techniques applied to obtain ceramic films and coatings

The review describes the techniques of film and coatings deposition using solution as a feedstock. After a short description of the processes of: (i) spray pyrolysis; (ii) solution precursor plasma spray (SPPS); and, (iii) solution precursor high velocity oxy-fuel, some properties of solution feedstock are discussed. The properties include the chemical composition of solutions used to synthesize typical coatings, which determine the flow of solution in a pipeline and its injection to plasma or flame. The chemical and physical phenomena occurring at spraying are described including atomization in gas, heating and vaporization of solution, precipitation of solution and formation of solid and heating of this solid. The molten solid is then splashed onto substrate. Typical coatings formed using solution precursor spraying and the applications of solution precursor thermally sprayed coatings are reviewed.

Introduction

The spray pyrolysis started to be used as early as in the 60ties of the last century [1]. The method was used to synthesize films of metal oxides, sulfides and selenides. The use of the torches generating combustion flames and plasma jets to process liquid solutions was first applied at the end of the 20th century [2]. The authors were inspired by the spray pyrolysis process. The torches generating high velocity flames have started to be used more recently [3]. The use of liquid feedstock in thermal spraying instead of usual solid powders was inspired by an interest in obtaining coatings having fine microstructure. Such coatings were easier to be sprayed with liquid feedstock (solution or suspensions) than with fine solid powders. In addition, thermal spraying is a well-established technology. The use of this technology for producing such coatings is industrially advantageous because of a dense network of job shops having expertise in coatings deposition, high productivity of coatings manufacturing and also a dynamic and well organized community of involved professionals.

Solution spray techniques 2.1 Spray pyrolysis

The technique of spray pyrolysis consists of atomizing the solution directed towards a hot substrate. The solutions of metal salts dissolved in organic or water are used generally to obtain coating of metal oxides. The atomized droplets impact the heated substrate to form lamellas and undergo a thermal decomposition. The substrate temperature is the main parameter, which influences the morphology and properties of coatings [4]. The final coating is a superposition of overlapping lamellas.

2.2

Solution precursor plasma spraying The d.c. torches generate arc which heats working gases to form plasma jets. The jet's temperature, in a conventional plasma torch, reaches 14 000 K and its velocity on the nozzle exit reaches 800 m/s. The torches used in SPPS process [5]: (i) conventional, one cathode one having radial introduction of solution, (ii) three-cathode torch having axial introduction of solution; and (iii) segmented anode torches. The injection of solution can be made by an atomizer or by a nozzle. The working gases used to generate the plasma is usually argon with a molecular gas such as hydrogen or nitrogen. The molecular gases are added in order to increase the thermal conductivity of the plasma. The conventional one-cathode torches have important fluctuations and instabilities depending on kind of working gas [6]. The molecular gases contribute considerably the voltage arc fluctuations in the kilohertz frequency range. The use of three-cathode torches without the use of molecular gases maybe a solution of the fluctuation problem [7].

2.3

High velocity oxy-fuel spraying The high velocity oxy-fuel torch is fuelled by a combustion hydrocarbon gases as ethylene, propylene and many others or hydrocarbon liquids as kerosene which is atomized prior to combustion. The combustion takes place in a chamber under high pressure, which varies from 0.3 to 4 MPa. The resulting jet issuing from a nozzle may reach temperature of about 3000 K and supersonic velocity by 2000 m/s. The solution can be injected radially, externally to the gun, directly into the HVOF flame or axially, directly into the combustion chamber [3].

3

Preparation of solutions The precursors used for solution spraying of metal oxides films and coatings are generally the salts of metals and very frequently nitrates and acetates. In this way it is possible to obtain coatings of multi-oxide compound such as e.g. YBCO or YIG. Some authors, like [8], used the nitrates as the initial products to start sol-gel procedure to obtain fine particles. Such additives as acetic acid modify the chemistry of solution and influence the morphology of coatings. For example the structure of TiO 2 was changing from cracked to crack-free after the introduction of acetic acid into the solution [4]. The salts are used as solutions in water or in ethanol.

4

Phenomena occurring at spraying The phenomena occurring during spray pyrolysis can be summarized by taking an example of spraying of metal salt solution onto a heated substrate. The most important of them is the deformation of the liquid into droplet and its thermal decomposition leading to the formation of the desired compound on the substrate [4]. The phenomena inflight include mainly aerodynamic breakdown accompanied by the chemical reaction that occurs on the substrate.

The solution precursor thermal spraying includes more phenomena occurring in-flight occurring before the impact with the substrate (see Fig. 1): a) aerodynamic break up; b) heating, vaporization and internal precipitation; c) internal pressurization and droplet breaking-up; d) solid particle sintering, heating; e) melting and evaporation from the melt.

Fig. 1. Phenomena which occur during a flight of a solution

droplet in a high temperature gas jet [9].

Injection

The solution may be transported to the injector with the use of a peristaltic pump or by a pneumatic transport. Its flow rate is controlled by the rotation of motor rollers or by the pressure acting on the fluid, respectively. The solution maybe atomized as in spray pyrolysis process or injected as a liquid stream to the jet or flame. The atomization is a process of which the external energy phase is needed to break up the liquid bulk. The liquid stream undergoes an aerodynamic break up and the ligaments and large droplets are formed. Then, depending on the drag force acting on the droplet and their surface tension the large droplet may disintegrate into smaller ones. The process is called secondary break up [10]. This process leads to the formation of finer microstructure.

4.2

Physico-chemical phenomena The droplets formed by successive break ups undergo rapid vaporization of solvent in hot temperature flame or jet. Consequently, the concentration of solute at the droplets surface increases, what may lead to precipitation in this region (because the concentration of solute reaches the super-saturation limit). The process is known as homogeneous precipitation [9] and was investigated for different thermal spray processes by e.g. [9,11] and many others. The formation of solid crust outside of heated liquid leads to an increase in internal pressure inside a shell. Under the action of the pressure the shell may get fractured and smaller droplets of liquid get free and the process of precipitation and shell formation may be repeated. The processes may lead finally to the formation of solid particle as shows it Fig. 1. The chemical phenomena occurring in flight of liquid droplet depend on the used metal salt precursor and on temperature. Generally, the sequence of chemical reactions includes evaporation of solvent, precursor pyrolysis and crystallization [12] 5

Possible applications of coatings The solution precursor spray techniques including spray pyrolysis and thermal spray techniques are at present intensively explored. The potential advantage of these technologies can be related to the purity of the initial liquid precursors, which enable chemically pure coatings to be synthesized. The particularity of processing of small nanoparticles renders it possible to obtain such phases as α-Al 2 O 3 , or TiO 2 crystallized as anatase obtained which are difficult to obtain using conventional thermal spray technique. Finely grained microstructure has another advantage, which is small thermal conductivity useful while synthesizing thermally insulating coating for TBCs. The capacity to coat large surfaces must be also underlined. Not but not least, the possibility of obtaining columnar microstructure, is useful for enhancing the thermal shock resistance of coatings. The properties of solution precursor sprayed coatings have been tested by now for the following applications: a) SOFC as electrodes, i.e. Ni-YSZ as anode, V 2 O 5 as cathode and YSZ as solid electrolyte; b) TBC made of YSZ and of ZrO 2 .La 2 O 5 ; c) optoelectronic applications of Eu:Y 2 O 3 or CZST coatings; d) dielectric coatings in electronics made of α-Al 2 O 3 ; e) photocatalysis of TiO 2 coatings; f) magnetic applications of YIG coatings; g) biomedical coating of TiO 2rutile coatings. The further perspectives of application these coatings in industry are related to the better control of the deposition process.