Fractal mechanism of basin of attraction in passive dynamic walking

Passive dynamic walking is a model that walks down a shallow slope without any control or input. This model has been widely used to investigate how humans walk with low energy consumption and provides design principles for energy-efficient biped robots. However, the basin of attraction is very small and thin and has a fractal-like complicated shape, which makes producing stable walking difficult. In our previous study, we used the simplest walking model and investigated the fractal-like basin of attraction based on dynamical systems theory by focusing on the hybrid dynamics of the model composed of the continuous dynamics with saddle hyperbolicity and the discontinuous dynamics caused by the impact upon foot contact. We clarified that the fractal-like basin of attraction is generated through iterative stretching and bending deformations of the domain of the Poincaré map by sequential inverse images. However, whether the fractal-like basin of attraction is actually fractal, i.e., whether infinitely many self-similar patterns are embedded in the basin of attraction, is dependent on the slope angle, and the mechanism remains unclear. In the present study, we improved our previous analysis in order to clarify this mechanism. In particular, we newly focused on the range of the Poincaré map and specified the regions that are stretched and bent by the sequential inverse images of the Poincaré map. Through the analysis of the specified regions, we clarified the conditions and mechanism required for the basin of attraction to be fractal.


Introduction
Passive dynamic walking is a model that walks down a shallow slope without any control or input [27], which is useful for investigating the mechanism of generating stable walking from a dynamic viewpoint. This has been widely used to examine how humans walk with low energy consumption [6,7,11,[24][25][26] and to provide design principles for energy-efficient biped robots [5, 9, 10, 21-23, 36, 37]. However, the basin of attraction is very small and thin and has a fractal-like complicated shape [2,29,33], which makes it difficult to produce stable walking. Furthermore, chaos appears in the walking behavior through a period-doubling cascade by increasing the slope angle [15], which makes producing stable walking even more difficult. Meanwhile, the basin of attraction shows a fractal-like shape, even without period doubling. In other words, the fractal-like basin of attraction appears even for a single attractor. Although this indicates that a different mechanism from the period doubling of the attractor induces a fractal-like basin of attraction, the mechanism is unclear.
In our previous study [29], we used the simplest walking model [13] for the analysis of passive dynamic walking and clarified the formation mechanism for the basin of attraction based on dynamical systems theory by focusing on the hybrid dynamics of the model composed of the continuous dynamics generated by the equations of motion during the swing phase with saddle hyperbolicity and the discontinuous dynamics generated by the impact upon foot contact. Specifically, we found that the fractallike basin of attraction is generated through iterative stretching and bending deformations by sequential inverse images of the Poincaré map for the collection of initial conditions from which the model can walk at least one step, which corresponds to the domain of the Poincaré map. However, whether the fractal-like basin of attraction is actually fractal, i.e., whether infinitely many self-similar patterns are embedded in the basin of attraction, is dependent on the model parameters, such as the slope angle [2,33]. The mechanism that determines whether the basin of attraction is fractal remains unclear.
In the present study, we improved our previous analysis in order to clarify the mechanism. In particular, we newly focused on the range of the Poincaré map, which corresponds to the collection of states after the model walked one step starting from the domain, and specified the regions that are stretched and bent by the sequential inverse image of the Poincaré map. Through analysis of the specified regions, we clarified the condition and mechanism required for the basin of attraction to be fractal.

Model
In the present study, we used a compass-type model (figure 1(A)) for the analysis of passive dynamic walking. This model has two legs (rigid links), the lengths of which are both l, connected by a frictionless hip joint. Here, θ is the angle of the stance leg with respect to the slope normal, and ϕ is the relative angle between the stance and swing legs. The mass is located only at the hip and the leg. The hip mass is M, and the leg mass is m. The leg mass is located at a distance b from the hip joint. In addition, g is the acceleration due to gravity. This model walks on a slope of angle γ without any control or input.

Structure of phase space by hybrid dynamics
In the present study, we focused on the simplest walking model, where m/M → 0 and b/l → 1 [13], because the dynamical characteristics remain almost unchanged [29]. This model is governed by hybrid dynamics composed of the continuous dynamics generated by the equations of motion during the swing phase and the discontinuous dynamics generated by the impact upon foot contact.
The equations of motion are given bÿ The impact upon foot contact yields the following relationship: where * − and * + are the state * just before and after the foot contact, respectively. The important property of this relationship is that the state just In addition, from (4), the state just after foot contact also satisfies 0 < θ + < π/2 ( 9 ) However, note that since the state just after foot contact is independent ofφ − , (5) generates no condition. This hybrid dynamic system determines the structure of the phase space, as shown in figure 2(A). Here, H is the section defined by the foot contact conditions (3)-(5) and forms a three-dimensional space in the four-dimensional phase space, and T is the jump in the phase space from the state just before foot contact to the state just after foot contact, as defined by (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). Moreover, U is the map from the start of a step to the next instance of foot contact. In other words, U is the map from T(H) to H, as defined by the equations of motion (1) and (2). The Poincaré map S is defined on the Poincaré section T(H) by This Poincaré map S represents one step, and an attractor of S represents stable walking. The basin of attraction of S is the main topic of the present paper. Here, S is parameterized by one parameter γ. In particular, S has an attracting fixed point at 0 < γ < 0.015, and there is a period-doubling cascade to chaos for 0.015 < γ < 0.019 [13].

Domain of Poincaré map and basin of attraction
We define D n (n = 1, 2, . . .) as the collection of initial conditions on T(H) from which the model walks at least n steps. This satisfies D n+1 ⊆ D n (figure 2(B)), which means that when the initial condition is in D n but out of D n+1 , the model will fall down at the n + 1th step. Since the Poincaré map S represents walking one step, S(D n ) indicates the state on T(H) after the model walked one step starting from D n . Since the model can walk at least n − 1 steps from S(D n ), the following condition is satisfied: Since the domain D of S on T(H) represents the collection of initial conditions on T(H) from which the model walks at least one step, D is identical to D 1 .
Using the inverse image of S, we can write D n = S −1 (D n−1 ). However, S −1 acts only on a part of D n−1 , as shown in figure 2(B), as clarified in the following section. First, the range R of S on T(H) is given by R = S(D 1 ) because D 1 is the domain of S, which corresponds to the collection of states after the model successfully walked one step starting from all states on T(H). This means that the state after each step must be in R unless the model falls down. The following equation is satisfied: We prove this below. First, since D n ⊆ D 1 , S(D n ) ⊆ R. Based on this consideration and (11), This contradicts d S(D n ). Therefore, this assumption is not satisfied. Since any state in D n−1 ∩ R is in S(D n ), we obtain (12). Therefore, instead of D n = S −1 (D n−1 ), we use In the same manner, instead of D n = S −n+1 (D 1 ), we use Since the model walks at least n steps from D n , D n approximates the basin of attraction as n increases. We confirmed the convergence by comparing D 100 and D 200 using 10 4 × 10 4 initial states. In order to clarify these geometric characteristics, we used θ +θ and θ −θ for the axes for γ = 0.001 and 0.013 in figures 3(B) and (D), respectively, as in our previous study [29]. Here, D 2 and D 3 are V-shaped, and D 3 has a thin slit (note that a V-shaped region indicates that the basin has one large slit). The basins of attraction have multiple slits and are complicated. In particular, while the basin of attraction has only a few slits for γ = 0.001, the basin has a number of slits, and self-similar patterns are embedded for γ = 0.013. In order to quantitatively clarify these properties, we examined how many gaps D 1 , D 2 , D 3 , and the basin of attraction have on the line θ −θ = 0.5. As a result, there are zero, one, two, and four gaps in D 1 , D 2 , D 3 , and the basin of attraction, respectively, for γ = 0.001. In contrast, there are zero, one, two, and infinitely many gaps in D 1 , D 2 , D 3 , and the basin of attraction, respectively, for γ = 0.013. We further investigated the number of slits in D 1 , D 2 , D 3 , and the basin of attraction for γ based on the gaps on θ −θ = 0.5 in figure 3(E). The number of slits in D n increases as n increases, and that in the basin of attraction increases exponentially as γ increases. Fractal structures appear in the basin of attraction over γ ≈ 0.0075.

D
approximates the basin of attraction as n → ∞, we begin with D 1 ∩ R to investigate the formation mechanism of the basin of attraction. As shown in figure 3(E), the numbers of slits in D 1 , D 2 , and D 3 remain unchanged for γ. Therefore, we assume that the mechanism is common for γ when n is small. We used γ = 0.013 to show the results below.
Since R = S(D 1 ) = T(U(D 1 )), we first examine U(D 1 ). In particular, the boundaries of U(D 1 ) are θ = 0, θ = −π/2, and 2θ =φ from the foot contact conditions (4) and (5), as shown in figure 4(A). In addition, since D 1 does not intersect with W cu , D 1 and U(D 1 ) are on the same side with respect to W cu , and U(D 1 ) also has a boundary near W cu (strictly speaking, D 1 intersects with W cu in a small range of 0 < θ < γ, but has no influence in the formation of D 1 ∩ R and so is ignored). Figure 4(B) shows the result for U(D 1 ) projected onto the θ-θ plane (figure 4(C) uses θ +θ and θ −θ for the axes in order to clarify the geometric characteristics).
Since T is a one-to-one mapping for [θθ], the boundaries of U(D 1 ) on the θ-θ plane become the boundaries of R by T, as shown in figure 4(A). Figure 4(D) shows the result of R obtained from (6) (figure 4(E) uses θ +θ and θ −θ for the axes to clarify the geometric characteristics). The boundaries a 1 b 1 and c 1 d 1 of D 1 ∩ R are obtained by applying T to the boundaries near W cu and 2θ =φ, respectively, of U(D 1 ). figure 4(D). is thin, as shown in figure 5(A), we extract a line segmentP 1Q1 from T −1 (D 1 ∩ R). However, sincė ϕ in T −1 (D 1 ∩ R) is not uniquely determined, we consider T −1 (D 1 ∩ R) as a quadrangular prism, the height of which is in theφ direction, as shown in figure 5(C). Then, line segmentP 1Q1 is considered to be a plane, which we call plane Z. We apply U −1 to the plane Z. Since U −1 is the map from H to T(H) and T(H) is a two-dimensional surface that has two constraint conditions (7) and (8) in the four-dimensional phase space, U −1 is applicable only to points in the plane Z that simultaneously satisfy the two conditions when the points are moved in the phase space in the time reverse direction using the equations of motion (1) and (2), as shown in figure 5(C). These points determineφ in T −1 (D 1 ∩ R). Figure 5(D) shows the result for the collection of the points in the plane Z indicated by the curveê 1f1ĝ1 to which U −1 is applicable. We obtained the curve e 2 f 2 g 2 by applying U −1 to this curveê 1f1ĝ1 , as shown in figure 5(E). (Figure 5(F) uses θ +θ and θ −θ for the axes to clarify the geometric characteristics.) Note that * 2 (except for D 2 ), such as e 2 , is used for D 2 . In order to obtain the curves e 2 f 2 g 2 andê 1f1ĝ1 in figures 5(D) through (F), we linearized the equations of motion (1) and (2) for θ and ϕ because the walking behavior appears around the saddle [γ 0 0 0]. Therefore, there are differences from the exact solution, as the approximately obtained curve e 2 f 2 g 2 is not inside D 2 (figures 5(E) and (F), see the appendix A for details).

Characteristics of
The curve e 2 f 2 g 2 , specifically the curve e 2 f 2 is bent to be V-shaped, as shown in figure 5(E). (Figure 5(F) uses θ +θ and θ −θ for the axes to clarify the geometric characteristics.) In order to examine where in D 1 ∩ R the curve e 2 f 2 g 2 is moved from by S −1 (= U −1 • T −1 ), we investigate where in D 1 ∩ R the curveê 1f1ĝ1 is moved from by T −1 . Since the curveê 1f1ĝ1 is in T −1 (D 1 ∩ R), the curve moves in D 1 ∩ R by T. Figures 4(D) and (E) show the result indicated by the curve P 1 Q 1 . This shows that when S −1 is applied to the curve P 1 Q 1 in D 1 ∩ R, two curves f 2 e 2 and f 2 g 2 are obtained in T(H). Since D 1 ∩ R is thin, as shown in figure 4(D),  the curve P 1 Q 1 approximates D 1 ∩ R. Therefore, the process to obtain the V-shaped curve e 2 f 2 g 2 from the curve P 1 Q 1 explains the process by which D 1 ∩ R is transferred to D 2 .
Figures 6(A) through (D) show a schematic diagram of the summary by which to obtain D 2 from D 1 . Specifically, D 1 ∩ R is extracted from D 1 ( figure 6(B)), and two regions are generated by S −1 , one of which is stretched and bent (figure 6(C) left), and the other of which is only stretched (figure 6(C) right). These regions are connected at the boundaries a 1 b 1 to form D 2 (figure 6(D)).
Next, we move to D 3 = S −1 (D 2 ∩ R). Since D 2 ⊆ D 1 , the deformation from D 2 to D 3 (figures 6(D) through (G)) is the same as that from D 1 to D 2 (figures 6(A) through (D)). Since D 2 is V-shaped (figure 6(D)), the extracted D 2 ∩ R is also V-shaped and has a large slit at the boundary c 2 d 2 (figure 6(E)). The large slit becomes slits at the boundary c 2 d 2 in D 3 by the deformation (figures 6(F) and (G)). Although figure 6(F) (right) has a slit, it is far from R and so is ignored (figure 6(G)).
When the slit generated in D n reaches, but does not penetrate, R, one slit is added in D n+1 , as observed in the process from D 2 to D 3 . In contrast, when the generated slit in D n does not reach R, the number of slits in D n+1 remains unchanged, as shown in figure 7. Moreover, since D n+1 ⊆ D n , it is possible that the slit becomes deeper as n increases to reach R and create a new slit. Therefore, when the generated slit in D n does not penetrate R, the number of slits of D n+1 increases by one or remains unchanged.

Appearance of a fractal
We consider the cases in which the generated slit in D n penetrates R for the first time at n = N, as shown in figure 8(A). (There may be multiple slits that do not penetrate D N ∩ R to the left and right of the generated slit, but because they do not affect the explanation below, they are not shown in figure 8.) By applying S −1 to D N in the same manner as in figure 6, a penetrating slit appears close to the outer edge of the V-shaped D N+1 , as shown in figure 8(D). In addition, since D n+1 ⊆ D n once a slit penetrates R, the slit penetrates R for n > N. Furthermore, the penetrating slit close to the outer edge of D N+1 generates a slit that penetrates R near the right edge of D N+1 ∩ R, as shown in figure 8(E). As a result, a penetrating slit also appears close to the inner edge of the Vshaped D N+2 , as shown in figure 8(G). Furthermore, Figure 8. Schematic diagram of the process by which to deform D n and generate penetrating slits after the generated slit penetrates R for the first time at n = N in (A). A penetrating slit is generated close to the outer edge in D N+1 in (D). A penetrating slit is generated close to the inner edge in D N+2 in (G). A penetrating slit is generated close to the generated slit in D N+3 in (J). the penetrating slit produces another penetrating slit in D N+3 near the slit generated by the large slit of D N+2 due to the V-shape, as shown in figure 8(J). This slit also penetrates R. These penetrating slits produce new penetrating slits near the edge, and the number of slits increases at an accelerated rate as n increases. As a result, a fractal basin of attraction appears. Figures 9(A) through (E) show D 4 to D 8 for γ = 0.013. At N = 5, the generated slit penetrated R for the first time ( figure 9(B)). After that, the penetrating slits were generated close to the outer edge at n = 6 (figure 9(C)), near the inner edge at n = 7 (figure 9(D)), and close to the generated slit at n = 8 (figure 9(E)) in that order. Therefore, infinitely many slits are generated and the fractal basin of attraction appears, as shown in figures 3(C) and (D).

No fractal appears
Next, we consider the cases in which no slits in D n penetrate R, even when S −1 is applied several times. In particular, suppose that n is so large that D n converges and also suppose that D n has a slit that does not reach R. In this case, since the number of slits does not change even when S −1 is applied, the basin of attraction does not have a fractal structure and has a finite number of slits. Figures 10(A) and (B) show D 50 and D 51 , respectively, for γ = 0.001. Here, D 50 has four slits, and the leftmost slit does not reach R. As a result, D 51 has four slits as in D 50 . In addition, D 50 and D 51 have no difference and are identical to the basin of attraction ( figure 3(B)), which we confirmed by comparing the regions using 10 4 × 10 4 initial states, and so converge. Therefore, the basin of attraction does not have a fractal structure.

Stability and basin of attraction
Bipedal walking has intrinsic instability due to saddle dynamics, and clarifying the mechanism by which walking can be stabilized is important. Passive dynamic walking is a useful model to examine the mechanism from a dynamic viewpoint. In order to clarify the stabilization mechanism, investigating both the stability and basin of attraction is crucial. However, while previous studies have focused on the stability by the eigenvalue analysis of the linearized Poincaré map around the fixed point on the Poincaré section [8,13,15,16,21,[35][36][37], the basin of attraction has not been well studied. This is partly because while eigenvalue analysis allows us to easily investigate the stability, no general analytical method has been provided for investigating the basin of attraction. We used an analytical approach based on dynamical systems theory to clarify a specific property embedded in the basin of attraction, which is useful to further investigate the characteristics of the basin of attraction in walking. While passive dynamic walking has no control or input, the use of control and input changes the dynamic characteristics of walking and also varies the stability and the basin of attraction [3,4,30,31,34]. We would like to improve and clarify our analysis in the future.

Initial-value sensitivity and convergence to attractor
The Poincaré map S represents walking one step, and slits are generated by applying the inverse image S −1 many times to the region from which the model walks at least one step. These slits come from the large slit of the V-shaped D 2 ( figure 6). When there are only a finite number of slits in the basin of attraction, the generated slit in D n does not reach R and is not used in D n+1 (figure 7). Therefore, these slits are not stretched much. In contrast, once the generated slit in D n penetrates R, the generated slits for n N are stretched greatly and create stripe patterns by producing penetrating slits, especially close to the basin boundary (figure 8). These penetrating slits become thinner as n increases. Since slits indicate a region in which the model will fall down, whether the model continues to walk or not becomes very sensitive around the basin boundary. Furthermore, since penetrating slits become thinner, two states located at different sides of the large slit of D 2 become closer as S −1 is applied to the two states many times. This means that two states located at different sides of a thin slit in the basin of attraction move away from each other as S is applied many times and the two states come to reach different sides of the large slit of D 2 .
When the basin of attraction is fractal, there are infinitely many penetrating slits close to the basin boundary. Therefore, when the model walks from an initial state near the boundary on the basin of attraction, there are numerous penetrating slits between the initial state and the attractor, and the model must traverse the slits for the state to approach the attractor. The model must walk at least the steps that are required to generate the penetrating slits by applying S −1 . Therefore, the model takes a long time to approach the attractor, depending on the initial state.

Limitations of our analysis
In the present study, we clarified that the fractal basin of attraction appears when the generated slit in D n penetrates R and that fractal basin of attraction does not appear when the generated slit in D n does not reach R for an n so large that D n converges. However, it is possible that the generated slit in D n reaches, but does not penetrate, R for so large n that D n converges. In this case, although it is not at an accelerated rate, the number of slits increases as n increases. While infinitely many slits appear in the basin of attraction, no penetrating slits are generated. Although our analysis does not exclude this possibility, our simulation results did not show such a case for any γ.
We used the simplest walking model for the analysis of passive dynamic walking, i.e., we assumed the extreme case m/M → 0 and b/l → 1 for the compass-type model [13]. Therefore, we did not explain the mechanism of the basin of attraction for general models of passive dynamic walking. However, the period-doubling cascade to chaos appears and the fractal basin of attraction is observed without the period doubling even when the extreme case is not assumed [2]. This suggests that similar mechanisms to those observed herein are embedded in general models of passive dynamic walking.

Biological relevance
The fractal appears in human walking, especially in the gait rhythm [17][18][19][20]. However, unlike passive dynamic walking, human walking is generated through the control. The basin of attraction of compass-type models used in passive dynamic walking is enlarged by the control and the number of slit changes [30,31,34]. However, the stance leg during human walking is almost straight and rotates around the foot contact point like an inverted pendulum [28]. In addition, the stance and swing legs are switched by the foot contact and lift off. Therefore, saddle instability and hybrid properties are inevitable in the gait dynamics, as in passive dynamic walking, and the stretching and bending deformation remains crucial for the formation mechanism of basin of attraction. In fact, our previous study [30,31] showed that even when a controller inspired by spinal central pattern generators [32] is incorporated in a compass-type model, the basin of attraction has slits due to the deformation.
Human walking is generated through the central nervous system and the body mechanical system. Fractal properties are reduced by aging and pathological disorders such as Parkinson's and Huntington's diseases [17,20]. A simple neuromechanical model demonstrated that fractal properties are reduced by changing the motor control model to emulate the pathological disorder [12]. These properties suggest that the neural system contributes to the fractal in human walking. In contrast, the body mechanical system also has potential to contribute to the fractal in human walking. Passive dynamic walking exhibits a chaos attractor depending on the model parameter [13,15] and shows a fractal basin of attraction even for the single attractor as shown in the present study. The steady state of a dynamical system with a single attractor never shows a fractal, but instead shows regular behavior, unless the system is disturbed. However, when the dynamical system is specific and has a fractal basin of attraction, fractal behavior can be induced by a disturbance or noise without fractal properties. In fact, the fractal appears in walking of compass-type models with a controller and noise without fractal properties [1,14]. Even for passive dynamic walking, the mechanisms for fractal and non-fractal basins of attraction clarified in the present study will provide useful insights for understanding human walking. The analysis of measured human data has limitations for elucidating the underlying mechanism in human walking, and physical models are useful to overcome the limitations.
where −Δ is used as the negative time to analyze U −1 (figure 2(A)). Θ(0), Θ(−Δ), and −Δ correspond to the state just before foot contact, the state just after foot contact, and the duration of a step, respectively. Θ(−Δ) gives the deformation of T −1 (D 1 ∩ R) by U −1 .
Since Θ(−Δ) is in T(H), the following equations are satisfied from (7) In addition, since Θ(−Δ) is in H, the following equation is satisfied from (3): In order to approximately solve (16), we linearize the equations of motion (1) and (2) The solution is obtained by where C 1 , C 2 , K, and φ are the integration constants (0 φ < 2π). Here,ψ on plane Z is obtained bẏ C 1 and C 2 are determined by the initial conditions of θ andθ, as follows: In contrast, K and φ are determined by the initial conditions of θ, ϕ,θ, andφ.