Julien Nespoulous

Christian Soize

Christine Funfschilling

Guillaume Perrin

Optimization of train speed to limit energy consumption

Keywords: Optimization under constraints, high-speed train, energy saving

The speed profile of a train plays an important role in energy consumption and resulting costs. The industrial objective of this work is thus to develop a method to reduce the energy consumed by a train over a journey by playing on the driver commands (traction and braking forces) while respecting punctuality constraints. A coupling between measured data and simulation is proposed to solve this optimization problem. First, a rigid body approach (Lagrangian formalism) is introduced to characterize the dynamics of each element of the train and their interactions with their environment. In particular, the aerodynamic (including the wind effect), traction, and braking forces are taken into account, and a special attention is paid to the vertical and lateral characteristics of the track as they play a key role on the train dynamics. Secondly, a model for energy consumption and recovery (thanks to dynamic braking) is introduced. Experimental measurements of a high-speed line are then used to estimate the parameters on which the two previous models are based and to validate their predictive capacities. The optimization problem under constraints is finally solved using an evolutionary algorithm where the constraints are implemented using an augmented Lagrangian formalism. The performance of the proposed method in terms of speed optimization and energy consumption reduction is compared to measurements associated with commercial trains.

Introduction

The electric energy consumed by a train greatly depends on the driver commands. Indeed, we can measure variations of the energy consumed of around 20% between two drivers on the same journey. To limit this, some crossing points indicate at which speed and at which time the driver should pass a position. Nevertheless, important variations of the energy consumed are still noticed. Thus, the objective of this work is to propose a strategy to help the driver to find the train controls that will allow him to best meet the punctuality constraints while consuming as little energy as possible. This would particularly be useful when traffic is disrupted.

Several papers have addressed related objectives. Most of them are based on the coupling of two simplified models: one for the train dynamics and one for the energy consumption. In terms of dynamic modeling, rigid body approaches combined with a Lagrangian formalism are mainly considered [START_REF] Bosso | Long train simulation using a multibody code[END_REF], as well as the projection of the dynamics equations on the longitudinal to the train axis [START_REF] Wu | Longitudinal train dynamics: an overview[END_REF]. As for the modeling of the consumed energy, more or less sophisticated CONTACT Julien Nespoulous. Email: julien.nespoulous@sncf.fr models can be found in the literature to fit with the reality [START_REF] Bosquet | Modélisation énergétique et identification des trains pour l'écoconception des lignes ferroviaires à grande vitesse[END_REF]. The same framework will be used in this work for physical modeling.

These models can then be used to minimize the energy consumed by the train during its journey. The main difficulty lies in the fact that the dynamics and the constraints are nonlinear. To speed up the resolution of the nonlinear dynamic problem, [START_REF] Wang | Multiple-phase train trajectory optimization with signaling and operational constraints[END_REF] proposes to use the pseudo-spectral method, which projects the continuous equations on orthogonal polynomials. The optimization problem is then solved thanks to a nonlinear programming algorithm. Another approach is to translate the nonlinear dynamic problem and constraints into a mixed integer linear programming problem, which is easier to numerically solve with a linear programming solver [START_REF] Wang | Optimal trajectory planning for trains using mixed integer linear programming[END_REF]. Genetic or dynamic algorithms [START_REF] Scheepmaker | Review of energy-efficient train control and timetabling[END_REF] were also proposed to manage these non-linearities. In this work, we have chosen to use these evolutionary algorithms because of their general and flexible nature. This type of energy optimization has been used to more precisely study different driving choices: the braking strategy and its impact on energy recovery [START_REF] Martinis | Models and Methods to Optimise Train Speed Profiles with and without Energy Recovery Systems: A Suburban Test Case[END_REF], but also the arbitration between delay-recovery, energy-saving or on-time strategy [START_REF] Wang | Multi-train trajectory optimization for energy efficiency and delay recovery on single-track railway lines[END_REF]. Finally, there are works integrating the minimization of rail and wheel wear. Indeed the modifications and the potential multiplication of braking and acceleration phases induced by the minimization of energy could prematurely degrade the system [START_REF] Conti | Energy and wear optimisation of train longitudinal dynamics and of traction and braking systems[END_REF].

In this work, we present a method that distinguishes itself from these different works for several reasons. First, the optimization method is applied to the driver commands and not directly on the speed trajectory. Thus, the method can directly be transposed to real high-speed trains. Then, the constraints combine both punctuality (applied on the end speed and position of the train) and physical constraints (applied all along the journey). The model includes both pneumatic and dynamic braking, which is able to recover energy. Those modifications allow a better description of the dynamic behavior of the train. Finally, we use an iterative method to limit the number of evaluations of the cost function, which would otherwise be far too long. A validation process has been carried out by comparing the simulation results with measured data in real traffic conditions.

The implementation of a real-time driver command must be based on the use of a robust model, a fast optimization algorithm, and an appropriate computer. The robust aspects with respect to uncertainties are in progress and are not presented in this paper. The optimization algorithm that is presented is already fast because it is based on the use of a surrogate nonlinear dynamical model (digital twin model) constructed from the full detailed dynamical model of the train. In the optimization algorithm, the computation devoted to each draw of the population of points (last paragraph of Section 3.2) is parallelized. For each draw the surrogate nonlinear dynamic model is used (second paragraph of Section 3.2) and requires, using one core of a CPU, about 30 seconds for the system under consideration (the one described in Section 2.1). Using a GPU instead of a CPU would yield a negligible calculation time with the GPU and consequently, the real-time optimization will be able to be performed without difficulty with an appropriate computer configuration (several GPU's in parallel).

The organization of the paper is the following. Section 2 describes the models that will be used to model the high-speed train longitudinal dynamics and its energy consumption. The optimization problem, its numerical implementation, and the results it can provide, are then presented in Section 3. At last, Section 4 concludes and proposes several prospects for this work.

High-speed train dynamics and energy consumption

This section presents the system of interest (a train running at high speed on a dedicated railway track) and proposes models to compute its dynamic behavior and the energy it consumes. The predictive character of these different models is then evaluated by confronting the simulation results with in-line measurements.

System description

In this work, we use a model of a TGV-Dasye that is a commercial high-speed train running on the French Rhin-Rhone high-speed line. This train is composed of eight passengers cars pulled by two engines at the head and tail of the train. Except for the two motor cars, bogies are shared between carriages. The engine bogies are called motor bogies to dissociate them from carrying bogies. Each bogie is composed of two wheelsets. These three bodies (cars, bogies, and wheelsets) are integrated in a rigid body framework to describe their dynamic behavior.

There are three kinds of interactions imposed on the system: the internal, the external, and the driver command interactions.

First, the rigid bodies are linked to each other by suspension joints. These joints can take different forms with linear (springs, dampers) but also nonlinear (bumpstops, bushes, pinlinks elements ...) behaviors. These behaviors are assumed well known and described by their behavior laws. We consider these suspensions as internal to the system because it consists in interactions between rigid bodies.

The external forces regroup all the interactions between the train and its environment. It includes the weight, the aerodynamic forces, and the contact forces. The weight is expressed on each rigid body thanks to the mass and inertia tensor. The aerodynamics can be applied to each moving rigid body, but we neglect it on bogies and wheelsets. In particular, the aerodynamic forces take into account the effect of wind, which has a big influence on the energy consumed by the train [START_REF] Quost | Modélisation de l'effet du vent sur les trains à grande vitesse[END_REF]. This wind speed is generally decomposed in two terms: the mean speed of wind v w (t), which is forecasted at different positions of the track, and the random fluctuations around this mean value (that can be modeled by a stochastic process). In this work, these fluctuations are neglected and we only keep the mean speed of wind. The mean speed of the wind is then projected on the longitudinal axis defined by the track and can be used to obtain a better approximation of the aerodynamic force. An example of the wind profile that has been used and its projection on the track axis is shown on Figure 1. To model the wheel/rail interaction, the Hertz model is chosen to describe the normal forces and the contact creep forces are approximated thanks to Kalker's wheel/rail theory (see [START_REF] Liu | Comparison of wheel-rail contact models in the context of multibody system simulation: Hertzian versus non-Hertzian[END_REF][START_REF] Spiryagin | Creep force modelling for rail traction vehicles based on the Fastsim algorithm[END_REF] for more details about these contact forces).

Finally, the command forces correspond to traction and braking, as their amplitude is controlled by the driver. Two types of braking are however distinguished, the dynamic and pneumatic braking, because they do not apply to the same subsystems and do not have the same impact on the energy consumed. In the same way as for the traction forces, the dynamic braking is applied to the motor wheelsets and can be simply described as the inversion of the motors. It is therefore particularly interesting as it is able to recover energy and inject it in the catenary. On the contrary, the pneumatic braking only dissipates energy but is applied to each carriage wheelsets.

The driver command u(t), which indirectly controls the train speed, can be modeled by a continuous time-dependent function representing the traction and braking normalized torque. This command is a mathematical modeling of the traction/braking manipulator used by the train drivers. The traction and braking forces capacity are limited by the motor and braking Wind profile (km/h) capacities. In the proposed model, the driver command u(t) takes its values in the interval [-1, 1]. The bound -1 (resp. 1) is the maximum value for the braking (resp. traction) forces, which cannot be exceeded in the model. These maximal traction and braking capacities depend on the train speed. These bounds correspond to longitudinal accelerations that cannot induce passengers comfort and train couplers problems. The jerk of the train is controlled by the time derivative of u. In this version of this work, no constraint has been introduced for limiting the jerk. Nevertheless, we have checked that for the optimal driver command found, this derivative stays sufficiently small. If u(t) is positive (respectively negative), the driver uses the motors (respectively brakes). The dynamic braking is prioritized over the pneumatic braking following the braking strategy of in-line trains. Thus, the pneumatic braking is only used for highamplitude braking (for example for emergency stops). For instance, at a time t, if u(t) = -0.8 the driver is braking using 80% of the maximal braking capacity. For u(t) < 0, the traction torque is equal to zero. By defining this control parameter instead of directly working on the train speed, as it is generally done, we get as close as possible to the behavior of the driver who can choose when to use motors or brakes. In addition to being novel, this approach will also facilitate the integration of the obtained results into the on-board computers.

Models for the train dynamics and the energy consumption

Dynamics equations

A nonlinear 3D dynamic model is constructed from the Lagrangian formalism using a multibody approach applied to each rigid body (8 passenger cars, 2 motor cars, 13 bogies, and 26 wheelsets). All the forces that have been listed in Section 2.1 are taken into account in the model (including aerodynamic, contact, and suspension forces) and are adapted to the TGV-Dasye train under consideration. As the objective of this work is to reduce energy consumption under punctuality constraints and speed limitations, the quantities of interest we need to extract from the train dynamics are the evolution of the speed and the energy consumption during the simulated journey, but also the arrival time and the final position. An important point is that all these quantities can be derived from the longitudinal dynamics of the train only, which invites us to directly consider the behavior of the whole train. In this work, we have chosen to study the longitudinal dynamics of the whole train whose model is derived from the complete 3D model briefly evoked above. This has been done by summing the Lagrangian of the different rigid bodies and then projecting the Lagrange equations on the longitudinal axis. It should be noted that one of the following step of this work will include the rail and rolling stocks damage in the optimization problem. This objective will impose the use of a complete multibody dynamic model that will be carried out. This is the reason why the complete multibody model is not described in this paper.

To obtain a model whose calculation time is relatively small, several hypotheses are now introduced to simplify the equations. First, it is assumed that the rigid bodies keep the same distance to each other. As a consequence, the internal interactions can be neglected, and the train can be considered as one element evolving at a longitudinal speed denoted by v(t). In particular, the angular speed is assumed to be negligible for the cars and bogies (only the wheelsets turn significantly), which considerably simplifies the expression of the total kinetic energy of the train. In practice, it appears that the rotation term represents approximately 4% of the longitudinal part. Thus, the kinetic energy T (v(t)) is often written, with a corrective term k corr = 1.04, as

T (v(t)) = 1 2 m tot k corr v(t) 2 , (1)
with m tot the total mass of the train. It may vary between two journeys until more or less 5% of a central value (because of different numbers of passengers in the train). As the total mass will be measured directly on in-line trains in the future, we assume it well known. Thus, we determine it for each series thanks to the experimental measurements.

Secondly, we assume there is no slip between wheels and rail, and that the wheel radius is approximately constant and equal for all the wheels. In this case, the traction and braking torques are equally distributed on each wheelset and can be converted in a global longitudinal force applied to the train. In addition, if we assume that the contact parameters have no influence on the driver maximum command force, the longitudinal components of the different braking at time t can be expressed as functions of the command driver u(t) and the train speed v(t),

F T (v(t), u(t)) = u(t) F T max (v(t)) , F B (v(t), u(t)) = u(t) F B max (v(t)) , (2)
where F T max and F B max are respectively the maximal traction and braking capacities. Empirical models based on experiments are chosen for these maximal capacities.

Then, the expression of longitudinal running resistance, which regroups both aerodynamic and contact forces is approximated with the Davis formula [START_REF] Davis | The tractive resistance of electric locomotives and cars[END_REF]:

F Davis (v(t), v w (t)) = A + B v(t) + C (v(t) -v w (t)) |v(t) -v w (t)| , (3)
in which A, B and C are experimental parameters specific to each train. In this expression, coefficient C refers to the aerodynamic air pressure effect, coefficient B mostly comes from the wheel flange friction, and coefficient A refers to the journal friction, the rolling and track resistances. In order to improve the running resistance in curve, a correction term is applied to each car α which is in curve under the form:

F Curve α (s α ) = m α g k Curve R(s α ) , (4)
with s α the longitudinal abscissa, m α the mass, g the acceleration of gravity, and R(s α ) the curve radius of the track associated with the α-th car of the train. As previously, the coefficient k Curve is determined experimentally and is specific to each train. Finally, the weight of the α-th car, projected on the longitudinal axis, is written as

F W α (s α ) = -m α g sin(θ(s α )) , (5)
where θ(s α ) is the declivity of the track.

With these hypotheses and approximations, the dynamics equation in the longitudinal axis can be written as,

m tot k corr dv dt (t) = F T (v(t), u(t)) -F B (v(t), u(t)) -F Davis (v(t), v w (t)) - α F Curve α (s α ) + α F W α (s α ) . (6)
Thanks to this equation, which describes the dynamics of the whole train, we can calculate the position and speed of the train all over the journey, and compute the final speed and final position of the train, which will be included as constraints in the optimization problem.

Energy consumption

From the speed of the train and the driver commands, we can estimate the motor power and the dynamic braking power used by the train. With an appropriate definition of efficiencies, it is then possible to calculate the energy consumed as well as the energy recovered by the system during the braking. Let η T and η R be the traction and recovery efficiencies. Attention is drawn to the fact that the definition of these two parameters is different. As an example, if we say that the traction efficiency is around 87%, this means that the mechanical power of the motor represents 87% of the electric power transmitted by the catenary (the rest is lost in the traction chain). On the opposite, the recovery efficiency is the part of mechanical power that is reinjected in the catenary (around 82%). The efficiency of the traction chain depends on the efficiency of each of its component. The behavior of these components has not been measured by the constructor for the high-speed train under consideration. Consequently, we have chosen to globally model the efficiency by integrating the whole traction chain. We have observed on other types of trains that the efficiency depends on the train speed and on the longitudinal traction force. As we have not access to the values of the efficiency model for the train used as reference, we have proposed an approximated model that directly depends on the mechanical power, and thus indirectly on the train speed and the traction force. We have verified that the result obtained after the identification is close to the efficiency obtained on other types of train and that the simulated energy consumed was close to the measurements. To take into account this dependence, the following simplified models are proposed for the efficiency coefficients,

η T (P T ) = a η P T + b η , (7)
η R (P B ) = c η P B + d η . (8)
In these two equations, P T and P B are respectively the mechanical traction and mechanical braking power, whereas a η , b η , c η and d η are constants that will also be chosen thanks to experimental measurements. However, traction and braking are not the only functions that consume energy. Indeed, a part of the energy provided by the catenary is transmitted to auxiliary tools of the train, such as the air conditioner and the cooling system. To take into account this additional consumption, we define the auxiliary power as a constant P aux . Finally, with motor, dynamic braking, and auxiliary power, the energy consumed by the train can be calculated along the journey.

Parameter identification

In addition to potentially affecting the quality of the model, the different assumptions presented in the former sections resulted in the introduction of a series of parameters that need to be carefully adjusted. The first step of the parameter identification consists in performing a sensitivity analysis in order to determine the most important parameters that will be identified, the other ones being fixed to their initial values provided by SNCF, the French railway company, for the TGV Dasye. For this step, eleven parameters have been considered: the resistant force coefficients, the traction and recovery efficiency constants, the auxiliary power, the corrective parameter for the resistant force in curve, and two corrective parameters related to the efficiency of the traction and of the energy recovery. The sensitivity analysis has thus been performed with respect to each one of this eleven parameters fixed to their initial values. The sensitivity of the energy consumed by the train is then performed using a sensitivity interval of ±10 % around their initial values. The results of this analysis show that the mots important parameters are:

• the resistant force coefficients A, B, and C,

• the traction and recovery efficiency constants a η , b η , c η and d η ,

• the auxiliary power P aux .

From this analysis, we also found that coefficient C and the two coefficients a η and b η associated with the traction efficiency η T played the most important role on the energy consumed during the journey. This is not surprising: the resistant force coefficient C has more impact than A and B because at high speed, the aerodynamic force is preponderant. Moreover, as the traction phases are much more frequent than braking phases, η T has a greater effect on the energy consumed than η R . In the following, these eight influent parameters are called the calibration parameters.

To identify these calibration parameters, measurements were carried out on the French Rhin-Rhône high-speed line (LGV RR), where the train of interest is assumed to run. The characteristics of this line, which are the relative altitude (integration of the declivity θ(s)) centered on the starting altitude of the journey, the curvature R(s), and the speed limitation v max (s) are described in Figure 2. The measurements of the position, the speed, and the energy consumed by the train on different portions of this track are detailed in [START_REF] Bosquet | Modélisation énergétique et identification des trains pour l'écoconception des lignes ferroviaires à grande vitesse[END_REF]. Wind forecasts have also been simulated by Meteo France near the track at each hour. As an illustration, a particular wind field and its interpolation along the track are shown in Figure 1. Based on this information, and keeping in mind that the driver commands were unfortunately not recorded during these measurements, the following procedure is now proposed for the estimation of the eight calibration parameters.

• The initial values of the parameters are fixed to the values provided by the train constructor and operator. • An inverse method is used consisting in calculating the driver commands in such a way that the simulated train trajectory fits with the measured one. In other words, we look for the traction and braking forces so that the simulated and measured speeds and positions are as close as possible. • Finally, the energy consumed is calculated by simulation and is compared to the measured value.

Comparisons between simulated and measured train speed and energy consumed along the track are shown in Figure 5. We can observe that the inverse method is quite accurate as it is possible to recover driver commands that allow to simulate the same trajectory than the measured one (left plot). Nevertheless, we clearly see on the right plot that working on the driver commands only is not sufficient as the associated simulated energy consumed (dashed line) is not equal to its measured value (dotted line). It should be noted that the energy consumed by the train depends on the traction force and on the speed of the train, which both depend on the driver commands. This energy is estimated by using the nonlinear dynamic equation and the energy consumption model presented in Section 2.2.

Values for the calibration parameters can however be found such that speed and energy consumed fit relatively well with the measurements. As the driver commands are not unique (they depend on the parameters), the solution of the identification may not be unique as well. To avoid overfitting, the chosen values for the calibration parameters correspond to a compromise between a good representation of the measured speed and energy consumed, and a proximity to the values provided by the constructor. The result of the identification is shown in Figure 3 (solid line). Here, the simulated and measured energy are very close. More precisely, by playing with the values of the calibration parameters, the mean error (calculated as the L 2 -norm between the measured energy consumed along the track and the simulated one) goes from 6.4 % to 3.0 % for 10 different measured journeys.

As an example, the identified values for the resistant force coefficients obtained by using the above methodology are

A = 2680 N , B = 31.2 N × s/m, C = 0.427 N × s 2 /m 2 .
It should be noted that the parameters-identification method that is proposed can be applied to other types of trains.

Model validation

To quantify the relevance of the identified values, comparisons have also been performed for measurements that have not been used during the identification procedure. As an illustration, Figures 4, and5 compare the energy consumed with initial and identified values for the calibration parameters for data not used during the estimation procedure. Looking at these figures, it can be seen that the predictive capabilities of the model are relatively satisfactory and have been improved compared to the case where the calibration parameters are fixed at their nominal values. More precisely, the mean error goes from 6.9% to 3.1% when optimizing these values. For most series (8 series out of 10), the error is under 3%, which is very encouraging in view of the simplifications introduced, but some measurements have an error slightly over 5%. These differences between simulation and measurements may be due to a shift between speed and energy measurements but can also come from the simplifications listed in the previous section, the presence of specific rolling conditions, such as rain, or uncertainties on the total mass of the train. As an example, running the simulation with a train full of passengers or empty results in variations of the error of around 2%. The auxiliary power, assumed constant, is also a possible source of error. Indeed, it can be dissociated in two parts. The power dedicated to the passenger comfort is approximately constant but the part transmitted to traction or braking auxiliary (cooling system) may vary during the journey. These aspects may be inspected in further studies to improve the quality of the modeling. train during a given journey, while verifying punctuality constraints. For this purpose, we will start by presenting the general constrained optimization problem. A number of simplifications will then be listed in order to allow the identification of interesting solutions in a reasonable computation time. The relevance of this simplified problem will finally be illustrated on two real traffic configurations.

Definition of the constrained optimization problem

Definition of the cost function. In this paper, we focus on the impact of the driver command {u(t), t ≥ 0} on the electric energy consumed by the train during a whole journey. To calculate this energy, we denote by P E the electric power, and we distinguish different configurations depending on the sign of u. On the one hand, if u(t) > 0, the train is in traction phase and the electric power can be written as a function of the traction force F T , the traction efficiency η T , the traction power P T , the auxiliary power P aux , and the train speed v:

P E (v(t), u(t)) = 1 η T (P T ) F T (v(t), u(t)) v(t) + P aux . (9)
On the second hand, if u(t) < 0, the train is braking, and a part of the dynamic braking power can be recovered:

P E (v(t), u(t)) = -η R (P B ) F D (v(t), u(t)) v(t) + P aux , (10)
where η R is the recovery efficiency and F D is the dynamic braking force. By writing t f the total time of the journey, the total energy consumed by the train is finally given by

f E (v, u) = t f 0 P E (v(t), u(t))dt . (11)
In some specific cases, it could be interesting to recover energy at special positions on the track or at specific periods (near cities, rush hour for instance). Thus, it might be interesting to express this energy as an integral in space rather than in time, and to weight certain portions of the track for energy recovery. This recovery criterion is however not considered in this study.

Definition of the constraints. So that the driver commands can be used in real configurations, several constraints have to be added to the optimization problem. The first ones concern the comfort of the passengers, and aim to limit acceleration and jerk. As we only consider the longitudinal problem, these constraints are however not restrictive (they could be more limiting when considering lateral and vertical axes). They are therefore not taken into account in the following of the paper. The second group of constraints relates to circulation safety. In our problem, these constraints mainly concern the speed of the train, which can be limited to take into account a high curve radius, a tunnel, a switch or a well identified state of degradation of the track. Let v max (s) be the speed limitation. The first constraint we will consider is written as

∀s , v(s) ≤ v max (s) , (12)
in which s is the curvilinear abscissa of the track. The last type of considered constraints are the punctuality constraints, which force the train to arrive at the right point x f at the right moment t f with the appropriate speed v f = 0. In this work, we propose to express these constraints as integrals on the train speed v and acceleration a:

t f 0 v(t)dt = x f , (13)
t f 0 a(t)dt = v f = 0 . (14)
Let C be the space of the driver commands allowing to respect the constraints defined by the equations ( 12), ( 13) and [START_REF] Hansen | The CMA evolution strategy: a tutorial[END_REF].

Definition of the constrained optimization problem. If we defined {(x(t; u), v(t; u), a(t; u)), 0 ≤ t ≤ t f } as the solution of the dynamic system associated with the driver command u = {u(t), 0 ≤ t ≤ t f }, the optimization problem is written as,

min u∈C f E (v(u), u) . (15)

Practical solving of the optimization problem

In order to numerically solve the constrained optimization problem defined by Eq. ( 15), we need to proceed to several simplifications.

Constant piecewise approximation.

Looking for a function u solution of the former optimization problem, means trying to solve an extremely difficult optimization problem because it is defined in a space of infinite dimension. In order to reduce the complexity of the problem, we first propose to limit the search to piecewise constant functions. To this end, we decompose the track into segments of the same length, this length being adapted to the available calculation capacity. The greater the number of segments, the better it is possible to adapt the train speed to the track design and constraints, but the higher the numerical cost. For instance, this length could be chosen around δ s = 200 m for journey of about one hundred kilometers.

Solving the dynamic equations. To solve the nonlinear system defined by Eq. ( 6), we use the Runge Kutta scheme at the order 4. The associated time step δ t is chosen equal to 0.5 s according to a convergence analysis. For this value of δ t , this scheme gives a good precision of the results in an acceptable calculation time.

Perturbation method. As the punctuality constraints can only be evaluated at the end of the dynamic calculation, most of the tested trajectories will not respect these constraints. In order to avoid a too high rejection rate, and thus very high computation costs, we propose in this paper to test only new driver commands under the form of a deviation around a trajectory verifying the punctuality constraints. In other word, if u is a driver command allowing to respect the punctuality constraints, it is assumed that there is a good chance that, for ∆u small with respect to u, the command u = u + ∆u will be close to respect the constraints. In that case, if v = v + ∆v and a = a + ∆a are the train speed and train acceleration associated with u, whereas v and a are the train speed and train acceleration associated with u, we have

t f 0 v(t) dt = x f , t f 0 a(t) dt = 0 , (16)
such that the new punctuality constraints are rewritten as,

t f 0 ∆v(t) dt = 0 , t f 0 ∆a(t) dt = 0 . (17)
In practice, if we randomly choose ∆u, even if the constraints have a chance to be almost verified, the chances are very low if not zero that they are completely verified. To be able to better compare two driver commands, we propose a double adaptation of the train trajectory. If the train stops before the arrival point, we stop the calculation to avoid negative speed and the constraint on the final position is not verified. On the opposite, if the train is assumed to stop after the arrival point, we stop the calculation at time t f . In this situation, both constraints on final speed and final position are not respected.

Optimization algorithm. In order not to limit the research to the immediate surroundings of the initial solution without drawing a huge number of driver command perturbations, we need to use the information given by the dynamic responses to guide the choice of the perturbations. This is how the CMA-ES method works, which we choose to solve the optimization problem under constraints. For this algorithm, to start the calculation, we need to define a starting point, a starting amplitude, and several hyperparameters (see [START_REF] Hansen | The CMA evolution strategy: a tutorial[END_REF] for further details about this method), the constraints being implemented as a penalization of the cost function with Lagrangian multipliers. Among these hyperparameters, we can mention the size of the population, the number of draws selected, and the adaptation coefficient. Then, at each iteration, the method follows several steps, which can be summarized as follows.

• Draw a population of points following the direction given by the covariance matrix around a starting point. • Estimate the cost function and the constraints for each generated point.

• Select the best points.

• Update the Lagrangian coefficients associated with the constraints.

• Update the covariance matrix and the starting point in order to find a better direction for the next draws.

For the considered application, this method has allowed us to find very interesting solutions with a relatively low number of iterations.

Presentation of the optimization results

The optimization procedure we have proposed in the former section is run for the optimization of the driver command under punctuality and safety constraints. The values of the cost function and constraints are kept as the algorithm is iterated. The simulation is run over the LGV Rhin-Rhône track to verify the quality of the optimal trajectory compared with the energy measurements. Figure 6 displays the cost and constraints distributions represented with box plots of the population with respect to the number of iterations. The final position and speed constraints are expressed relatively to the expected ones. This figure gives a visual interpretation of the random search strategy. In fact, during the first iterations, we see that the standard deviation is important. This means that the method searches roughly a direction to explore. After few iterations, the energy consumed by the train begins to decrease: the method has found an interesting energy-saving direction. However, the punctuality and speed limitation constraints are not respected yet. Then, around 50 iterations, the Lagrangian penalization takes more importance and a new direction is defined. This has an impact on the energy consumed that increases slightly but also on the constraints that are closer to be respected. This method continues until the standard deviation reaches the threshold defined as a hyperparameter. We choose it in order that the solution obtained at the end of the optimization calculation is converged.

The quality of the optimal solution is then verified by comparing it with measurements. Thus, left Figure 7 shows a comparison of the speed between the measured and optimal trajectories. It can be observed that the optimal trajectory reduces the maximum speed of the train when it is possible but brakes later to stay on time. Indeed, reducing the maximal speed of the train leads to a reduction of the resistant force, thus a reduction of the energy consumption, but the punctuality constraint imposes that the measured and simulated trajectories arrive on the same time. In order to reduce the energy consumption but still to arrive on time, the optimal solution is to decrease the maximal speed of the train and brake later. This is not clearly shown in Figure 7 but the measured speed of the train is around 25 km/h during almost 1 km at the end of the journey whereas the optimal solution brakes at the last moment, this is why both trajectories respect the punctuality constraint. If the punctuality constraint has been tighter, braking later would not have been sufficient and the solution would have been to increase slightly the maximum speed. Right Figure 7 displays the associated energy consumed. It can be seen that the solution is effectively better and consumes around 20.9% less energy. Thus, the method is interesting to reduce the energy consumed by trains. It can also be observed that, during the first part of the journey, the energy consumed soars, which corresponds to the acceleration parts. Then, during the second part, the energy consumed still increases but slightly, which coincides with the constant speed part of the journey. The energy consumed decreases, which means that the train is braking and recovering energy. The trajectory and energy consumed may be more precise by choosing a lower discretization path δ s , which should improve the energy saving. We can observe this spatial discretization path on the optimal result because the solution is not smooth and seems to be defined piecewise. Finally, if we compare the trajectory with the declivity of the track, it can be understood why the speed is not constant once the train reaches its maximal speed. Indeed, in term of energy, it is not efficient to accelerate during a slope if we need to brake when going down. Thus, before the kilometer 60, it can be seen that the train accelerates and takes advantage of its inertia to climb the slope without accelerating more. This is shown in Figure 8 that displays the optimal trajectory of the optimization run on another journey (LGV est) in the left figure and the associated altitude of the track in the right figure. Here, the arrival time is set arbitrarily so that the train needs to reach its maximal speed for most of the journey to be on time. Thus, we can highlight the efficiency of the maximal speed constraint, which is well respected. Despite this, we still observe that before each positive declivity (kilometers 31, 37 or 68), the optimal trajectory have important speed, which confirms the previous observations. The punctuality and the speed limitation constraints are respected as we can see on the trajectory (left plot in Figure 7). In this work, the comfort and security constraints have not been taken into account but the Lagrangian penalization offers the capacity to add as many constraints as we need. The only impact will be on the time calculation of the method, which will need more time to converge to an optimal trajectory.

Conclusion

In this paper, a model has been proposed to describe the high-speed train dynamics among the longitudinal axis. The parameters have been identified and the model has been validated thanks to measurements carried out on a real high-speed line. Then, the optimization problem has been defined and numerically solved. Finally, the obtained results have been analyzed and compared it to in-line train.

The method allows the energy consumed to be reduced. It could be used to define an optimal trajectory on each journey to help drivers following the most energy saving trajectory. This application of the method takes even more importance when it is known that the automatic

  5.25 5.50 5.75 6.00 6.25 6.50 6.75 7

Figure 1 .

 1 Figure 1. Wind profile on the track.

Figure 2 .

 2 Figure 2. Rhin-Rhône high-speed line description.

Figure 3 .

 3 Figure 3. Speed (left figure) and energy consumed (right figure) of the train as a function of its position for the measurements, for the initial and the identified values of the parameters (identification process).

Figure 6 .

 6 Figure 6. Random draws cost and constraints (LGV Rhin-Rhone).

Figure 7 .

 7 Figure 7. Optimized solution (LGV Rhin-Rhone).

Optimization of the driver commands

Based on the simplified models presented in Section 2, the objective of this section is to present an algorithm to optimize the driver commands in order to minimize the energy consumed by the train is a technology being developed today. It will need a proper trajectory definition as the driver will probably no longer be in the train.

To complete the method, we may introduce other objectives. The optimal trajectory obtained as a result does not take into account other important stakes of railway world. As an example, the wear of wheels and rolling stokes might be fasten if all trains follow the optimal trajectory. This multi-objective aspect will be inspected in further works. The uncertainties of the model parameters is also another further development as we saw that they impact the model quality. Taking it into account allows a better modeling of the physical system and thus improves the optimal trajectory we may find.


Johan Tordsson
email: johan.tordsson@elastisys.com

Erik Elmroth
email: erik.elmroth@elastisys.com

Odej Kao
email: odej.kao@tu-berlin.de

Li Wu
email: li.wu@elastisys.com

Causal Inference Techniques

Keywords: Performance diagnosis, Microservices, Causal inference, Experimental evaluation, Self-healing

  L'archive ouverte pluridisciplinaire

I. INTRODUCTION

Microservice architecture is a popular paradigm for designing large-scale applications because of its benefit of accelerated delivery time, resilience to failures, and flexible development and deployment [START_REF] Balalaie | Microservices architecture enables devops: Migration to a cloud-native architecture[END_REF]. A microservice-based application is composed of loosely-coupled services, which can be developed with heterogeneous technology stacks and communicate through lightweight protocols [START_REF] Newman | Building microservices: designing fine-grained systems[END_REF]. However, performance issues are inevitable in mcroservice due to its largescale services, complex dependencies, and frequent updates.

With the aid of monitoring tools and anomaly detection in microservices, a large number of measurements can be observed and unexpected behaviors can be detected. However, automatic diagnosis of root causes from the observed data for is difficult for the following challenges. 1) A large volume of anomalous metrics: a performance issue is likely to result in a large number of components emitting anomalous metrics, as the anomaly tends to propagate horizontally across inter-dependent services or resources as well as vertically across multiple layers. To illustrate the phenomena, Figure 1 shows that metrics from multiple data resources (microservice (MS), container, and operating system (OS)) become abnormal when a service exhausts the allocated CPU resource. To complicate further, the number of services can be very large in microservices (e.g., Uber reports 4000 services deployed [START_REF] Lin | Microscope: Pinpoint performance issues with causal graphs in micro-service environments[END_REF]), amplifying the volume of anomalous metrics. 2) Heterogeneous anomaly symptoms: a performance issue tends to manifest heterogeneous anomaly patterns, due to different types of metrics, frequently updated and polyglot services, and the control mechanisms in the system (e.g., load-balancing, auto-scaling, restart, etc). As shown in Figure 1, a CPU hog issue has different anomaly patterns in the metrics of duration, request rate, and resources. 3) A wide range of root causes: root causes of a service performance degradation can be external (e.g., hardware failures) or internal issues (e.g., software bugs), and the causes keep evolving with the frequent updates of microservices. Therefore, historical failure data is unable to contain instances of all root causes, failing the historical data based diagnosis methods. A variety of approaches have been proposed in the literature to locate performance issues in microservice systems (briefly discussed in Section V). One of the dominant approaches is to formulate performance diagnosis as a causal inference problem, which constructs an anomaly propagation graph based on causal inference for locating the root causes. For example, Microscope [START_REF] Lin | Microscope: Pinpoint performance issues with causal graphs in micro-service environments[END_REF] applies the PC algorithm (named after its authors, Peter and Clark) and Loud [START_REF] Mariani | Localizing faults in cloud systems[END_REF] uses Granger causality to obtain the graphs. However, only few CI methods are studied in the existing approaches, and their performance against the challenges lying in microservice performance diagnosis is not well understood.

To this end, we select six representative CI methods from three categories including Granger causality, causal network learning algorithms, and structural equation models (SCM), and evaluate their performance on locating root causes against several dimensions of the challenges in microservice performance diagnosis (the background of CI techniques is introduced in Section II). The root cause diagnosis is based on observable metrics collected from a set of fault injection experiments, where three types of anomalies are injected into multiple different services of two widely used microservice benchmarks. For comparison, we developed a non-CI method that leverages the domain-knowledge of service dependencies. To the best of our knowledge, this work is the first such comprehensive study to be reported on evaluating CI techniques for microservice performance diagnosis. The empirical results presented in this work are useful to understand how well CI methods perform for root cause analysis so as to make a proper decision.

To sum up, our contributions are following threefold: 1) A method for evaluating the performance of CI techniques on microservice performance diagnosis. In this method, we formulate three research questions and implement multiple anomaly scenarios to showcase several dimensions of the challenges, including the large number of metrics, heterogeneous anomaly patterns, and a wide range of root causes (Section III). 2) A comprehensive comparison of six representative CI methods on microservice performance diagnosis, including the accuracy of coarse-grained and fine-grained causes localization, computation overhead, and scalability. We also compare CI methods against a non-CI method based on domain-knowledge (Section IV).

II. BACKGROUND

A. Coarse-grained and fine-grained diagnosis

Given a collection of services S and their performance metrics M in a microservice system, we use M (s) to denote metrics for service s, and m (s) i for an individual metric (e.g., duration, container CPU utilization, server filesystem utilization, etc) for service s. Then, a performance diagnosis problem can be formulated as: giving a set of k metrics M = {m

(s1) 1 , m (s1) 2 , m (s2) 1 , ..., m (sn) k
} in a system with n services, once a service performance degradation is detected, how can we pinpoint the root cause m (src) rc that first identifies the anomaly? We define a coarse-grained diagnosis as when metrics are from service-level (e.g., duration) and only faulty service s rc can be identified, and a fine-grained diagnosis as when metrics from multiple data resources (e.g., service, container, and server) and the detailed explanation of service's abnormality m (src) rc can be identified. The main objective of applying CI methods to performance diagnosis is to identify the causal graph with giving (anomalous) metrics, showing the anomaly propagation paths. For example, Figure 2(a) shows four metrics (m 1 -m 4 ) collected from a microservices system, and Figure 2(b) shows a potential causal graph, where m 1 is the source of the anomaly.

B. Overview of causal inference methods

There is an extensive study on inferring causal relations from observational data with CI methods in [START_REF] Nichols | Causal inference with observational data[END_REF]. Here we give a brief introduction of causal inference techniques.

Granger causality: Granger causality (GC) is one of the most common methods for inferring the causal relationship between time-series based on temporal precedence. A time series X = {x 1 , x 2 , ..., x t , ...} is said to Granger-cause another time series Y if including information about the past of X significantly increases the prediction accuracy of the current value y t of Y in comparison to predicting it based on the past values of Y alone. For example, metric m 1 is said to G-cause m 2 in Figure 2(c) if lag 1 is detected between them.

GC [START_REF] Granger | Investigating causal relations by econometric models and cross-spectral methods[END_REF] was initially formalized as a prediction model using vector AutoRegressive (VAR). However, the used regression model assumes linearity and stationarity, and the accuracy of the causal inference is highly affected by predefined model order (which specifies how many previous time points are taken into account in the regression). There are some extensions of the basic GC concept, such as modeling nonlinear dependencies with transfer entropy, including multi-variable in GC, etc. Nevertheless, GC is limited to lagged causal dependencies and has known deficiencies in contemporaneous effects in sub-sampled time series [START_REF] Spirtes | Causal discovery and inference: Concepts and recent methodological advances[END_REF].

Causal network learning algorithms: One of the most advanced theories with widespread recognition in discovering causality is the Causal Bayesian Network (CBN), which uses a Directed Acyclic Graph (DAG) to represent the causal relationships. There are two main approaches for learning a DAG: score-based approaches and constraint-based approaches. GES [START_REF] Chickering | Optimal structure identification with greedy search[END_REF] is one score-based approach. It starts with an empty graph, then greedily searches a Markov equivalent class that maximizes a score function. The PC algorithm [START_REF] Spirtes | Causation, prediction, and search[END_REF] is one of the state-of-the-art constraint-based approaches. It starts with a complete, undirected graph and recursively deletes edges based on conditional independence tests. For example, In Figure 2(d), the edge between metric m 1 and m 4 is deleted as m 4 ⊥ m 1 |m 3 . Once the skeleton is obtained with Learning a DAG from data is highly challenging and complex as the number of possible DAGs is super-exponential to the number of nodes. When applied to a high dimensional dataset, the runtime of the PC algorithm is exponential to the number of nodes in the worst case. But if the true underlying DAG is sparse, which is often a reasonable assumption, the runtime is reduced to polynomial run. There are some other methods, like FCI [START_REF] Spirtes | Causation, prediction, and search[END_REF], which accounts for unobserved direct common causes and can still partially identify which links must be causal. PCMCI [START_REF] Runge | Detecting and quantifying causal associations in large nonlinear time series datasets[END_REF] incorporates time-order as a constraint (causes precede effects) and utilizes a set of causal orientation rules to identify causal directions.

Structural causal model (SCM):

The key idea of SCM is to use the invariance of structural equations without committing to a specific functional form [START_REF] Pearl | Causal inference in statistics: An overview[END_REF]. It assumes that the value of each variable is a deterministic function of its direct causes in the graph and some independent factors, like unmeasured disturbances. A general SCM is defined as Equation 1, where Y is the effect, X is the direct cause of Y, N y is an independent factor, f y is a function which links a variable to its direct causes can be any whatsoever. Y = f y (X, N y ).

(1) Based on the concept of structural causal models, a set of variants are proposed to capture the underlying cause-effect relations. LiNGAM [START_REF] Shimizu | A linear non-gaussian acyclic model for causal discovery[END_REF] assumes a linear function as: Y = βX + N y , where X is independent non-Gaussian noise N y . LiNGAM learns β such that X and Y -βX are independent by applying independent component analysis (ICA). In our example Figure 2(e), the effect metric m 2 can be defined as a function of causal metric m 1 with an Non-Gaussian noise N m2 . ANM [START_REF] Hoyer | Nonlinear causal discovery with additive noise models[END_REF] extends the LiNGAM model and assumes that effect is a function of cause: Y = f y (X) + N y , where X ⊥ ⊥ N y . The function is trained to map between m 2 and m 1 , then it tests whether X and Y -f y (X) are independent using kernel independence tests.

III. METHODOLOGY

In this section, we introduce our experiment design which includes three research questions (Section III-A), CI-based performance diagnosis framework (Section III-B), selected CI methods (Section III-C), and fault injection experiments (Section III-D).

A. Experimental Design

To understand the performance of CI techniques on locating root causes in microservices, we first formulate three research questions, which demonstrate several dimensions of the challenges (i.e, numerous metrics and heterogeneous anomaly patterns) in microservice performance diagnosis; To cover the challenge of a wide range of root causes, we conduct a set of experiments by injecting three types of anomalies into multiple different services. Next, we evaluate the performance of CI methods according to the three research questions based on the metrics data collected from microservice benchmarks. Lastly, we study the computation overhead and scalability of selected CI methods. The details of our experiment design are shown in Table I. The three research questions are formulated as follows:

RQ1: How do CI techniques perform on pinpointing the faulty service that initiates the performance issue (coarsegrained diagnosis)?

To get localize the faulty service that initiates the performance issue, we apply CI methods on response times of all services to construct the anomaly propagation among services. The number of metrics is the same as the number of services, and the heterogeneity of anomaly patterns in response times is relatively low. In order to understand the performance of CI methods against different numbers of metrics in a low-level heterogeneity, we compare the performance of CI methods under three types of feature reduction, including no feature reduction, feature reduction, and ideal feature reduction. Moreover, we compare CI methods to a non-CI method which uses the domain knowledge of service dependency graph.

RQ2: How do CI techniques perform on locating the culprit metric that contributes to the faulty service's abnormality with giving the faulty service (fine-grained diagnosis with giving the coarse-grained cause)?

We assume the faulty service is known, which can be satisfied with our previous work [START_REF] Wu | MicroRCA: Root cause localization of performance issues in microservices[END_REF] and other faulty service localization methods [START_REF] Lin | Microscope: Pinpoint performance issues with causal graphs in micro-service environments[END_REF], [START_REF] Chen | Causeinfer: Automatic and distributed performance diagnosis with hierarchical causality graph in large distributed systems[END_REF], then apply CI methods on metrics exposed by the faulty service to identify the culprit metric that leads to the service's abnormality. The number of metrics is relatively low, but the heterogeneity among them increases as different types of metrics are likely to have different anomaly patterns. We use all observable metrics to pinpoint the culprit metric that initiates the anomaly. The number of metrics is n times larger than the number of metrics in RQ2, where n is the number of services. Meanwhile, the heterogeneity of metrics is in a high level as diverse anomaly patterns exist not only in different types of metrics but also in heterogeneous services. With a high-level heterogeneity, we also compare the performance of CI methods against feature reduction.

B. Performance Diagnosis Framework

Figure 3 shows the framework for performance diagnosis with CI techniques. The input of the framework is metrics corresponding to the research questions (metrics for research questions are listed in Table I) with an optional feature reduction. Next, each CI method is employed to construct an anomaly propagation graph. Finally, it ranks the nodes in the graph and outputs a ranked list of potential root causes, where the top items have the highest probability to be the root causes. Based on the input metrics, the output root causes can be coarse-grained and fine-grained. In our evaluation, we implement feature reduction with an unsupervised machine learning algorithm called Birch clustering [START_REF] Gulenko | Detecting anomalous behavior of black-box services modeled with distance-based online clustering[END_REF] for identifying abnormal metrics, and utilize the PageRank graph centrality algorithm, which is commonly used in the state-of-the-art performance diagnosis approaches, to rank the root causes. Notably, we evaluate CI techniques with a basic ranking method as our focus is on comparing the CI techniques. For the actual application of such techniques, the performance can be improved with advanced ranking methods [START_REF] Wu | MicroRCA: Root cause localization of performance issues in microservices[END_REF], [START_REF] Meng | Localizing failure root causes in a microservice through causality inference[END_REF].

C. Causal Inference Techniques

The objective of this work is to understand the performance of the class of CI techniques on microservice performance diagnosis. To achieve this goal, we select six representative CI methods from three categories, namely Granger causality, causal network learning algorithms, and structural equation models (SCM), based on their popularity, assumptions, and availability of implementations. Notably, we use the pure form of these methods rather than the tailed ones used in existing approaches [START_REF] Chen | Causeinfer: Automatic and distributed performance diagnosis with hierarchical causality graph in large distributed systems[END_REF], [START_REF] Meng | Localizing failure root causes in a microservice through causality inference[END_REF] as we aim to understand the performance of different types of CI methods on microservice performance diagnosis in a general way, rather than to achieve a high accuracy of diagnosis. The chosen CI methods are as follows:

• Granger causality (GC): GC [START_REF] Granger | Investigating causal relations by econometric models and cross-spectral methods[END_REF] is a popular CI method that has been adopted by many performance diagnosis systems [START_REF] Mariani | Localizing faults in cloud systems[END_REF], [START_REF] Thalheim | Sieve: Actionable insights from monitored metrics in distributed systems[END_REF] to build the causal graph. We use χ 2 as the Granger causality test. • PC algorithm with partial correlation (PC-corr): PCcorr is a version based on PC algorithm [START_REF] Kalisch | Estimating high-dimensional directed acyclic graphs with the pc-algorithm[END_REF] that uses Fisher's z-transformation of the partial correlation to test the (conditional) independence. • PC algorithm with kernel (PC-kernel): PC-kernel is a version of PC algorithm that uses kernel-based independence criteria [START_REF] Gretton | Kernel methods for measuring independence[END_REF] which can deal with non-linear causal-effect relationships and non-Gaussian noise. We use the Hilbert-Schmidt independence criterion to test independence. • Greedy Equivalence Search (GES): GES [START_REF] Chickering | Optimal structure identification with greedy search[END_REF] is a scorebased CI method. We use "obs" which is an 0 -penalized Gaussian maximum likelihood estimator based on BIC score in the implementation. • Causal Additive models (CAM): CAM [START_REF] Bühlmann | Cam: Causal additive models, high-dimensional order search and penalized regression[END_REF] is an SCM method which identifies causal relations by fitting an additive SEM with Gaussian error, where the causaleffect relationships can be non-linear. • Linear Non-Gaussian Acyclic Model (LiNGAM):

LiNGAM [START_REF] Shimizu | A linear non-gaussian acyclic model for causal discovery[END_REF] is an SCM method that assumes linear causal relationships but with non-Gaussian disturbance.

We conclude the chosen CI methods and their hyperparameters in Table II. In addition, we implement a non-CI method named Corr for comparisons, which uses the service dependency graph directly and incorporates the anomaly symptoms by setting weights with correlation coefficients.

D. Benchmarks, Injected Faults, and Evaluation Metrics

The evaluation is based on metrics data collected from fault injection experiments on microservice benchmarks.

Benchmarks: We select two representative microservice benchmarks named sock-shop 1 and train-ticket 2 (one of the largest microservice benchmarks), which simulate an ecommerce website for selling socks and a train ticket booking system, respectively. Sock-shop consists of thirteen microservices and train-ticket has forty-one microservices. Both of them are polyglot (e.g., Java, golang, Node.js, etc) and intercommunicate using REST over HTTP. Comparing to sockshop, train-ticket has longer propagation paths and exposes a larger number of metrics.

Injected Faults: We inject three types of faults that are commonly used in the evaluation of the state-of-the-art performance diagnosis approaches [START_REF] Lin | Microscope: Pinpoint performance issues with causal graphs in micro-service environments[END_REF], [START_REF] Mariani | Localizing faults in cloud systems[END_REF] to different microservices in the benchmarks: (1) Latency, we use a traffic control tool named tc 3 to delay the network packets; In our experiment, we inject 200 milliseconds delay to each microservice; (2) CPU hog, we use stress-ng 4 , a tool to load and stress compute system, to exhaust CPU resources. For non-compute intensive microservices, like payment in sock-shop, we exhaust its CPU heavily with 99% usage and 95% for other microservices; (3) Memory leak: we use stress-ng to allocate memory continuously. For memory-intensive microservices, like carts and orders in sock-shop, a high memory usage introduces a high CPU usage, therefore we only provision 1VM to exhaust memory resource, and 2 VMs for other microservices.

In order to inject the performance issues to microservices, we update the existing Docker images by deploying the above faults injection tools and trigger the anomalies by running the corresponding commands inside the container. Each anomaly lasts 1 minute and the system has 5 minutes to cold down before another injection. To increase the generality, we repeat the injection process 3-5 times for each fault and service.

Evaluation Metrics: We quantify the performance of CI methods on diagnosing the root cause with two metrics: precision at top k (PR@k) and average precision at top k (AP@k). These two evaluation metrics are defined as follows:

1 Sock-shop -https://microservices-demo.github.io/ 2 Train-ticket -https://github.com/FudanSELab/train-ticket 3 tc -https://linux.die.net/man/8/tc 4 stress-ng -https://kernel.ubuntu.com/ cking/stress-ng/

• Precision at top k denotes the probability that the top k results given by an algorithm include the real root cause, denoted as P R@k. A higher P R@k score, especially for small values of k, represents the algorithm correctly identifies the root cause. Let R[i] be the rank of each cause and v rc be the set of root causes. We now define P R@k for a set of given anomalies A as:

P R@k = 1 |A| a∈A i<k (R[i] ∈ v rc ) (min(k, |v rc |)) (2)
• Average Precision at k (AP@k) quantifies the overall performance of an algorithm with an average of PR@k:

AP @k = 1 k 1≤j≤k P R@j. (3)

IV. EXPERIMENT RESULTS

In this section, we present the fault injection experiments setup, the results of the three research questions, and the discussion about the computation overhead and scalability of the selected CI techniques.

A. Testbed

We create a testbed using Kubernetes where we deploy a microservice benchmark and a set of monitoring tools, as shown in Figure 4. In our Kubernetes cluster, there is one master node and multiple worker nodes (4 for sock-shop and 6 for train-ticket), including one worker node is dedicated for data collection and the rest for the microservices. Each worker node in the cluster consists of 4 vCPU and 15 GB memory. In addition, one virtual machine (6 vCPU and 12 GB memory) outside the cluster runs the load generator. In our deployment, we limit the CPU resource to The results of PR@1, PR@3, and AP@5 for coarse-grained diagnosis against feature selections.

and set the replication factor to 1 for the microservice where faults are injected.

Load Generator: We develop a closed-loop load generator using Locust5 , a distributed, open-source load testing tool that simulates concurrent users in an application for each benchmark. We customize the workload to reflect the behavior of real users. For example, in sock-shop, more requests are sent to the entry points front-end and catalogue, and fewer to service shopping, carts, user and orders. In total, we provision 450-600 queries per second to both benchmarks.

Data Collection: We collect metrics from multiple layers of the testbed with a set of monitoring tools, including service mesh Istio6 for microservices, Cadvisor 7 for containers, and node-exporter8 for server nodes. Meanwhile, Prometheus9 is deployed to collect metrics from above tools and is configured to scrape metrics every 5 seconds, which is the highest monitoring frequency feasible in our testbed with the existing versions of the used tools.

B. RQ1: How do CI techniques perform on coarse-grained diagnosis?

We evaluate the performance of CI techniques on locating the faulty service that initiates the performance issues based on the response times of all services. We first take all the metrics without any feature reduction as the input of CI techniques, and obtain the diagnosis result in terms of PR@1, PR@3, and AP@5 shown in Figure 5(a). Overall, all the methods cannot identify the faulty services well (with a maximum of 27% in PR@1). This is because services are inter-dependent not only in abnormal status but also normal status. These normal dependencies in metrics introduce extra causal links in the graph, which renders the graph centrality algorithm unable to get the precise root cause.

Next, we apply a feature reduction which results in an F1-score of 0.6 in anomaly detection. The diagnosis result shown in Figure 5(b) shows that CI technique with feature reduction improves on average 77.4% in PR@3, outperforming this scenario with no feature reduction by 31.6%. In addition, we can see that constraint-based methods achieve a higher performance than non-CI method Corr, like PC-corr achieves an improvement of 19% in terms of PR@1 over Corr. This is because feature reduction with an F1-score of 0.6 includes false positives, resulting spurious propagation paths in Corr. In contrast, CI methods can eliminate parts of the spurious paths, thus increasing the accuracy.

Finally, we evaluate the performance of CI methods with an ideal feature reduction, where only the expected anomalies, including the faulty service and its upstream services, are used to construct the anomaly propagation graph. The performance of all methods is shown in Figure 5(c). Overall, all the methods achieve a higher performance when only true-positives are detected as anomalies, on average, achieving 70.7% in PR@1, improving 39% over feature selection with an F1-score of 0.6. Particularly, both Corr and PC-kernel achieves 100% precision.

We note that the results reported here can be improved if advance methods are applied according to the state-of-the-art methods, like using metrics from multiple layers to score the abnormality of services and assigning attributes to edges and nodes before ranking that have been done in our previous work MicroRCA [START_REF] Wu | MicroRCA: Root cause localization of performance issues in microservices[END_REF].

Summary: Neither CI methods nor the non-CI method based on service dependency graphs can identify the faulty services well without feature reduction, as causal relations exist among both normal and abnormal metrics. However, their performance can be improved if feature selection can be applied before causal inference. In particular, PC-kernel and the non-CI method Corr achieve 100% in precision when only true-positives are detected as anomalies. More importantly, constraint-based CI methods, like PC-corr, are more robust than knowledge-based non-CI method to false positives.

C. RQ2: How do CI techniques perform on fine-grained diagnosis with giving the faulty service?

To address this research question, we apply CI methods to the faulty service's relevant metrics and evaluate their performance on locating the culprit metrics that contributes the service abnormality. Table III shows the performance of each CI method on locating the culprit metric on different types of anomalies in terms of PR@1, PR@3, and AP@5 when no feature reduction is applied. We can see that Granger causality (GC) performs better than other methods on latency issues. This is because latency anomalies propagate to other resources through service performance degradation, which can be reflected as time-lags among metrics, which satisfies the assumption of GC.

Regarding CPU hog issues, most of the CI methods identify the root cause with a higher performance except the CAM method. This is because CPU issues reflect themselves as strong correlations between cause and effect metrics and weak correlations with non-causal metrics. Figure 6 shows the maximum and minimum coefficient of determination r 2 of linear regression (r_max, r_min) and Pearson correlation coefficients (p_max, p_min) between the culprit metric and other metrics. We can see that there is a significant difference between the maximum and minimum values (95.7% of p_max and 12% of p_min). Due to these obvious differences between causal and non-causal metrics, it is easier to diagnose the root cause. Meanwhile, due to this linearity between cause and effect metrics in CPU hog issues, the CAM method performs poorly for this type of anomalies as it models nonlinear causal relations. Conversely, we can see that latency and memory leak issues yield smaller differences between causal and non-causal metrics, resulting in lower performance for those methods that assume linearity.

For memory leak issues, we can see that LiNGAM and GES perform better than other methods in terms of PR@3. This is because memory leak issues manifest themselves in multiple resource metrics (e.g., memory utilization and CPU utilization), which introduces contemporaneous effects that can be difficult to identify for other CI methods. LiNGAM is a model that can capture such effects, thus it performs better for memory leak issues.

Summary: Anomaly symptoms in service relevant metrics vary with the type of anomaly, which violate or satisfy the assumptions of a specific CI method. To exemplify, Granger causality performs well on latency issues, SCM method LiNGAM performs well on instantaneous effects in memory leak issues, and most of the CI methods, particularly the PC-algorithm, perform well for correlated anomaly symptoms caused by CPU hogging. On average, PC-algorithm has a better performance on the three types of anomalies.

D. RQ3: How do CI techniques perform on fine-grained diagnosis with all observational metrics?

We address this research question by evaluating the performance of CI on identifying the culprit metrics from all observed metrics. Figure 7(a) shows the performance of culprit metric diagnosis of each CI method in terms of PR@1, PR@5, PR@10, and AP@15, where no feature reduction is applied. In this experiment, we use a larger k = {1, 5, 10, 15} in the evaluation metrics comparing to RQ1 and RQ2 (k = {1, 3, 5}). This is because the number of potential root causes is much larger, and a small k cannot differentiate the performance of CI methods well. We note that LiNGAM is not included when the number of metrics exceeds the sample size, for which the method is not defined.

G C PC -c or r PC -k er ne l G ES

C A M 0.0 0.2 0.4 0.6 0.8 1.0 PR@1 PR@5 PR@10 AP@15 The results of PR@1, PR@5, PR@10, and AP@15 for fine-grained diagnosis against feature reductions. We can see that most of the methods cannot locate the culprit metric well, with an average of 5.6% in PR@1, and 20.6% in PR@5. This is because (1) higher-dimensional features exaggerate the difficulty to infer a precise anomaly propagation graph; (2) anomalies propagate across services and also their metrics, resulting in many instantaneous effects that are difficult to discover; (3) due to the heterogeneity of services, their anomalous metrics commonly manifest diverse symptoms. It is hard to capture various causal relations with one CI method; (4) normal patterns of metrics also introduce unexpected inference.

Next, we evaluate the performance of CI methods with the feature reduction. Figure 7(b) shows the performance of CI methods on reduced metrics. Overall, all the CI methods achieve an improvement after the feature reduction with achieving an average 5.4% in PR@1 and 29.1% in PR@5 more than no feature selection. Particularly, GC, PC-corr, PC-kernel, and LiNGAM achieve an average of 57% in terms of PR@5. With the feature reduction, some unexpected and uncorrelated patterns can be removed from metrics, thus reducing the inference to causal graph learning. Additionaly, we notice that LiNGAM achieves the highest rank in PR@5 with 60%, which indicates that LiNGAM has a higher performance in identifying fine-grained causes in microservice systems.

Lastly, we compare the performance of CI methods on different groups of metrics corresponding to our three research questions. Table IV shows the performance in terms of PR@3 of each CI method. We can see that the performance of CI methods on dealing with metrics of all services (RQ3) is lower than others, as metrics in RQ3 include diverse anomaly patterns not only from service-relevant metrics but also from heterogeneous services. In comparison, most CI methods achieve a higher performance in identifying the culprit metric from metrics exposed by faulty service.

Summary: It is difficult to identify fine-grained causes from all metrics, which is a mixture of many normal and abnormal metrics (20.6% in PR@5). The performance of CI methods can be improved with the feature reduction, but it does not rank the root cause in the top 5 of the list well (49.7% in PR@5). Instead, a drill-down approach that identifies the faulty service first, then pinpoints the culprit metric that attributes to the faulty service, is more promising to obtain the fine-grained causes.

E. Computation overhead and scalability

Lastly, we discuss the overhead of the selected CI methods and their scalability to the number of microservices and metrics, which are important to an online and scalable microservice performance diagnosis.

Overhead: We analyze the computation overhead of each CI method, varying with the number of metrics it handles. Figure 8(a) shows the computation time of the six evaluated CI methods (running on the same Intel Core i7 version server with 8-core CPU and 16 GB memory) against the number of metrics (from 10 to 30, and with 5 minutes of data gathered from each metric, in total of 60 data points). We can see that the computation time of CAM and PC-kernel increase dramatically with the number of features, and GC also becomes slightly slower when the number of metrics increases.

Scalability: We evaluate the scalability of CI methods on the one the largest benchmark, named train-ticket, which has longer propagation paths and a larger number of metrics.

We conduct the same fault injection experiments on the service station in benchmark train-ticket, which has the longest propagation paths. Then we apply four computation efficient CI methods, including GC, PC-corr, GES, and LiNGAM on metrics corresponding to RQ1 and RQ3.

We evaluate the performance of each CI method by measuring the rank of root cause against different number of metrics, which are produced by changing the threshold of the feature reduction. Figure 8(b) and (c) show the average rank of root cause (lower is better) against a different number of service-level metrics and metrics of all services. We can see from Figure 8(b) that the performance of Corr with the service dependency graph by far outperforms all CI methods. This is because the service station has a high in-degree (upstream services), which results in high ranking accuracy in Corr's graph centrality algorithm. Meanwhile, we can see that LiNGAM identifies the culprit service in the top 5 often as the number of metrics increases. However, from Figure 8(c) we can see that all CI methods have increasing ranks of root causes with the increasing number of metrics, which is due to the same issue of diverse anomaly symptoms mentioned in Section IV-D.

Summary: The computation overhead of some CI methods, such as PC-kernel and GES, increases drastically with the number of metrics. When metrics have similar patterns, like service-level metrics, LiNGAM can identify the faulty service well even when the number of microservices is high. However, identifying the root causes from all metrics in the system is difficult due to the high-level heterogeneity.

V. RELATED WORK

In recent years, many approaches have been proposed to diagnose problems in computer network [START_REF] Steinder | A survey of fault localization techniques in computer networks[END_REF], clouds [START_REF] Ibidunmoye | Performance anomaly detection and bottleneck identification[END_REF] and microservices. Overall, these approaches employ machine learning [START_REF] Gan | Seer: Leveraging big data to navigate the complexity of performance debugging in cloud microservices[END_REF], [START_REF] Wu | Performance Diagnosis in Cloud Microservices using Deep Learning[END_REF], pattern recognition [START_REF] Brandón | Graph-based root cause analysis for service-oriented and microservice architectures[END_REF], graph-based [START_REF] Chen | Causeinfer: Automatic and distributed performance diagnosis with hierarchical causality graph in large distributed systems[END_REF] methods with observability from the system, such as logs [START_REF] Zhou | Logsayer: Log pattern-driven cloud component anomaly diagnosis with machine learning[END_REF], requests execution tracing data [START_REF] Zhou | Latent error prediction and fault localization for microservice applications by learning from system trace logs[END_REF], and metrics [START_REF] Mariani | Localizing faults in cloud systems[END_REF], to diagnose the root cause. Here we review the graph-based approaches based on CI techniques as follows.

Graph-based approach can not only identify the root cause but also provide visibility of the issue with a visualized graph. This kind of approach identifies the root cause by constructing a causal graph from observational data, then inferring and ranking the causes based on the graph. Microscope [START_REF] Lin | Microscope: Pinpoint performance issues with causal graphs in micro-service environments[END_REF] and Mi-croRCA [START_REF] Wu | MicroRCA: Root cause localization of performance issues in microservices[END_REF] pinpoint a coarse-grained root cause by modeling the anomaly propagation through service dependency graph and service co-location.Both leverage the ground truth of the service invocation paths and service deployment by parsing network information or metrics from multi-layer in MSA.

In addition, graph-based approaches are widely employed to locate root cause at a fine granularity [START_REF] Mariani | Localizing faults in cloud systems[END_REF], [START_REF] Chen | Causeinfer: Automatic and distributed performance diagnosis with hierarchical causality graph in large distributed systems[END_REF], [START_REF] Meng | Localizing failure root causes in a microservice through causality inference[END_REF], [START_REF] Thalheim | Sieve: Actionable insights from monitored metrics in distributed systems[END_REF], [START_REF] Qiu | A causality mining and knowledge graph based method of root cause diagnosis for performance anomaly in cloud applications[END_REF] by using causal inference methods to construct the causal graph. Sieve [START_REF] Thalheim | Sieve: Actionable insights from monitored metrics in distributed systems[END_REF] and Loud [START_REF] Mariani | Localizing faults in cloud systems[END_REF] construct causal graphs among key performance indicates (KPIs) with Granger causality tests. However, Granger causality tests assume time-lag and linearity or nonlinearity are embedded in the cause-effect metrics, which is quite restrictive for observational metrics in microservice systems. Causeinfer [START_REF] Chen | Causeinfer: Automatic and distributed performance diagnosis with hierarchical causality graph in large distributed systems[END_REF] and [START_REF] Qiu | A causality mining and knowledge graph based method of root cause diagnosis for performance anomaly in cloud applications[END_REF] discover the causal graph with PC algorithm based on conditional independence test. Even though Causeinfer uses a cross-entropy G-squared independence test to overcome the limitations of strict assumptions that other independence tests require, contemporary effects are difficult to identify with CI tests. Moreover, the causal-effect relations are not sparse among metrics which would add more overhead to the PC algorithm. A variant of the PC algorithm which considers the time-order of metrics is used in MicroCause [START_REF] Meng | Localizing failure root causes in a microservice through causality inference[END_REF] to identify a causality graph of metrics. However, it would inherit both the limitations of timelag based method, and the high computation burden of the PC algorithm. To illustrate these issues and compare methods, we in this paper provide a benchmark for different types of CI techniques on identifying the root causes of different types of anomalies. We note that our experiment results can give some insights for improving the efficiency of CI methods on performance diagnosis.

VI. CONCLUSIONS

In this paper, we evaluate the performance of CI techniques for microservice performance diagnosis using six representative CI methods from three categories. Our evaluation based on a set of fault injection experiments shows that 1) CI methods need to incorporate anomaly symptoms through anomaly detection or anomaly score to achieve high performance. In addition, constraint-based CI methods are more robust to false positives in anomaly detection than the domain knowledgebased non-CI method; 2) The heterogeneity in service relevant metrics depends on the type of anomaly. This makes the performance of CI techniques vary with the types of anomalies, among which PC-algorithm performs well for correlated anomaly symptoms. 3) To identify fine-grained root causes that indicates both the faulty service and its abnormality, a drilldown approach that identifies the faulty service first then looks into the detailed information from its relevant metrics, is better than inferring from all metrics. 4) With an increasing number of services or metrics, the runtime of CI method GES and PCkernel tends to increase drastically, and the precision of finegrained diagnosis decreases. However, CI method LiNGAM can still perform well in identifying the faulty service that initiates the performance issue.

Based on our findings, we propose the following research directions to improve the efficiency of the application of CI techniques for microservice performance diagnosis:

1) Automatic search or online tuning of CI methods: It is hard to cover all types of anomalies with one specific CI method due to the heterogeneous anomaly pattens. Therefore, an automatic selection or online tuning of CI methods based on the characteristics of the observational data can help improve the performance of locating root causes from a wide range.

2) Learning causal graph with soft intervention: The constructed causal graph can be improved with the data collected from soft intervening on microservice systems, such as restart or scale-out. However, before applying any interventions, a risk estimation of the intervention is required.

3) Employing expert knowledge as much as possible: Collaborating CI methods with expert knowledge or ground truth, like service dependency, can eliminate inference from irrelevant or duplicated features, thus increasing the accuracy of causal graph learning and performance diagnosis.

Figure 1 :

 1 Figure 1: Multiple metrics become abnormal due to the anomaly propagation from a CPU hog issue. Moreover, anomaly patterns are different in metrics due to the heterogeneity and dynamics in microservices and also the control mechanisms in the system.

1 Figure 2 :

 12 Figure 2: Overview of causal inference techniques. conditional independence tests, edges are oriented with Dseparation.Learning a DAG from data is highly challenging and complex as the number of possible DAGs is super-exponential to the number of nodes. When applied to a high dimensional dataset, the runtime of the PC algorithm is exponential to the number of nodes in the worst case. But if the true underlying DAG is sparse, which is often a reasonable assumption, the runtime is reduced to polynomial run. There are some other methods, like FCI[START_REF] Spirtes | Causation, prediction, and search[END_REF], which accounts for unobserved direct common causes and can still partially identify which links must be causal. PCMCI[START_REF] Runge | Detecting and quantifying causal associations in large nonlinear time series datasets[END_REF] incorporates time-order as a constraint (causes precede effects) and utilizes a set of causal orientation rules to identify causal directions.Structural causal model (SCM):The key idea of SCM is to use the invariance of structural equations without committing to a specific functional form[START_REF] Pearl | Causal inference in statistics: An overview[END_REF]. It assumes that the value of each variable is a deterministic function of its direct causes in the graph and some independent factors, like unmeasured disturbances. A general SCM is defined as Equation1, where Y is the effect, X is the direct cause of Y, N y is an independent factor, f y is a function which links a variable to its direct causes can be any whatsoever. Y = f y (X, N y ).(1) Based on the concept of structural causal models, a set of variants are proposed to capture the underlying cause-effect relations. LiNGAM[START_REF] Shimizu | A linear non-gaussian acyclic model for causal discovery[END_REF] assumes a linear function as: Y = βX + N y , where X is independent non-Gaussian noise N y . LiNGAM learns β such that X and Y -βX are independent by applying independent component analysis (ICA). In our example Figure2(e), the effect metric m 2 can be defined as a function of causal metric m 1 with an Non-Gaussian noise N m2 . ANM[START_REF] Hoyer | Nonlinear causal discovery with additive noise models[END_REF] extends the LiNGAM model and assumes that effect is a function of cause: Y = f y (X) + N y , where X ⊥ ⊥ N y . The function is trained to map between m 2 and m 1 , then it tests whether X and Y -f y (X) are independent using kernel independence tests.

Figure 3 :

 3 Figure 3: Performance diagnosis framework with CI techniques.

Figure 4 :

 4 Figure 4: An overview of our experimental testbed.

  Figure5: The results of PR@1, PR@3, and AP@5 for coarse-grained diagnosis against feature selections.

Figure 6 :

 6 Figure 6: and Pearson correlation coefficient between culprit metric and other metrics of different types of anomaly.

( a )

 a No feature reduction G C PC -c or r PC -k er ne l G ES C A M Li N G A M

Figure 7 :

 7 Figure 7: The results of PR@1, PR@5, PR@10, and AP@15 for fine-grained diagnosis against feature reductions.

Figure 8 :

 8 Figure 8: Computation overhead, the rank of root causes against the number of metrics or services.

Table I :

 I Research questions, metrics, benchmarks, feature reduction, non-CI comparison, and results

	Research question	Metrics	Benchmark	Feature reduction	Non-CI comparison	Results
	RQ1	response times of all services	sock-shop	N	Y	Section IV-B
				Y
				Ideal
	RQ2	metrics of the faulty service	sock-shop	N	N	Section IV-C
	RQ3	metrics of all services	sock-shop	N	N	Section IV-D
				Y
	Overhead	metrics of RQ3	sock-shop	Y	N	Section IV-E1
	Scalability	metrics of RQ1	train-ticket	Y	Y	Section IV-E2
		metrics of RQ3			N
	RQ3: How do CI techniques perform on locating the culprit
	metric from all available metrics (fine-grained diagnosis)?

Table II :

 II The chosen causal inference methods and their hyperparameters

	Method	Category	Assumptions	Hyperparameters	References
	GC	Time-lag	Linearity	AIC criteria for lag selection;	package statsmodels [22]
				maxlag = 30s; χ 2 test for granger
				causality; p-value = 0.05.
	PC-corr	Constraint-based	Gaussianity; Linearity	CItest = gaussian, alpha=0.01	Causal discovery toolbox [23]
	PC-kernel	Constraint-based	Non-Gaussianity; Non-linearity	CItest = hsic_gamma, alpha=0.01	Causal discovery toolbox [23]
	GES	Score-based	Gaussian; Linearity	score=obs	Causal discovery toolbox [23]
	CAM	SCM	Gaussianity; Non-linearity	score=nonlinear, cutoff=0.001, sel-	Causal discovery toolbox
				method=gamboost, pruning=False,
				prunmethod=gam
	LiNGAM	SCM	Non-Gaussianity; Linearity	nonparametric; max_iter=1000	ICALiNGAM in python
					package lingam [24]

Table III :

 III The performance of CI methods on fine-grained diagnosis with giving faulty service against different type of anomaly.

	Methods	GC	PC-corr	PC-kernel	GES	CAM	LiNGAM	average
	Latency
	PR@1	0.24	0.15	0.15	0.21	0.09	0.21	0.175
	PR@3	0.56	0.41	0.44	0.41	0.38	0.41	0.435
	AP@5	0.55	0.45	0.49	0.44	0.36	0.44	0.455
	CPU Hog
	PR@1	0.23	0.63	0.73	0.17	0.2	0.17	0.355
	PR@3	0.70	0.80	0.83	0.73	0.5	0.73	0.715
	AP@5	0.66	0.85	0.87	0.63	0.45	0.63	0.68
	Memory Leak
	PR@1	0.32	0.11	0.11	0.18	0.21	0.18	0.19
	PR@3	0.57	0.57	0.50	0.61	0.39	0.61	0.54
	AP@5	0.56	0.51	0.49	0.54	0.46	0.54	0.517

Table IV :

 IV The performance of CI methods on different ranges of metrics in terms of PR@3.

			RQ		Metrics		GC	PC-corr	PC-kernel	GES	CAM	LiNGAM	Average
			RQ1	response times of all services	0.47	0.62	0.47	0.25	0.41	0.25	0.45
			RQ2	metrics of culprit service	0.61	0.59	0.59	0.58	0.42	0.58	0.56
			RQ3	metrics of all services with filtering	0.45	0.4	0.43	0.23	0.15	0.26	0.32
	Computation time (s)	10 20 30	GC PC-corr PC-kernel	GES CAM LiNGAM
		0
		10	15 Number of metrics 20	25	30

Locust -https://locust.io/

Istio -https://istio.io/

Cadvior -https://github.com/google/cadvisor

Node-exporter -https://github.com/prometheus/nodeexporter

Prometheus -https://prometheus.io/

ACKNOWLEDGMENT

This work is part of the FogGuru project which has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 765452. The information and views set out in this publication are those of the author(s) and do not necessarily reflect the official opinion of the European Union. Neither the European Union institutions and bodies nor any person acting on their behalf may be held responsible for the use which may be made of the information contained therein.