Mehdi Benallegue

Jean-Paul Laumond

Metastability for High-Dimensional Walking Systems on Stochastically Rough Terrain

Biped walking is a complex task, but usually with a natural limit-cycle behavior when walking on an even ground. However, perturbations during walking can make the robot fall. Several works addressed the issue of measuring the robustness to disturbances, and most methods study the effect of a single perturbation. But when walking, the disturbances can be multiple, such as walking on rough terrains. The metastability is a concept that helps studying the case of multiple disturbances. The performance measure is the expectation of the time during which the robot can keep balance. However, until today, only two methods permit to measure this metrics: the discretization of all the state-space and the Monte-Carlo sampling. The former one cannot address high dimensional state-space and the latter is too much time-consuming when the falls are rare. We propose here a solution for walkers with high-dimensional states, even when falls are very rares. The novelty of this method is to rely on the property that limit-cycle walkers may return to the limit-cycle several times before to fall. This method is even extended to the cases of bifurcations or chaos. We illustrate the performance of the approach with a simulated high-dimensional actuated walking system.

I. INTRODUCTION

Legged walking is a complex, non-linear, and highly constrained task. Nonetheless, the simplifications made by researches in biomechanics, neuroscience and robotics enabled synthetic discoveries on elementary models. One among the most important of these results is the limit-cycle property of stable walking on flat terrain, whether for humans or artificial walking systems. The limit-cycle describes a balanced state of walking, and enables to evaluate some characteristics such as energetic efficiency.

The most prominent models for artificial energy-efficient biped walking are the passive-dynamic robots [START_REF] Mcgeer | Passive dynamic walking[END_REF]. These robots are able to produce stable gaits on shallow slopes without any other source of energy but the gravity. Also, the gait of these walkers have a natural convergence to a limit cycle that depends only on the natural passive dynamics of walking. These models have been extended to flat ground by adding small actuators while keeping most of the passive dynamics and the limit cycle properties [START_REF] Collins | Efficient bipedal robots based on passive-dynamic walkers[END_REF]. These robots are known to have the closest energy-efficiency to humans [START_REF] Mcneill | Walking made simple[END_REF].

However, by making these simplifications and idealizing the environment, the roughness of real grounds is not considered. Nevertheless, the robustness to unexpected disturbances is among the most important differences in walking perfor-mances between humans and robots. For example, robots show very high sensitivity to disturbances due to ground unevenness [START_REF] Hobbelen | Limit cycle walking[END_REF], while humans show impressive capabilities to reject large disturbances. Actively controlled robot have usually specific criteria to respect to keep "balance", such as keeping the center of pressure strictly inside the support convex polygon (see [START_REF] Wieber | On the stability of walking systems[END_REF] for an overview of balance criteria of walking systems). Violating the balance criteria leads generally to an unavoidable fall. On the contrary, for most passivitybased limit-cycle walkers, there is no formal definition of a balance criterion; because the joint positions are not actively controlled, and the natural dynamic passivity leads to an unpredictability on the outcome of the perturbation. Therefore, there is usually no way to decide whether a state is balanced or not, except to wait for a stable gait or a fall.

The improvement of limit cycle walking systems in terms of disturbance rejection requires to be able to quantify and to compare the capabilities of the robots facing perturbations. Several methods were created to evaluate these performances. Some of these methods evaluate the range of perturbations that can be handled by the walker, such as measuring the size of the basin of attraction of the limit-cycle [START_REF] Schwab | Basin of attraction of the simplest walking model[END_REF], or simply quantifying the biggest perturbation that can be absorbed by the walker. Other methods evaluate the rate of convergence to the limit cycle after small perturbations (the largest Floquet multiplier [START_REF] Mcgeer | Passive dynamic walking[END_REF] or the Gait Sensitivity Norm [START_REF] Hobbelen | A disturbance rejection measure for limit cycle walkers: The gait sensitivity norm[END_REF]).

These metrics cannot show the walker performances against repetitive disturbances such as on rough terrains. That's because the two first ones measure the effect of a single perturbation, and the latter ones are local and do not model the non-linearity of the system when getting away from the vicinity of the limit-cycle.

Byl and Tedrake proposed an alternative method to estimate the balance performance of a walker, using the concept of metastability [START_REF] Byl | Metastable Walking on Stochastically Rough Terrain[END_REF]. It consists at considering that, on rough terrains, the probability of falling tends to 1 as time goes to infinity, for any walking system. Therefore, one good way to evaluate the robustness of a walker is to evaluate the expectation of the number of steps that the robot can achieve before to fall, which is also called the mean first-passage time (MFPT). However, the proposed method to compute this expectation is based on a Markov chain model, and requires a discretization of all the state space. This can be possible for low dimensional robots but becomes computationally unrealistic for complex walking systems. An alternate way to compute the MFPT is to try a Monte-Carlo sampling method, which has the advantage of a high accuracy in the evaluation MFPT and an arbitrary precision in the state description, but requires very long simulations when failures are rare [START_REF] Byl | Metastable Legged-Robot Locomotion[END_REF]. Making failure as rare as possible is the objective of the optimization in the design of walking systems. The Monte-Carlo method is then very restrictive.

Our contribution is to extend the metastability study to high dimensional limit-cycle walking systems, while keeping reasonable simulation time. The novelty of our method relies on the property that limit-cycle walkers may return to the limitcycle several times before to fall 1 . We show that this property performs noticeably well to estimate MFPT, especially for very rare falls.

In the next section, we present the metastability concept and its application to limit-cycle walking systems subject to stochastic perturbations. Then we introduce the algorithm to compute MFPT and we show how it can be extended to bifurcations and chaos. In Section III, we present an example of a high-dimensional robot that has an important variability in its limit-cycle, and we show how we measured MFPT. In SectionIV, we conclude our paper with a discussion about the advantages and the limitations of the approach.

II. METASTABILITY ANALYSIS OF WALKING

A. Limit-cycle walking: a metastable system

The metastability is the property of some stochastic dynamic systems to keep a specific behavior for long periods, but being guaranteed to leave this state after a sufficiently long time [START_REF] Talkner | Discrete dynamics and metastability: Mean first passage times and escape rates[END_REF]. The system reaches then what is considered as a failed state (see Figure 1). This class of systems can be found in several fields: such as chemical reactions, crystalline structures, brain wave patterns, etc. Computing MFPT is a usual way to estimate the stability of such kind of systems [START_REF] Muller | Rates and mean first passage times[END_REF].

A Fig. 1.

The perturbations can take away the system from the metastable state 'A'. If the system escapes its basin of attraction, it may fall in a deeper "failed" basin. When the failed state is absorbing, the system never gets back to the metastable state.

Walking on rough terrain is also a metastable system: the limit-cycle is an attractive balanced state, and the fall is considered as an absorbing failed state. An absorbing state is a state from which the system cannot escape.

In order to study the balance performances of a limitcycle walker, it is usual to take profit from the periodicity of walking, by considering the time-discretization constituted by a Poincaré map. The Poincaré map describes the intersection of the state vector with a manifold while evolving in time. For walking case, the manifold is provided naturally from a cyclic event during the gait, and is generally defined by 'the impact instant of the swing leg on the ground'. Let's note R the set of all possible states at impact. The state dynamics between the impacts can be described with a discrete-time dynamic system:

x k+1 = f (x k , w k ) (1)
Where x k ∈ R is the state vector of the walker at the kth impact; w k represents the perturbations that happen during the k-th step and f is the function that describes the dynamics during the step.

This function associates to each state x k ∈ R and a perturbation w k the return state x k+1 , which belongs to R only when the walker does not fall during the step. The fall is a state where there is a contact between a non support limb and the ground. Let's define then a "fallen state"

×

x that gathers all possible states of a fallen walker. This simplification permits to define the Poincaré manifold as R * = R ∪ × x, and to extend the function f by defining that

∀w, f ( × x, w) = × x (2)
which means that the walker cannot get up once it has fallen, making the fall the absorbing failed state. The walker has an attractive limit cycle, which is often a periodic one-step trajectory which intersects R at a given state x (pronounced "x-breve"). In contrast, in some cases, there are multiple period bifurcations and even chaotic behaviors during the limit-cycle [START_REF] Goswami | A study of the passive gait of a compass-like biped robot: symmetry and chaos[END_REF]. In such cases, the limit cycle intersects R in a set of states instead of one. These cases will be discussed in Subsection II-C. For now, we consider x as the unique state, which corresponds to the fixed point of the system when no perturbations occur (see Figure 2):

f (x, 0) = x (3)
This walking system subject to perturbations can be considered as an infinite state machine with two specific states which are the limit-cycle x state and the absorbing fallen state × x (Figure 3).

B. Computing the mean first passage time

In our context, we define MFPT as follows: starting from a given initial state x 0 , MFPT is the expectation of the number of returns that the walker can make before falling. In other terms it is m

= E(K) where K ∈ N s.t. x K = × x and x K+1 = × x.
The initial state x 0 may take any value, but let's consider it, for now, as equal to the limit-cycle state x 0 = x.

In fact, this definition of MFPT does not provide the actual time of the first passage, but gives the number of steps instead. This can seem misleading, because the steps do not necessarily Fig. 2. A 3D line representing the state of the limit-cycle walker with no perturbations. The intersecting plane is the Poincaré manifold, the thickest ring is the limit cycle and the seperate vertex is the fall-state. When starting from a state in the basin of attraction of the limit-cycle, the state intersection with the Poincaré manifold converges to x. When starting from a position outside this basin of attraction, there is a fall. Fig. 3. The limit-cycle walker subject to perturbations as a state transition system. The cloud represents the set of states that are neither fall nor limitcycle state. last the same time, and it can be possible that, in some cases, few steps last longer than more ones. However, we prefer this definition because it enables to compare different robots. Indeed, a small robot has faster dynamics than a taller one, and it may have shorter MFPT while showing better (scaled) disturbances rejection.

A stochastically rough terrain is a terrain with a texture, for which the unevenness follows a probability distribution; this is considered as a stochastic perturbation for a robot designed to walk on flat surfaces. As stated before, we assume that on a stochastically rough terrain, the probability that the robot falls tends to one as the number of steps goes to infinity:

lim k→∞ P (x k = × x) = 1, (4)
And if the robot does not fall, the probability that the robot comes back, at least once, to the limit-cycle tends to one as the number of steps grows. In other terms for any neighborhood V (x) of x in R:

lim k→∞ P (∃k 0 ∈ N * x k0 ∈ V (x)|x k = × x) = 1 (5)
The proposed algorithm relies on these two statements to consider the walking process as a Bernoulli trial: (i) either the robot falls or (ii) it comes back to the limit cycle. Coming back to the limit cycle can be approximated by checking if the state belongs to a given neighborhood V (x) of x which is small enough to be considered as a good approximation. This trial has a probability of fall p f and a probability of limitcycle return (1 -p f ). The robot stops at any falling trial.

Walking then follows a geometric distribution: a Bernoulli process which stops at the first fail. The expectation of the number of returns is then:

r = 1 -p f p f (6)
The mean number of returns does not give the mean number of steps. But if we know m r the mean number of steps of a return and m f the mean number of steps of a fall, MFPT can be estimated by:

m = r m r + m f (7)
The estimations of p f , m r and m f can be obtained by repeating trials. Each trial starts from x, and stops either at x or at × x. We denote pf , r, mr and mf the estimations of p f , r, m r and m f respectively.

Let's evaluate the number of walking steps which are required in order to obtain a good approximation of this value. Suppose that, in order to have the desired precision, we want to make a ∈ N trials. If we use the Monte-Carlo method, we will have to wait for n M C = a. m steps to obtain the evaluation of m. But if we use our method, we will have to wait only for

n = a mf pf + a mr (1 -pf ) = a m pf (8)
which can be dramatically less than n M C when the falls are rare, and never more than n M C . Note that when p f is equal to zero we have m = ∞. In that case, Monte-Carlo sampling is unsuitable, while our method takes n = a. mr (from ( 8)).

The method efficiency relies on the fact that the robot has a non-negligible probability to come back to the limit-cycle state before to fall. However, statement (5) may be wrong for special kinds of perturbations. For example, if perturbations create a new fixed point x , then the walker never reaches x when starting from x . In that case, the previous developments are still valid, but our method may have less advantages compared to Monte-Carlo sampling because the value of (1 -p f ) may be too small. In fact, our method can be applied with the state x having any value, but the biggest possible return probability (1 -p f ) permits to have the fastest sampling.

So, what happens if we start from another state x 0 which is not the limit-cycle state? Then, the statements ( 4) and ( 5) can still be reasonably considered true. In that case, we still can take profit from them and use both x 0 and x as starting state. That requires then a twofold sampling, the first starting from x 0 and the second starting from x. Beside the values p f , m f and m r , we need, starting from x 0 , to estimate also (i) the probability to fall without reaching the limit-cycle p f,0 , (ii) the mean number of steps needed to fall m f,0 , and (iii) the mean number of steps needed to reach the limit-cycle m r,0 . The new MFPT is then estimated by: m = pf,0 mf,0 + (1 -pf,0 )( mr,0 + r mr + mf )

where pf,0 , mf,0 and mr,0 are the estimates of p f,0 , m f,0 and m r,0 , respectively.

C. Bifurcations and chaos

For some cases of limit-cycle walking systems, there exists multiple period bifurcations which can multiply the duration of a limit cycle. The bifurcations happen for passive walkers on steep slopes [START_REF] Goswami | A study of the passive gait of a compass-like biped robot: symmetry and chaos[END_REF], or for limit-cycle robots on repetitive terrains, such as treadmill belts. For some cases the bifurcation spreads until creating what is called a chaotic behavior [START_REF] Goswami | A study of the passive gait of a compass-like biped robot: symmetry and chaos[END_REF]. In this case the fixed point is not stable anymore. Instead, there is a set of states {x 1 , ..., xl } such that: ∀i ∈ {1, ..., l -1}, f (x i , 0) =x i+1 [START_REF] Muller | Rates and mean first passage times[END_REF] and f (x l , 0) =x 1 [START_REF] Schwab | Basin of attraction of the simplest walking model[END_REF] (see Figure 4(a)). When subject to perturbations, the steps can lead from any state to any other, and for a period long enough eventually fails (see Figure 4(b)).

Bifurcations may be handled by three different methods. Method 1-To consider only one limit state: Let's consider, without loss of generality, that x = x1 . In this case, Property (5) holds, and we can apply the method shown in the previous subsection. However, the return probability (1 -p f ) can be very low for high number of limit-cycle states. Therefore the advantage to Monte-Carlo method is less interesting.

Method 2-The Markov chain model: Another solution is to consider the walking system as a stochastic finite-state machine (a Markov chain model) with an absorbing state × x and the other states x1 , ..., xl . When the walker is not in any of these states, it is considered in transition (see Figure 4

m i = l j=0 t ij (m ij + m j ) (12)
with m 0 = 0 being MFPT for the fall-state. This provides the following single-shot computation:

m = (I -T ) -1     l j=0 t 1j m 1j . . . l j=0 t lj m lj     ( 13
)
where T is T with removed first row and first column. This model does not rely on the bifurcation properties of the limit-cycle of walker. In fact, it holds for any set of states x1 , ..., xl . The method can be, for example, adapted to create partial discretization of the state space, e.g. around a limit-cycle. This would enable to reduce the fall probability p f , but with a loss of accuracy because of the discretization. In that sense, our method constitutes a compromise between the method of Byl and Tedrake [START_REF] Byl | Metastable Walking on Stochastically Rough Terrain[END_REF] (full state-space discretization) and Method 1. represents the approximation of the set of limit-cycle states.

Method 3-The limit kernel method: Method 2 shows several advantages in the case of bifurcations. However, if the number of the limit-cycle states is too large, for example in chaotic behaviors, the required matrix can take a lot of memory (theoretically infinite). We then suggest to make an approximation X of the set of all state-space which is reached by this limit behavior, and to use Method 1 with X as limit state.

We call X the limit kernel. The method creates a loss in the accuracy of MFPT estimation, which is likely an overestimation of its actual value. That is because a non limitcycle state x could lie inside X, but non limit-cycle states are usually less balanced than limit-cycle ones; thus the balance of the state x is overvalued. Note that a too big kernel leads to overestimate MFPT while too small kernel to a higher fallwithout-return probability p f and makes MFPT longer to be estimated. In addition, X must be defined in such a way that we must be able to know whether a given state x is inside X or not.

In the following section we present an example of a limit kernel approximated by a multidimensional ellipsoid using the covariance matrix of the limit-states set. The example does not show a true chaotic behavior (a chaos due to the dynamics of the system itself), but the rounding and the timediscretization of our numerical simulations create a big set of limit-cycle states. Then we show that it is possible to reach high accuracy in MFPT estimation with a relatively loosefitting approximation.

III. EXAMPLE: A 4-DOF 2D ROBOT

A. Robot geometry and actuation

The walking robot model we consider is depicted in Fig 5 .  It moves in the sagittal plane. The robot bodies consist of a rigid upper body with two masses, one around the middle of the segment (torso) and another at the top (head); and two rigid legs with a mass on the segment. Each leg is equipped with a prismatic spring-damper located between the mass and the ground. We call toe, the bottom of the mobile part of the leg. The hip-toe length l p at the rest position of the spring is denoted by l p,0 . The masses distribution and limbs proportions taken for the simulation are anthropometric (Table I).

TABLE I THE PARAMETERS VALUES OF THE SIMULATED ROBOT

The robot is considered to be always walking, that means that there is always at least a contact between a toe and the surface of the ground. If this condition is not satisfied, the robot is considered as fallen. The so-called stance leg is the leg in contact with the ground. The other leg is called the swing leg. When a leg is in a swing phase, its toe comes back instantly to the rest position of the spring, so the length of the swing leg is constant with value l p,0 . We denote by l p the length of the stance leg.

The state vector of the robot is: x = θ, θ, η, φ, α, α, l p , l p where θ is the time-derivative of θ.

The simulations are performed using the dynamic simulator ODE with a time discretization of 1 KHz. A special simulation time sample is set to 0.1 ms at the double support instant in order to increase the physical realism of the impacts. Impacts are considered inelastic and contacts are considered perfect with no slipping.

The toe of the stance leg follows a passive spring-damper behavior. The force applied to the ground has a proportionalderivative (PD) expression:

f t = -K toe,p (l p -l p,0 ) -K toe,d . l p ( 14
)
where K toe,p is the elasticity of the spring and K toe,d is the damping factor. At the instant of raise-up of the stance toe, a velocity controlled impulsion is applied to the ground to give a propulsion to the robot (Figure 6). The force is automatically computed by the dynamic simulator to provide the desired velocity:

f t = h( lp,r ) (15)
where lp,r is the given desired velocity and h is the speedreferenced controller. The force f t is applied during one timestep of simulation. The inter-leg joint is also controlled. The actuator is a pure torque generator controlled by a simple PD that generates a torque toward a reference angle:

τ h = -K hip,p (θ -θ r ) -K hip,d . θ (16)
where K hip,p is the proportional gain, θ r is the reference angle and K hip,d is the derivative gain. In order to preserve the natural dynamics of the legs the values of the gains are small. The trunk is stabilized to keep an upright orientation. Among the several methods for controlling the upper body, we choose the solution that generates a torque between the stance leg and the trunk. By doing so, we leave the passivity of the swing leg creating a more efficient motion [START_REF] Mcgeer | Stability and control of two-dimensional biped walking[END_REF]. The trunk torque is actuated by a pure torque generator, controlled by a PD that brings back the torque to vertical:

τ t = -K trunk,p .α -K trunk,d . α (17)
where K trunk,p is the proportional gain and K trunk,d is the derivative gain.

The values of all these parameters appear in Table II.

B. The limit kernel

When starting the simulation on a flat ground from a state in the basin of attraction, the robot reaches a limit behavior within few steps. Figure 7 shows a phase plot of the inclination of the left leg over the time during the limit-cycle. The limit cycle seems to be unique, but there is a variability in the state vectors x k when meeting the Poincaré manifold: the thickness of the curve is the overlay of many curves. This variability is probably due to the time-discretization and rounding errors. This means that there is no real asymptotically stable fixedpoint of the system, and a we have very large set {x 1 , ..., xl }.

Left leg angle (rad) Left leg angular velocity (rad/s) Fig. 7. Phase plot of the limit-cycle of the left leg inclination during several walking steps, each cycle here has a duration of 2 steps.

Our idea is to study the distribution of these states in the multidimensional space. We made a rough sampling of the state space (8000 initial states) and we ran a simulation from each of them. The simulations for the that lie into the basin of attraction reached the limit cycle. In order to be sure to catch most of the possible dynamics of the limit-states, we recorded the states for each successful simulation, between the 400th and the 500th step, obtaining 3000 values of xi .

To analyze the obtained values of the limit cycle, we have used principal components analysis (PCA). In Figure 8, we see that the three first vectors carry the most important part of the variability. Within the directions of these vectors, the distribution is regular and reaches its maximum around its mean value. This suggests that a multidimensional ellipsoid can fit well the points cloud. An ellipsoid is an easy shape to define: a state is inside the ellipsoid if and only if its distance to the center is under a threshold d 0 . What we mean by distance here is an anisotropic metrics defined by a semi positive definite matrix. We propose to use the distance defined by the inverse of the covariance matrix of the set of states (cov({x 1 , ..., xl })) -1 . However, the covariance matrix can be non-invertible. For example, the value φ is always at zero during double support, creating a zero column/row in cov({x 1 , ..., xl }). So, the solution is to use the Moore-Penrose pseudo-inverse of the covariance matrix C = (cov({x 1 , ..., xl }))

+ . The function of distance to the center is then defined by:

d(x) = (x -x) t C(x -x) ( 18
)
where x is the mean value of {x 1 , ..., xl }. Note that even if the pseudo-inverse solves the problem of non-invertibility, this new distance ignores completely the null-space of the covariance matrix. That means that whatever the value of φ, the distance is constant. As a consequence, a special treatment must help to decide whether a state is in the Kernel or not. We choose to follow the same choice as Byl and Tedrake [START_REF] Byl | Metastable Walking on Stochastically Rough Terrain[END_REF]: we ignore φ in the MFPT estimation. This makes our limit kernel approximation oversized.

Figure 9 shows the distribution of the values of d(x i ). Most of the distances to the center are less than 20 while all of them are less than 40. In order to show that a relatively loose-fitting shape can be enough to estimate correctly MFPT, we take the threshold d 0 = 1000.

C. Simulation results

The model of perturbations we propose is to change at each step the slope of the ground (see Figure 10). The slopes follow a Gaussian white noise. The variance of the Gaussian defines the texture of the ground. This kind of perturbations has the advantage to not require any scaling, regarding the size of the robot, because this perturbation is equivalent to a step down or a step up with a height which is proportional to the walkingstep length. We ran the algorithm described in Section II-C with the ellipsod as a limit kernel and starting with x 0 = x and a = 100000 trials were performed for each slope variance. The results are shown in Figure 11(a). The robot has infinite estimated MFPT until 0.05rad of standard deviation for the ground slope. We know that this estimation is inaccurate, but it suggests that the true MFPT is much higher than the biggest MFPT we can evaluate with simulation techniques. We qualify this value as being virtually infinite. While the variance increases, MFPT decreases rapidly to reach less than 10 steps.

More importantly, Figure 11(b) compares the number of steps required by Monte-Carlo sampling and by our method. Considering that we simulate 100 steps in 4 seconds, the total time required to get these results when using our method is of 48 hours, while using Monte-Carlo sampling takes 319 hours. Such time costs have been estimated without considering the 0.05 rad case. including this case is out of the scope of Monte-Carlo method, while it takes only 1.3 hour for our method to estimate the virtually infinite MFPT.

In order to evaluate the accuracy of our method, we also ran a classical Monte-Carlo simulation for the case 0.1rad. MFPT estimation was 22.66 steps, while it was 22.75 steps with our method. This shows an overestimation of 0.37% only but with 50% less simulation time.

We have seen through this example a relatively loose fitting limit kernel which enabled an accurate and a much faster estimation of MFPTs, especially for rare fall cases. Beside the gain over Monte-Carlo, our method has several advantages compared to other measures of robustness to disturbances, as discussed below. MFPT and the number of steps required for the estimation. Both are in logarithmic scale.

IV. DISCUSSION

This section presents the benefits of our method over the existing approaches as well as its limitations.

In contrast to the previous robustness estimators, the metastability based estimator enables to analyze the effect of repetitive disturbances for which Property (5) holds. Indeed, the Floquet multiplier and the gait sensitivity norm methods study the convergence to the limit-cycle after a small perturbation. These approaches can be extended to the presence of repetitive disturbances but they remain local around the limitcycle. They cannot show under which conditions the robot would fall. Furthermore these methods become less suitable in the case of bifurcations, and unable to address the case of chaos.

Some other methods, such as searching the largest admissible perturbation and measuring the size of the basin of attraction, are more global and do not need the presence of a stable limit-cycle. However, they cannot be adapted to repetitive disturbances.

MFPT is then a good metrics to measure the tolerance to rough terrains. The method by Byl and Tedrake [START_REF] Byl | Metastable Walking on Stochastically Rough Terrain[END_REF] to estimate MFPT needs a discretization of all the state-space. Let's take the example of 8-dimensional state, which may be reduced to 7-dimensions if we ignore φ. Suppose that we want to discretize each dimension with 20 values, which is very rough, the discretization would take 20 7 = 1.28 billion of states. Again, the computation requires the inversion of a square matrix of 20 7 × 20 7 = 1.6384 × 10 18 cells. Our method does not have any scaling problem with the increase of the complexity of the walking system. Furthermore, our method permits an arbitrary accuracy in MFPT estimation. That's because the only states which are approximated are the limit-cycle states, and we can choose any tight approximation, with, in return, a loss in the performances of the method, but still being better than Monte-Carlo method.

However, as acknowledged, our method is not relevant if the limit-cycle changes because of the perturbations. To overcome this drawback, a solution can be to identify the new limitcycle and to change x for the new one. Also our method fails when important chaos cannot be approximated correctly enough. Nevertheless, we keep in mind that our method can never be worse than a Monte-Carlo approach that remains the only alternative.

Finally, the biggest disadvantage of our method is that it is more time-consuming than approaches based on Floquet multipliers and gait sensitivity norm. But there is virtually no other solution for such nonlinear system than trying several times or multiple possibilities in order to have a synthetic, global estimator about walking performances on textured terrains.

V. CONCLUSION

This paper has shown how to take advantage of the existence of the limit cycle to efficiently estimate metastability properties of a walking robot. The method considers the walker as a Markov chain, with a special absorbing state representing the fall, and other states representing the limit-cycle behavior. When the walker is neither fallen nor in a limit-cycle state, it is considered in transition. With this model, the computation of MFPT becomes straightforward.

On the other hand, the concept of metastability is particularly adapted to the model of stochastically rough terrain. The given example have shown that the assumption which is the basis of our approach is nicely verified with this model of perturbations: the robot comes back regularly to the limit cycle, and this may happen even for high variances of the terrain inclinations.

Finally, MFPT is a good criterion to measure the robustness of walking systems to repetitive disturbances. It can be used either (i) to estimate the tolerance of a walking system to given disturbances and to give an estimation of the distance that can be traveled before falling, (ii) to use it for comparing different robots, even if they have different dynamics, (iii) to optimize the dynamical structure of a robot in order to maximize MFPT, or (iv) to compare different control policies for a same robot.

  (c)). Compared to the Markov chain model presented by Byl and Tedrake [3], one transition does not represent one walking step but potentially several ones. In this case, we need to sample differently: starting from each states xi and xj we estimate (i) the transition distribution (the probability t ij = P (x k+1 = xj |x k = xi )) and (ii) the mean number of steps m ij required to go from xi to xj . For notation simplification, we denote the fall-state × x by x0 . This gives the the matrix M = (m ij ) and the probability transition matrix T = (t ij ) whose first line is 1 0 • • • 0 because of the absorbing fall-state. Let's define m = m 1 , • • • , m l t being the vector of all MFPTs of the limit-cycle states x1 , ..., xl . Vector m is easy to compute, by considering the relationship:

  Fig. 4. A simple example of three-states bifurcation. The cloud represents the set of states that are neither fall nor limit-cycle state. The doubled arrows represent transitions of the Markov chain model, and each transition between xi and xj takes the mean number of steps m ij . The shape labeled ' X' represents the approximation of the set of limit-cycle states.

Fig. 5 .

 5 Fig. 5. A sketch representation of the actuated 4 DOF model

Fig. 6 .

 6 Fig. 6. The actuation and passivity of toes

Ktoe, p

 p = 5000 n/m K toe,d = 200 n.s/m lp,r = 1 m/s K hip,p = 1 n.m/rad K hip,d = 0.15 n.m.s/rad θr = 0.3 rad K trunk,p = 30 n.m/rad K trunk,d = 15 n.m.s/rad

  Fig. 8.The PCA of the limit-cycle states values. We see in (a) that the three first components carry most of the information. In (b), (c) and (d) we see histograms of their distribution.

Fig. 9 .

 9 Fig.9. A histogram of the repatition of d(x i ). There is no state at almost 0 distance because the multidimensional space is very small for small distances.

Fig. 10 .

 10 Fig. 10. The perturbations model. Each step has a different slope inclination. The inclinations follow a Gaussian white noise.

  (a) MFPT estimations for different ground texture. (b) The number of simulated steps required for the estimation.

Fig. 11 .

 11 Fig. 11.MFPT and the number of steps required for the estimation. Both are in logarithmic scale.

TABLE II THE

 II ACTUATION/PASSIVITY PARAMETERS OF THE ROBOT

Section III presents an example where this property appears.

ACKNOWLEDGMENTS

This work is partially supported by grants from the EU Project CLONS (Closed-Loop Neural Prosthesis for Vestibular Disorders, Grant Agreement 225929). We deeply thank Alain Berthoz for his support and our stimulating discussions.