Common formation mechanism of basin of attraction for bipedal walking models by saddle hyperbolicity and hybrid dynamics

In this paper, we investigate the mathematical structures and mechanisms of bipedal walking from a dynamical viewpoint. Especially, we focus on the basin of attraction since it determines the stability of bipedal walking. We treat two similar but different bipedal walking models (passive and active dynamic walking models) and examine common mathematical structure between these models. We find that the saddle hyperbolicity and hybrid system play important roles for the shape of the basin of attraction in both models, which are quite common for more general bipedal models and important for understanding the stability mechanism of bipedal walking.


Introduction
In this paper, we study bipedal walking using the mathematical models. Especially, we focus on the basin of attraction of stable walking. In a bipedal walking model, a limit cycle corresponds to stable walking and the size and shape of the basin of attraction of the limit cycle determine the robustness of walking for various noises and disturbances. The study for the stability of walking will contribute to designing a biped robot and walking support system, as well as to understanding human walking.
The study using a mathematical model often has a question as to whether the result of analysis using the model is essential for the phenomenon or model specific. One method to answer this question is to compare various models and to find a common mechanism, which may suggest an essential feature of the phenomenon beyond mathematical models. In this paper, we compare two models (passive and active dynamic walking models) and search for a common mechanism of formation of the basin of attraction of their stable walking.
The passive dynamic walking was proposed by McGeer [10], which walks down a shallow slope without any actuator or controller. To investigate the linear stability, the simplest walking model was introduced by Garcia et al. [6] and the basin of attraction was computed by Schwab and Wisse [16] (we use this model as a passive dynamic walking model). They showed that the basin of attraction is very small and thin, and it has a fractal-like shape. However, why the basin of attraction has such a shape remains as an open question. In [13], we introduced some new ideas about a mechanism of forming the shape of basin of attraction and showed that the saddle hyperbolicity of the upright equilibrium point plays an important role to form the basin of attraction.
In this paper, in addition to the passive dynamic walking model, we take an active dynamic walking model. Different from the passive dynamic walking model, it walks on even ground using an actuator controlled by a phase oscillator, which is inspired by central pattern generator [2,17]. This paper shows that the basic formation mechanism of basin of attraction of the passive dynamic walking model in [13] is common to that of the active dynamic walking model, despite such differences in their models.
On the study of bipedal walking, the following two facts are important: • Bipedal walking is a typical hybrid system due to foot contact and foot off • The center of mass of the human body moves like an inverted pendulum (the inverted pendulum mechanism [9,14]) Modeling bipedal walking should reflect these facts. More specifically, during human walking, there are two states called the single support phase and the double support phase. On the single support phase, one leg (called stance leg) supports the body and the other leg (called swing leg) swings from back to front. When the swing leg contacts the ground, the state switches to the double support phase, where the both legs supports the body. The former stance leg lifts off the ground and the state switches to the single support phase again. Physical condition between the single and double support phases is quite different and these states are dominated by different equations of motion. This means that the dynamical system of bipedal walking is a hybrid system. When we observe human walking more closely, we find that the stance leg is almost straight, and it rotates around the foot contact point like an inverted pendulum. Therefore, the center of mass is at its highest position during the midstance phase and at its lowest position during the double support phase. In contrast, the locomotion speed is lowest during the midstance phase and highest during the double support phase. This means that humans produce efficient walking through the pendular exchange of potential and kinetic energy while conserving mechanical energy [3][4][5]. This is called the inverted pendulum mechanism [9,14], and inverted pendulums have been widely used as the simplest model for the movement of the center of mass, when investigating the underlying mechanism in human walking [1,7,8,11,12].
In the present study, we aim to clarify the mechanism that determines the geometric characteristics of the basin of attraction for our bipedal walking models by considering the theory of dynamical systems and focusing on the saddle point or the saddle periodic orbit that is inherent in the governing dynamics related to the inverted pendulum. We show that the basin of attraction is quite thin, and that its reason is common between these two models. In [13], we showed that the saddle property (hyperbolicity) of the upright equilibrium point and the λ-lemma, one of the most basic properties of the hyperbolicity of dynamical systems, are important for the thin basin of attraction in the passive dynamic walking model.
In this paper, we show that the saddle-center periodic orbit of the active dynamic walking model plays a similar role. Indeed, we show that the basic formation mechanism of the basin of attraction is common between the two models, although the detailed structure is rather different. We have already shown that Poincaré sections and center-stable/center-unstable manifolds play important roles for such a shape to be formed by the passive dynamic walking model [13]. In this paper, a similar formation mechanism is also found in the active dynamic walking model.
Because the saddle property is embedded in general locomotor systems (that are not limited to our models), our results may contribute not only to elucidating the stability mechanism in our bipedal walking models, but also to improving the understanding of the stability mechanism in human walking and thus to producing design principles for the control of biped robots and walking-support systems.

Model
In this paper, we use a simple compass-type bipedal walking model (Fig. 1). This model has two legs (rigid links), each having the length l, that are connected by a frictionless hip joint. Let the angle of the stance leg with respect to the slope normal be θ 1 , and the angle between the stance leg and the swing leg be θ 2 . The mass is located only at the hip and the feet; the hip mass is M and the foot mass is m. g is the acceleration due to gravity. This model walks on a slope of angle γ and is controlled by the input torque u at the hip.
In this model, the typical walking behavior is as follows. A new step starts when both feet are on the slope, just after the swing leg makes contact with the slope. The The passive dynamic walking model walks on a slope (γ > 0) without any input torque (u = 0), while the active dynamic walking model walks walks on even plane (γ = 0) with input torque (u = 0) front leg is the stance leg, and the other leg is the swing leg. The stance foot is fixed on the slope, and the stance leg rotates freely without friction. The stance and swing legs move as a double pendulum. The swing leg swings forward, and the swing foot contacts the slope. In this model, the collision is assumed to be fully inelastic (no slip, no-bound). The swing leg immediately becomes the new stance leg, and vice versa (double support duration is infinitesimal).
The passive dynamic walking model describes walking down the slope (γ > 0) without any torque (u = 0) by balancing the energy dissipation due to foot contact with the energy generation due to the gravitational potential energy. In contrast, the active dynamic walking model walks on even ground (γ = 0) controlled by the input torque (u = 0).
We note that since the legs are rigid links, the swing leg collides with the slope when the stance leg is nearly vertical. We can avoid this foot scuffing by adding complications to the model, such as passive knees. However, for simplicity, we ignore the foot scuffing and allow the leg to pass through the slope in our models.

Equations of motion for the single support phase
The configuration of the mechanical model is described by two variables (θ 1 , θ 2 ). The equations of motion are given by a Lagrangian equation: For the passive dynamic walking model, we take the torque u = 0. For the active dynamic walking model, we use a phase oscillator for generating a torque, whose phase is φ, to control the model. The oscillator phase follows the dynamics where ω is the frequency and we simply determine the input torque u by In the passive dynamic walking model, the phase space is four dimensional with variables (θ 1 , θ 2 ,θ 1 ,θ 2 ). On the other hand, in the active dynamic walking model the phase space is five dimensional with (θ 1 , θ 2 ,θ 1 ,θ 2 , φ).
After appropriate rescaling, we have the following non-dimensionalized equations.

Foot contact
The swing foot contacts the slope when the following conditions are satisfied: Conditions (3) and (4) are used to ignore the foot scuffing when the swing leg moves forward. We assume that foot contact is a fully inelastic collision (no-slip, no-bound) and that the stance foot lifts off the slope as soon as the swing foot hits the slope. The relationship between the state just before foot contact (θ − 1 , θ − 2 ) and the state just after foot contact (θ + 1 , θ + 2 ) is as follows: Due to the collision, the angular velocities discontinuously change just at the moment of foot contact. We assume that the stance leg does not interact with the slope when it lifts off and that the input torque does not work at the instant. From these assumptions, the conservation of angular momentum yields the following relationships: .
Since the roles of the legs are swapped at the collision so that θ 2 varies as (5), we change the oscillator phase: Note that the state just after foot contact depends only on . From the fact, in the passive dynamic walking model, the state just after foot contact forms a two-dimensional submanifold in the four-dimensional phase space since φ is ignored. In the active dynamic walking model, the state just after foot contact forms a three-dimensional submanifold in the five-dimensional phase space.

Structure of phase space by hybrid dynamics
The models are hybrid systems composed of the continuous dynamics during the single support phase and the discontinuous dynamics at foot contact. The hybrid dynamics determines the structure of the phase space, as shown in Fig. 2a. H is the section of foot contact defined by the conditions (2)-(4). T is the jump in the phase space from the state just before foot contact to the state just after foot contact, defined by the relationship (6). Therefore, the image of T , T (H ), is the region representing all states just after foot contact and a new step starts from T (H ). U is the map from the start of a step to the foot contact. In other words, U is the map from T (H ) to H , defined by the equations of motion (1). The Poincaré map S is defined by S = T •U : T (H ) → T (H ) on the Poincaré section T (H ). This Poincaré map represents one gait cycle, and an attractor of the Poincaré map represents stable walking.
Note that the structures of the phase spaces of the passive and active dynamic walking models are quite similar, but the number of dimensions is different. In the passive dynamic walking model, the phase space is four dimensional, the section is three dimensional, and T (H ) is two dimensional. On the other hand, in the active dynamic walking model, the phase space is five dimensional, the section is four dimensional, and T (H ) is three dimensional. We also consider the sequence of inverse images of D, S −n (D) (n = 1, 2, . . .). The region S −n (D) indicates the collections of initial conditions on which the model takes at least n + 1 steps. The sequence approximates the basin of attraction, and we investigate the mechanism by which the shape of the basin of attraction is formed from the geometric structure of these inverse images.
In the following sections, we numerically compute D, S −1 (D), S −2 (D), . . .. Numerical computation of the region D is rather straightforward. We take many initial points in T (H ), numerically integrate the equations of motion from these initial points and check whether the orbit can take one step or falls down. We use fall down threshold at θ 1 = ±π/2. T (H ) is two dimensional and parameterized by two variable θ 1 ,θ 1 in the passive dynamic dynamic walking model, and T (H ) is three dimensional and parameterized by three variable θ 1 ,θ 1 , φ in the active dynamic dynamic walking model. Therefore we use these variables to compute and show the numerical results. We can compute S −1 (D), S −2 (D), . . . in a similar way. Because S −n (D) for sufficiently large n is considered to be a sharp approximation of the basin of attraction, we compute S −n (D) for n = 50 and n = 200 and take them as the basin of attraction if the two numerical results are same 1 .
In this paper, we choose the following parameters to analyze the models. In the passive dynamic walking model, we use A 0 = 0, β = 0, γ = 0.011, where we consider a limit case in which the foot mass is much smaller than the hip mass as in [6]. With this parameter, the corresponding Poincaré map has a unique attracting fixed point at (θ 1 ,θ 1 ) ≈ (0.214, −0.212) on the Poincaré section. In the active dynamic walking model, we use A 0 = −0.29135, β = 0.15, γ = 0, ω 0 = 1.1191. We use these parameters so that the corresponding Poincaré map has a unique attracting fixed point at (θ 1 ,θ 1 , φ) ≈ (0.144, −0.179, 0.0535) on the Poincaré section. More complicated cases for the passive dynamic walking model (for example, the Poincaré map has a chaotic attractor) were analyzed in the paper [13].

Center-stable and center-unstable manifolds
The equations of motion (1) for the passive dynamic walking model have an equilibrium point (θ 1 ,θ 1 , θ 2 ,θ 2 ) = (γ , 0, 0, 0), where the legs remain upright. The equilibrium point is deeply related to the geometric structure of the basin of attraction, as explained in Sect. 3. The eigenvalues of the linearized equations of motion at the equilibrium point are ±1 and ±i, and the equilibrium point is a saddle-center with one stable direction, one unstable direction, and two neutral directions.
The equations of motion (1) for the active dynamic walking model (at the selected parameters) have a 2π/ω 0 -periodic orbit 2 . We can find this periodic orbit by the Newton method. The time-2π/ω 0 map of the equations of motion from φ = 0 to 2π can be considered as a map from R 4 to itself, which has a saddle-center fixed point at (θ 1 ,θ 1 , θ 2 ,θ 2 ) ∼ (−0.041, 0.000, 0.891, 0.000). The eigenvalues of the linearized matrix of the map are approximately 316.8, 0.003157, 0.4021±0.9156i: one is greater than 1, another is positive and less than 1, and the other two are complex conjugate and on the unit circle, which shows that the periodic orbit is of the saddle-center type.
The equilibrium point for the passive dynamic walking model and the periodic orbit for the active dynamic walking model have a codimension one center-stable manifold (W cs ) and a codimension one center-unstable manifold (W cu ). Fig. 3a, b show these structures.

Basin of attraction in passive dynamic walking model
In this section, we briefly summarize results on the passive dynamic walking model given in [13] for comparison with the active dynamic walking model. Figure 4a shows the domain D and the basin of attraction B on T (H ). Both D and B are very thin in the space of (θ 1 ,θ 1 ). Fig. 4b is a zoom-in view. To clearly see the geometrical details, we use θ 1 +θ 1 and θ 1 −θ 1 for the axis in Fig. 4c. The intersection of the center-stable manifold W cs and T (H ) is shown by a green line in Fig. 4b, c. We showed that D had the following properties: a b Fig. 3 Phase diagram of the passive (a) and active (b) dynamic walking models. The passive dynamic walking model has a saddle-center equilibrium point (open dot in a) and the active dynamic walking model has a saddle-center periodic orbit (black arrow in b). a The center-stable and center-unstable manifolds of the equilibrium point are represented by the green and blue curves, and b the center-stable and centerunstable manifolds of the periodic orbit are represented by the green and blue surfaces. b Some orbits on the center-stable and center-unstable manifolds are drawn by the green and blue arrows. The red arrows a, b represent stable walking (attracting periodic orbit). The curved and straight red arrows represent U (motion in single support phase) and T (jump at foot contact), respectively

Geometric characteristics
• Two boundary curves of D are almost parallel, and one of them is very close to W cs .
We also showed the following properties for B: • B is located inside D and is thinner than D; • B is V-shaped; • There are fractal-like slits in B and a stripe pattern in the cusp of the V-shaped region.  To investigate how these geometric characteristics of B are generated from D, we calculated the inverse images of D, S −n (D) (n = 1, 2, . . .). Figure 5 shows D, S −1 (D), and S −2 (D), which showed the following: , is V-shaped, and has a slit.

The formation mechanism of thin domain
The domain D is very thin, as shown in Fig. 4a. This fact is related to the so-called "λ-lemma", one of the most important theorems in the theory of dynamical systems. From this theorem, we can say the following (see the textbook by Robinson [15] for the exact statement of the theorem and its proof): • A region intersecting the unstable manifold of a saddle equilibrium point moves toward the stable manifold of the saddle under the backward flow (time-reversal of the flow); • When the region comes close to the stable manifold under the backward flow, the region becomes thinner due to the hyperbolicity of the dynamics near the saddle. Figure 6 illustrates how a region X moves and is deformed into thinner regions Y and Z under the backward flow. Because the equilibrium point of our model is a saddle-center and not a saddle equilibrium point, we cannot apply this theorem directly. However, a region intersecting the center unstable manifold near the equilibrium point moves close to the center stable manifold under the backward flow of the equations of motions in a way similar to Fig. 6 in the stable and unstable directions. As shown in Fig. 2b, D is obtained by the intersection of T (H ) and the backward orbit whose initial point is in H . Therefore, D becomes thin along the center stable manifold, as shown in Figs. 4A, B, and 6. This explains why the domain is very thin.

The formation mechanism of V-shaped region and slits
Since the sequence of the inverse images S −n (D) (n = 1, 2, . . .) converges to the basin of attraction, it is important to clarify how the geometric structure of the inverse images is formed, and hence the shape of the basin of attraction. In particular, we show why the inverse images are V-shaped and why they have slits and stripe patterns. First, S −1 (D) is V-shaped in the thin region D, as shown in Fig. 5. The formation of this shape is explained as follows (Fig. 7a, b). D is a thin region along the center stable manifold, as described in Sect. 3 Fig. 7a). Since D is a thin strip, T −1 (D) is also a thin strip. Since U −1 is given by the backward flow of the equations of motion, T −1 (D) is strongly expanded along the direction of the stable manifold and strongly contracted along the direction of the unstable manifold by U −1 , as shown in Fig. 6. In addition, the flow of θ 1 (t) becomes slow near the equilibrium ) is moved and deformed by the backward flow to S −2 (D) = U −1 (T −1 (S −1 (D))) point due to the saddle hyperbolicity and this causes nonuniform expansion by U −1 near W cu . Therefore, S −1 (D) becomes V-shaped, as shown in Fig. 7b.
We can also explain why the inverse image S −2 (D) has a slit in the same way (Fig. 7c, d). We can also give a similar explanation for the stripe pattern, which is formed by the repeated expansion of nested regions.
This explains how the slits and stripe patterns in the basin of attraction are formed. The relative positions of T −1 (D) and the center unstable manifold and the hyperbolicity near the saddle determine the geometric characteristics of the basin of attraction.

Basin of attraction in active dynamic walking model
In this section, we investigate the basin of attraction for the active dynamic walking model in comparison with the passive dynamic walking model. We will see the same mechanisms as in Sect. 3 for this model.   Figs. 10b and 11b. In addition, these figures also explain why Fig. 8a, b have different shapes. The difference of the relative position between W cu and T −1 (D) induces the difference between these two figures. Figures 10c, d and 11c, d explain S −2 (D) of these two slices. The repeated application of this mechanism forms the "horn-like" shape of the basin of attraction. These two types of structures (and also other types of structures) coexist in the three dimensional space T (H ).

Geometric characteristics
These arguments show that the two different walking models share the common formation mechanism of the basin of attraction studied in this paper.

Transition between "two horns without head" and "two horns with head"
When the slicing position changes from φ = −0.342 to φ = 0.0535, the slice of the basin of attraction changes from "two horns without head" to "two horns with head". Figure 12a, b show the slices of the basin of attraction at φ = −0.0780 and φ = −0.0720. Between these two slices, the horn-like region merges and the geometric structure changes drastically. As mentioned in Sect. 4.2, the difference between Fig. 8c and d can be explained by the geometric change of S −1 (D). Therefore, we focus on the geometric change of S −1 (D) when the slicing position changes. Figure 13 shows the change of S −1 (D) when the slicing position changes gradually. From these figures, we observe that two horns approach each other from Fig. 13a and b, two horns are merged between Fig. 13b, c, and S −1 (D) becomes larger from Fig. 13c, d. a b c d The change of the shape of S −1 (D) comes from the transition from Fig. 10a to Fig. 11a. Between Fig. 13b and c, the boundary of S −1 (D) becomes tangent to the center unstable manifold, and the shape changes drastically. S −2 (D), S −3 (D), . . . also change similarly as the slicing position changes, and finally the shape of the basin of attraction changes drastically as in Fig. 12a, b.
We created a movie 3 about the structural changes of D, S −n (D) for n = 1, . . . , 5 and the basin of attraction. You can see the details of the change of the shapes of these regions from the movie.

Conclusion and future works
In the present study, we clarified the formation mechanism of the basin of attraction for two bipedal walking models by focusing on the intrinsic hyperbolicity in the governing dynamics and based on the viewpoint of the theory of dynamical systems. We showed that the formation mechanism of the basin of attraction for the passive dynamic walking model in [13] is not specific to that model, but is applicable also to the active dynamic walking model studied in this paper.
The thin basin of attraction of the passive dynamic walking model is closely related to the one-dimensional instability of the upright equilibrium. The thin basin of attrac-tion of the active dynamic walking model is also closely related to the one-dimensional instability of the periodic orbit. In these models, there is a codimension one centerstable manifold. For the both models, the one-dimensional instability comes from the saddle hyperbolicity of an inverted pendulum. Although the present study focuses on these two models, our analysis strongly suggests that the results are not specific to them, but are widely applicable to more general bipedal walking models due to the intrinsic saddle hyperbolicity of bipedal walking.
The detailed structures of the basin of attraction are, however, rather different between the passive and active dynamic walking models. We showed that the difference came from the relative position of the center-unstable manifold and T −1 (S −n (D)) (n = 0, 1, . . .). The formation mechanism of basin of attraction by the relative position of the center-unstable manifold and inverse images can explain various types of shapes of basin of attraction. Furthermore, we can use the mechanism to deform the basin of attraction by regulating the relative position of those objects under additional controllers and support systems. Therefore, the present study will hopefully contribute to improving the understanding of the stability mechanism in human walking and to producing design principles for the control of biped robots and walking support systems. This will be a subject of our future work.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.