Passive walking biped robot model with flexible viscoelastic legs

To simulate the complex human walking motion accurately, a suitable biped model has to be proposed that can significantly translate the compliance of biological structures. In this way, the simplest passive walking model is often used as a standard benchmark for making the bipedal locomotion so natural and energy efficient. This paper is devoted to a presentation of the application of internal damping mechanism to the mathematical description of the simplest passive walking model with flexible legs. This feature can be taken into account by using the viscoelastic legs, which are constituted by the Kelvin–Voigt rheological model. Then, the update of the impulsive hybrid nonlinear dynamics of the simplest passive walker is obtained based on the Euler–Bernoulli’s beam theory and using a combination of Lagrange mechanics and the assumed mode method, along with the precise boundary conditions. The main goal of this study is to develop a numerical procedure based on the new definition of the step function for enforcing the biped start walking from stable condition and walking continuously. In our previous work, it was shown that the period-one gait cycle can be produced by adding the proportional damping model with high damping ratio to the elastic legs. The present paper overcomes the practical limitations of this damping model and similarly demonstrates the stable period-one gait cycles for the viscoelastic legs. The study of the influence of various system parameters is carried out through bifurcation diagrams, highlighting the region of stable period-one gait cycles. Numerical simulations clearly prove that the overall effect of viscoelastic leg on the passive walking is efficient enough from the viewpoint of stability and energy dissipation.


Introduction
Studying passive bipedal walking is a fascinating field in the research of legged walking robots for its achievement in human locomotion mimicry. Passive dynamic walking is one of the locomotion strategies which can settle into a stable steady and symmetric gait on a gentle slope without requiring any active control using gravity and dynamic nature of the swing leg [1]. It is clear that the passive walking gait is highly energy efficient, and thus, investigating these natural gaits may lead us to gain a deeper understanding of human walking.
From the mathematical point of view, the complex impulsive hybrid nonlinear dynamics of the passive walking is influenced greatly due to changes of the biped robot model parameters. Therefore, optimal designs of passive walkers are one of the controversial topics that have attracted a lot of interest. Two wellknown passive bipeds with simple mechanical structure, the compass-like biped robot and the point-foot walker were then conceived by Goswami et al. [2] and Garcia et al. [3], respectively. They have numerically shown the existence of stable period-one gait cycles on a range of steeper slopes and then the cascade of period-doubling bifurcations as a route to chaos by varying the slope angle and even the structural parameters of the models [4][5][6]. After that, many follow-on robotic scholars have performed studies on the most popular passive walkers with knees [7], with a torso on hip [8], with different foot shapes [9,10], with stabilizing arm [11], by changing walking ramp surface into stairs [12] and also by considering the friction forces [13]. Early considerable research works conducted under the passive dynamic walking are mainly an analysis of its dynamic characteristics, namely limit cycle, Poincaré map, stability, bifurcations, chaos, control and so on, as in [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30], just to mention a few. It is important to note that the literature review of the passive dynamic walking theory is too long and here, in accordance with the aim of the paper, these references are mentioned such as the examples of recent studies in this field. Among the aforementioned works, Gritli and his co-workers have played a key role in the continuous study and the development of passive bipedal robots in recent years [14][15][16][17][18][19][20][21][22][23][24][25]. They revealed the existence of the period-three route to chaos in a compass-gait biped robot through a cyclicfold bifurcation [14,15]. Furthermore, they used the developed expression of the hybrid Poincaré map to analyze chaos and bifurcations and also to compute the spectrum of Lyapunov exponents in the compass-gait model under the OGY-based state-feedback controller [16][17][18][19][20][21]. Several forms of gait bifurcations, including the sub-critical and super-critical period-doubling bifurcations, the saddle-saddle bifurcation, the saddle-flip bifurcation and the Neimark-Sacker bifurcation were found in their works. As reported in [25], two control techniques based on tracking a desired trajectory via the feed-forward-plus-PD controller have been proposed to provide a remarkable energetic efficiency compared to the OGY-based control method. Moreover, in an effectual way, some research works focus on the analysis of walking behaviors of bipeds under a leg length asymmetry [26,30]. The interested readers are also referred to the papers in [31,32] for an extensive treatment of design and dynamic analyses of various passive and semi-passive biped walkers. It is worth noting that many excellent theoretical and experimental results are achieved in these literatures.
In the strive for simulating passive elastic mechanisms utilized in human walking, a group of researchers has decided to use compliant elements at different positions of the rigid biped such as torsional spring in joints [33], ankle [34] and torsional stiffness in the hip joint [35]. An extended model of the compass-gait biped with series-elastic ankle actuation was investigated in [36] which exhibited asymptotically stable walking on flat ground. The spring added to the ankle joint of the stance leg is stimulated toe push-off like human. Deng et al. [37] studied the walker's specifications by adding a fully passive torso to a passive biped via torsional springs. The presented results showed two stable gait cycles with different properties and also improvements in gait performances. A chaotic particle swarm optimization algorithm was applied by Wu et al. [38] to the optimal design of passive walker with hip torsional stiffness and damping at multi-variable level. A modified model of Garcia's basic passive biped with a sole foot and a combination of elastic and viscous elements at both hip and ankle joints was presented in [39]. By studying nonlinear dynamics and chaotic behavior of this biped model, the authors observed two possible patterns of motion similar to multi-pattern humanoid gait and also investigated the influence of system parameters on the stability of motion. These examples have been proven that adding elastic elements to the biped based on biomechanical compliance has advantages over the rigid link walker and acceptable results of flexibility in passive bipedal walking; however, it requires more sophisticated nonlinear mathematics to describe the impulsive hybrid passive dynamic walking and increases the designing leg mechanism complexity as anticipated.
One of the main motivations behind the study of the passive walkers is to obtain additional insights into the biomechanical design principles of human bipedal walking and their implementation in robotic platforms. In addition to the benefits of the local compliant elements in hard-bodied bipedal robots, there is still a gap between the structure of walkers and the human anatomy from a technical robotic design perspective. In contrast to the conventional rigid legs, the intrinsically flexible nature and extensible materials of the physical structure can mimic the properties of the biological tissues and muscles of the human leg and then the passive walker appears more realistic. Therefore, using lightweight flexible legs can help the walker generate natural and very efficient walking pattern, adapt to uneven terrains and increase adaptability. Nevertheless, there are few researches dedicated to the flexible passive walkers having elastic legs instead of rigid links. For instance, Shen et al. [40] have studied walking dynamic responses, complicated contact modes and impact-induced waves of a new compliant biped robot on an elastic walking surface. Great attention of this work has to be placed on the effects of sliding contact such as dynamic self-locking phenomenon on walking stability of the compliant walker. The other research work in [41], our earlier study, developed a modified flexible model of the Garcia's simplest passive walker with elastic legs. The unwanted vibrations induced by these undamped legs exhibited non-periodic gait cycle. After adding damping, we found proper initial conditions for period-one gait cycles and then investigated the influences of the walker's parameters through bifurcation diagrams.
It is observed that in the above cited references, the energy dissipation of the pure elastic legs was either ignored or considered as a simple viscous damping, while, in practice, due to the internal damping forces, the vibrations are damped. It turns out that it is highly desirable to use passive damping technology for flexible legs with simple computational tasks. The use of viscoelastic materials in structures such as flexible manipulators [42] has led to efficient designs by minimizing undesired vibrations. Hence, the simple, inexpensive and widely available viscoelastic materials are the best biomechanical models like artificial muscles to get a good design of highperformance passive walkers. In robotic platforms and prosthetic systems, when the goal is human-like walking behavior, the leg mechanism should simulate the compliant properties of muscle actions. In this respect, a relevant further step in the direction of improving human bipedal locomotion would be to investigate the effective physiological compliance of the robotic leg. To the best of authors' knowledge, no previous study has been made to utilize inherently dissipative materials in known passive walkers. Consequently, the present paper aims to implement the intrinsic viscoelastic properties of human muscles in flexible passive biped robots in order to enhance the natural look of walking and deal with the vibrations of the former elastic legs. Following our recent study in [43], the main contributions of this paper lie in modeling the flexible simplest passive walker, numerically finding proper initial motion state for its periodone gait cycles and discussing the effect of specific parameters on walking manner, but by proposing the viscoelastic legs.
The remaining part of this article is organized as follows. Section 2 presents the mechanics model, the corresponding formulation as well as the process of obtaining suitable initial conditions for stable walking. In Sect. 3, the results of numerical simulations and comparisons are given and discussed. Finally, Sect. 4 is devoted to the concluding remarks.

Walking pattern analysis
The simulation of one step comprises swing phase, a heel-strike collision and a set of discrete update laws. Relying on the author's study in [41], the goal of this section is to give a brief overview of the simulation procedure for one walking step of the flexible simplest passive walker with viscoelastic legs. Besides many similarities with elastic legs, using the viscoelastic legs will lead to have new assumptions. The following procedure is fully documented in [41] and the interested readers are referred to this publication for a comprehensive study. Figure 1 illustrates the Garcia's simplest walking model [3], a standard benchmark for studying periodic gaits and the modified flexible model. The new introduced model consists of two identical compliant legs, connected by a frictionless hip joint of point mass M that allows for rotational motion in the 2-D plane. The significant physical parameters and using variables for the numerical simulation are shown in Fig. 1.

Physical description
The legs are considered to have homogeneous material properties with same length L, mass density q, moment of inertia I and constant circular crosssectional area A. Using the classical Euler-Bernoulli beam theory and including the structural viscoelastic effect, the slender legs of the walker are modeled based on the bending mechanism of flexible links. The key kinematic assumptions made for the convenient analysis of deformations are as follows.
• For each leg, the centerline of cross section is inextensible. • During deformation, the cross section area is assumed to remain planar and normal to the deformed axis of the leg. • Shear deformation, rotary inertia, the nonlinear geometric effect due to bending and warping are negligible. • Each leg undergoes small deflection.
• The flexible legs with Young's modulus E and viscosity coefficient g obey the Kelvin-Voigt rheological model.
The total deformation of each leg is composed of rigid rotation due to the angles h and / and transverse displacement of neutral axis. The position vectors in the local coordinate systems denoting an arbitrary point P of the cross sections of two legs are represented by where w st x 1 ; t ð Þ and w sw x 2 ; t ð Þ are the transverse displacements of the neutral axis (at point x and time t) compatible with the bending of the flexible legs. In this manuscript, the subscripts st and sw indicate the stance leg and the swing leg, respectively.
The whole movement process is a repetition of consecutive gaits constrained in two phases: a swing phase and an impact phase. Such kind of walking is started by launching the robot with a wise choice of initial conditions and the angle of the ramp. During a step, the stance leg is modeled as a hinge, connected to the floor. The swing leg is moving freely as the other end of a double pendulum. At about midstance, the swing leg will experience an impact with the ground (scuffing), which is inevitable for a walker with equal leg lengths. However, the scuffing problem of the swing leg is neglected in the theoretical study and there is no slip and no bounce at the contact point. Practically, leg shortening measure mechanisms such as a prismatic-jointed knee [43] would solve the scuffing problem without affecting the walker dynamics. After this short through-pass, the second time that the swing leg reaches floor level is regarded as heelstrike, the end of the step. The former swing leg makes a fully inelastic collision and becomes the new stance leg. Instantaneously, the former stance leg loses ground contact, and a new step begins. Similar to the traditional rigid simplest passive walker, the flexible leg of this modified model is in point contact with the surface during the walking process.

Derivation of the equations of the swing motion
The continuous dynamics of the passive walker as it travels down a shallow slope describes the bipedal motion between collisions and hence it is represented by a set of nonlinear differential equations. These swing equations of the viscoelastic Euler-Bernoulli legs, a set of coupled partial differential equations (PDE) and ordinary differential equations (ODE) along with the corresponding boundary conditions are derived using the extended Hamilton principle. First, the kinetic and potential energy and also the dissipated energy due to the damping effects should be obtained.
The position vectors of point P located at distances x 1 and x 2 from the local frame origins of the legs with respect to the contact point of the stance leg in the inertial coordinate system are given by The subscript l denotes the value of the w st x; t ð Þ at ¼ l, i.e., w st l; t ð Þ. The total kinetic energy of the system can thus be expressed as The potential energy of each flexible leg arises from two sources, gravity and internal bending moment. From the mechanics of materials, bending moment in an arbitrary cross section is approximated by Then, the simplified form of the total potential energy of the system can be formulated as The constitutive law of the classical Kelvin-Voigt model for the viscoelastic material is where _ e is the extensional strain rate. Using Eq. (6) and defining the bending moment in terms of stress, the non-conservative force of the viscoelastic beam is Therefore, the total virtual work done by this nonconservative force for the legs may be written as Now, the generalized Hamilton principle takes the form where L ¼ T À U. Substitution of the variation principle of Eqs. (3), (5) and Eq. (8) into Eq. (9) leads to yield the nonlinear differential equations of motion as And the leg's boundary conditions are described as at x 1 ¼ 0 As expected, according to the equations in [41], viscoelastic legs affect only deformation equations as seen in Eqs. (12) and (13). The resolution of the resulting continuous dynamics containing four coupled PDE and ODE is a challenging mathematical problem to solve analytically and computationally. To overcome, the assumed modes method (AMM), one of the popular formulation techniques for converting the PDE into ODE, is adopted here to perform the discretization of the developed equations.
We are interested now in using a combined Euler-Lagrange formulation and AMM to obtain the temporal equations of the swing motion. Accordingly, the approximate deflection subject to transverse vibrations of an elastic beam is expanded in a finite series of the form as [44] where u i x 1 ð Þ and w i x 2 ð Þ are the admissible functions that must satisfy the simply supported-free boundary conditions [45] relying on the geometric boundary conditions in Eq. (11) and n is the number of modes assumed in the approximation. Indeed, the stance leg and the swing leg are fixed to the contact point on the walking surface and the hip joint, respectively, while the other end of both legs is free. v i t ð Þ and p i t ð Þ are the time-dependent modal generalized coordinates and n is the number of modes.
By inserting the discrete formulation of determined kinetic energy and potential energy in Eqs. (3) and (5) as well as the virtual work of non-conservative damping force in Eq. (8) into the Lagrange's equation [2], which may be written as below, the reduced equations of motion written in terms of the generalized coordinates q are T is the vector of generalized coordinates whose variables are functions of dimensionless time M is the inertia matrix, F represents a global force vector, a vector of gravity forces and also includes the centrifugal and coriolis terms and F v is the vector of viscous structural damping term. Equation (18) results in two sets of nonlinear time varying coupled, second-order ordinary differential equations in the matrix form. The first set of these equations is associated with the rigid rotation of the flexible legs defined by h and /, and it is quite identical for both elastic and viscoelastic; the other set is associated with the elastic deformation defined by v i and p i and it slightly differs from the given equations of the pure elastic legs in [41].
To transform the governing equations into dimensionless form and generalize the research results, the following dimensionless parameters are introduced as Under the AMM and retaining a finite number, n ¼ 1 mode, upper triangular dimensionless components of the inertia matrix and the nonzero vectors F and F v in Eq. (20) are derived in ''Appendix A.''

Transition rules at impact conditions
In order to describe what happens to the walker during heel-strike, transition rules must be formulated to prevent trespassing of the surface slope during motion. Heel-strike collision is modeled as a rigid plastic impact with immediate and full contact to the floor. When the swing leg hits the surface, an instantaneous transition occurs where the previous stance leg loses contact and the stance and swing leg switch. After this transition, the biped will be in the exact same configuration as shown in Fig. 1a with the exception that the generalized coordinates h and / have switched places. Moreover, an update of the angular velocities should occur to prevent the walker from falling through the floor. These changes in leg position states and their velocities are algebraic equations which called transition rules as given below where the -and ? signs denote the time instant right before and right after impact, respectively. Following our recent study in [41], the above procedure of heel-strike collision happened in the rigid passive walker will be repeated for the modified models with linear elastic and also viscoelastic legs. But applying the flexibility effects of the legs even small deflection causes changes in transition rules, because these passive walkers are very sensitive to small variations in their modeling. Due to the impulsive nature of the impact forces, the robot configuration remains unchanged during the collision with the ground. Thus, the transition rule of positions of the viscoelastic legs is written as follows where I n and 0 n are the n Â n identity matrix and zero matrix, respectively; 0 is the 1 Â n zero vector and 0 T represents the transpose of it.
Since the impact forces are the only external forces affecting the biped, the angular momentum of the whole mechanism and the former stance leg (new swing leg) around the impact point and the hip is conserved before and after the collision. Presented below are two expressions for the small mass elements of the viscoelastic legs that exploit the angular momentum balance to derive the update of the angular velocities as given in [41].
The above position vectors in Eq. (22) are the distances of mass elements to the mentioned centers of rotation. Given the number of unknown velocities, there is still a need for exploring two other momentum expressions that specify the role of deflection (generalized coordinates v i and p i ) in transition rule of velocities. As stated before in [41], these novel expressions which are the time integrals of Eq. (10-3) and Eq. (10-4) during the time instances just before and just after collision used to update the time derivatives of w st x 1 ; t ð Þ and w sw x 2 ; t ð Þ for each leg. According to switching of the viscoelastic legs at the moment of impact, the conservation law of momentum of these expressions is defined as Then, by employing the AMM and using the orthogonality of the shape functions, Eqs. (24), (25) can be written in the discrete form as to calculate the transition rule of the velocities by following relation And finally, by merging the transition rules in Eqs. (23) and (27), we obtain the jump equations like so where 0 2nþ2 ð Þ is the 2n þ 2 ð ÞÂ 2n þ 2 ð Þzero matrix and also J s is a global transition matrix. In order for these jump equations to correctly be applied when the swing foot strikes the slope, the configuration of the robot that results in an impact must be determined. This configuration corresponds to the geometric collision condition known as the impact surface or switching surface. In contrast to the rigid model, a configuration of the flexible walker that leads to impact with the slope must satisfy the following relation The switching surface C is detected when the swing leg reaches the walking surface, while it is in front of the stance leg and is also moving downward. It is transparent that the transition rules at the impact collision are still valid by replacing the elastic legs with viscoelastic beams.

Hybrid dynamics
Using the continuous dynamics in Eq. (20), together with the impact surface defined by Eq. (29) and the algebraic transition rules in Eq. (28), a complete mathematical model of the walking biped can be illustrated as a diagram in Fig. 2.

Methodology for finding period-one gait cycles
In this subsection, we will focus on exploring proper initial conditions of period-one gait cycles of the simplest passive walkers with viscoelastic legs by an effective numerical method as proposed in our previous works [41,46].
Almost all the existing bipedal walkers based on passive dynamic walking start walking only through certain initial conditions; otherwise, the robot will fall down. Finding stable walking motions relies on these appropriate initial conditions which are located in a quite narrow region called the basin of attraction. Although this fractal-like shaped area limits the performance of the passive walker from the viewpoint of stability and control, this shortcoming can be remedied by offering some solutions [7,47].
Because of nominally periodic sequence of steps, an effectual tool for studying the impulsive hybrid nonlinear dynamics of walking is based on interpreting a step as a Poincaré map. This map is a stride function [48] relating the state during one part of a step with the state during the same part of the next step, reading From a conventional choice of initial conditions, the above mapping evaluated at the Poincaré section C just after the heel-strike collision. Gait cycles are the fixed points of the linearized Poincaré map which numerically computed from a well-chosen estimate of the initial conditions through the Newton-Raphson iterative scheme. The cyclic stability of these periodic solutions is examined by evaluating the eigenvalues of this linearized map (Jacobian matrix) at the fixed points. The walking stability is preserved if all the eigenvalues of the Jacobian matrix are inside the unit circle. It is known that a general way to identify and analyze stability of fixed points is based on implicitly numerical computation of the linearized Poincaré map, but it will be a difficult task to implement this classic technique successfully for more complex models, especially with high computational cost. In order to bypass these difficulties, proposing new strategies such as shooting method [14], explicitly defined Poincaré map [22,23] and using boundary value problem [46,49] is one of the recent directions for passive bipedal research to discover unforeseen walking patterns.
Here, following our approach in [41,46], a step function (gait cycle) will be created in a numerical process based on the implicitly defined Poincaré map as the conception of a boundary value problem (BVP). This means that the continuous dynamics for a single walking step conceptually forms a BVP in which the boundary conditions correspond to the states at the beginning and end of the step. Indeed, these boundary conditions are related through transition rules at the impact. For a desired one-periodic gait, despite the existence of unanticipated walking patterns, the relevant BVP form of the nonlinear swing equations is reformulated as the following state space model with the state vector x ¼ q _ q ½ T and the step period T, where the initial conditions for a periodic motion Þas well as the time period T are unknown. Applying finite difference method, one of the simplest and fastest ways to discretize highorder equations, it is feasible to convert the nonlinear Fig. 2 Hybrid dynamics of the flexible passive walker with viscoelastic legs differential equations of the swing phase into algebraic equations. In fact, our linearization methodology lies in replacing the derivatives in the continuous dynamics by central finite difference approximations as below The solution of these difference equations for one step in the time interval 0; T ½ is a discretized trajectory consisting of a series of the phase-space points x n t ¼ x t n t ð Þ, which is divided into n t equally subintervals of width Dt ¼ T N . After substituting relations of Eq. (32) in Eq. (31) and linearizing it using the firstorder Taylor approximation, a step-to-step calculation procedure can be obtained for n t time intervals by the following system: where f s is the linearized form of the ordinary firstorder differential equations in Eq. (31). Then, employing the jump relations in Eq. (28) at the end of a single walking step and adding the impact surface expression in Eq. (29) to pinpoint the exact instant of impact contact, our step function makes the following simplified impulsive hybrid dynamics of a walking step: where t ¼ T À and t ¼ T þ refer to the time instances just before and just after impact, respectively, and In a nutshell, the produced multidimensional step function contains the following aspects: 1) nonlinear algebraic equations obtained from finite difference solution (discretization) of continuous dynamics, 2) end-of-step (heel-strike) detection and 3) transition rules. From Eq. (34), it is clear that the appropriate initial conditions x 0 þ and the correspondent time period T are the unknown variables or the zeros of the step function. They can be numerically calculated by using the capabilities of MATLAB software in the Optimization Toolbox (root finding by Fsolve) if a good starting point is supplied. Now, using the favorable initial conditions and the impulsive hybrid dynamics in the state space form, we have sufficient tools to simulate a cyclic walking motion as an Initial Value Problem (IVP) by defining the resulted x 0 as initial values. The capabilities of Matlab software are employed in order to utilize the numerical differential solvers for solving ordinary differential equations (ODEs). The numerical integration of the IVP is performed using the default solver ODE45 with desired accuracy and event function is used to detect zero-crossing of certain algebraic expression in the solver, in this case it notifies when the impact surface in Eq. (29) crosses zero. When ode45 receives this notification, the integration is stopped and the transitions rules in Eq. (28) are applied to calculate the updated positions and velocities of the biped. Hence, as this trend continues, the period-one gait cycle of the flexible passive walker can be simulated.
It is worth to mention that the presented method provides two features: quick calculation and easy implementation in particular for complicated passive models like this flexible biped robot.

Simulation results and discussion
This section will show, through theoretical simulation results, the efficiency of the assumed viscoelastic legs for the simplest passive biped robot. In the first set of simulations, proper initial conditions and then periodone gait cycles are determined by using the suggested method for a flexible passive walker with the default parameter values. Note that the structural parameters of the simplest passive walker can be freely defined, since the biped robot under study in this paper is not based on a physical machine. However, for comparison purposes, the parameters assigned to the walker are similar to those found in [3,41] and are given in Table 1. Also, we will use just one flexible mode for each leg, i.e., n ¼ 1. The next sets of simulations are obtained using the bifurcation diagrams, which demonstrate the effects of the slope angle and the dimensionless mass on the dynamic behavior of the flexible models in Table 1 compared to each other.

Typical walking gait
Before any gait cycle can be analyzed, first the nonlinear algebraic equations defined as a step function in Eq. (34) must be solved iteratively, whose converged solutions depend upon the initial guess. It is important to note that a reasonable choice of this starting guess, meaning within the neighborhood of the desired initial conditions, can provide a suitable solution. In this way, the exact initial conditions resulted for the simplest passive walker with pure elastic legs in our previous work [41] (as seen in Table 1) are used as the initial guess of the step function (Fsolve function in Matlab) of the walker with viscoelastic legs. After this stage, it is time to find walking patterns numerically by application of the IVP (using ODE45 in Matlab) as described in Subsection 2.2 and started with the proper initial conditions from Table 1. Figures 3 and 4 show these simulation results including the period-one gait cycle, the time evolution of angular velocities and the variations of the time-dependent generalized coordinates of the flexible walker with viscoelastic legs. If one can get stable gait cycle from proper initial conditions, there exist fixed point and a set of proper initial conditions around it. The set of proper initial conditions is called basin of attraction. Thus, there is a relation between this area and the initial condition of a Pure elastic g ¼ 0 ½0:0917; 27:38 Â 10 À6 ; 0:1835; À46:68 Â 10 À7 ½À0:1079; 0:00468; À0:0016; À0:00504 limit cycle. Because the basin of the attraction of these bipeds is small, there is a constraint on choosing initial conditions. In other words, the larger of basin of attraction gives more stable walking and can resist against to larger disturbances. By assuming flexible legs, this region becomes smaller, which makes it more difficult to produce stable walking, and thus, in the next level, using active control is necessarily to investigate the basin of attraction. We noted previously in [41] that non-periodic gait cycle was found for the elastic legs with unwanted vibrations. But it is evident from Fig. 3 that using viscoelastic material significantly reduces the vibrations of pure elastic legs and the period-one symmetric steady gait cycle can be obtained. Looking at Fig. 4, it is clear that the flexural deflections of the viscoelastic legs are not symmetric during walking pattern. In fact, the resulted symmetric gaits arise here from asymmetric viscoelastic legs with different deformations. Notice that although the passive walking pattern is simulated for small transverse displacements, we cannot ignore the influence of small deflections on the impulsive hybrid dynamics of the flexible walker, as seen in [41].
The phase portrait diagram and the loci of the eigenvalues of the Jacobian matrix of the step function are also plotted in Fig. 5. The local stability of the closed orbit in Fig. 5a can be analyzed by checking the eigenvalues of the Jacobian at the obtained fixed point.
As anticipated, all of the eigenvalues are in the unit circle as seen in Fig. 5b, and thus, this phase portrait shows convergence to the stable limit cycle of the passive biped robot with viscoelastic legs. Recall that the corresponding Jacobian matrix is calculated by simulating one step for a small perturbation on each of the states of the fixed point.
As explained in [41], we used proportional damping (in the form of Rayleigh damping) to damp the vibrations of the pure elastic legs regardless of its limitations. The results showed that by increasing the damping ratio, the vibrations completely damp out and the flexible biped robot behaves almost like corresponding rigid passive walker for the same slope angle. The obtained initial conditions for a sample of these damped elastic legs with high damping ratio are given in Table 1. Figure 6 reveals the variations of the leg angles for the flexible cases in Table 1 compared to the similar rigid model. As seen, for the first few steps, the leg angles of the flexible models with damping technology are fairly close to the rigid model and the differences grow as time goes on. Nevertheless, the proportional damping model requires high damping ratio for the same performance compared to the viscoelastic material. Our further investigation demonstrates that assuming low damping coefficient (g\0:001) leads to non-periodic gait cycles for the walker with viscoelastic leg. This result is quite obvious, so we refuse to present the corresponding figure. Even in this case, as clearly depicted in Fig. 7, the viscoelastic leg provides a better robustness against vibrations than the similar elastic leg with low damping ratio (f ¼ 10 À4 ) and finally the angular velocities of both flexible legs will be shown fewer oscillations.
Besides the practical limitation of a simple viscous damping (proportional damping), in some special applications such as artificial muscles, the intrinsic viscosity of the structure is the dominant property of the system. Therefore, we are interested in the use of commercially available viscoelastic materials for the compliant legs and also gain insight into the application of passive dynamic walking in building an efficient biped robot.

Analysis with bifurcation diagrams
This part is devoted to examine the effectiveness of the essential parameters: the ground slope c and the dimensionless mass M on the walking motion. By varying these two parameters in the dynamic equations of passive walking, the bifurcation diagrams of the models with viscoelastic and high damped elastic legs as well as the corresponding rigid model can be drawn for certain other parameter values in Table 1. In the present investigation, the bifurcations are observed from the stance leg angle, the angular velocities of both legs and the step period. According to our previous work [41], these diagrams including the initial conditions, the step periods and the walking patterns which are obtained by applying continuation method for analysis the effect of bifurcation parameter. Assuming flexible legs will add considerable complexity to the modeling and consequently to the simulations. Although the concept of bifurcation and chaos for the passive biped robots is a vast field of study, we focused on the variations in the region of period-one gaits for the first theoretical study of flexible passive walkers. Of course, providing further insights on walking behaviors such as multi-periodic and chaotic motions can be considered as a possible future research.
As seen in Garcia's study [3,4], changes in the slope angle give rise to period-doubling bifurcations leading to chaos. Our numerical simulations in Figs. 8, 9, 10 and 11 for the simplest passive walking with two identical viscoelastic, high damped elastic and rigid legs also display this scenario. These figures represent the comparison between the models from the occurrence of the first period-doubling bifurcation and also the region of existing period-one stable gaits. According to these simulation results, using flexible legs decreases the highest slope for stable period-one gaits from 0:0121 (of the rigid model) to c\0:0625. In contrast, both flexible models exhibit the same walking characteristics and there is not much difference between their maximum reachable slope angles for period-one motions because of the high damping ratio of the elastic legs. In other words, the damped elastic legs will be as efficient as the viscoelastic legs if high damping ratio is considered. Generally, it appears that assuming viscoelastic legs expands the region of period-one gait cycles compared to the damped elastic model. In a similar manner with the high damped elastic legs [41], the displayed bifurcation diagram in Fig. 8 confirms that the flexible model with viscoelastic legs follows the walking pattern of the basic rigid model on shallow slope angles. The results in Figs. 9 and 10 show that the angular velocity of the stance viscoelastic leg generates the same behavior as that obtained by means of the rigid stance leg, whereas the angular velocity of the swing viscoelastic leg has such a large difference with the rigid swing leg. Apart from that, the angular velocities of both viscoelastic legs decrease as the slope increases. From Fig. 11, it can be seen that the inclination angle has a small influence on the step period of all models. Moreover, the steps of the passive walker with viscoelastic legs are shorter and faster than the rigid biped.
In the sequel, we investigate the effect of mass ratio M through the bifurcation diagrams in Figs. 12 and 13. As depicted, by increasing the mass ratio, the walking pattern and its properties remain almost constant, but the first bifurcation point is shifted backward along with the increase in slope angle. The main cause of this remarkable change is increasing the leg deflection at the steeper slopes. Despite these identical results for both viscoelastic and elastic legs with high damping ratio, using viscoelastic material still leads to the development of the stable period-one region.
In an effort to improve the accuracy of the assumed modes discretization method, the number of mode shapes used in the viscoelastic model (as defined in Table 1) increases and henceforth we computed the bifurcation diagrams for two assumed mode shapes ðn ¼ 2Þ. The graphs in Figs. 14, 15 and 16 illustrate that by applying the higher mode shapes, the style of walking is preserved and only the first bifurcation point is slightly shifted forward as the slope angle varies. Apparently, this approximate increase in the stable period-one region is even beneficial to the walking stability, in contrast with the other studied flexible passive models in Table 1.
Furthermore, in order to reveal the sensitive dependence of the gait attractors in the phase space on nearby initial conditions, the Maximal Lyapunov Exponent (MLE) is calculated by employing the capabilities of MATLAB software as the slope varies and hence used to distinguish between a variety of steady-state solutions, such as periodic motions and chaotic gaits. It is clear that the positive MLE indicates a chaotic attractor and if it is negative, we have periodic gaits [18]. Notice that the objective of this paper lies mainly on the generation of the period-one passive gaits and the production of the perioddoubling route to chaos for the analysis part. Therefore, the task of further investigation chaos or prediction of the walking dynamics behavior can be performed in future works. However, the computed Maximal Lyapunov Exponents for the studied passive walker with viscoelastic legs in Table 1 are shown in Fig. 17. Comparing with Fig. 14, it can be seen that the region where all these values are negative corresponds to period-one gait and that with positive value correspond to chaotic motions for slopes higher than 0:0625.

Conclusions and future works
In this paper, following our previous work in [41], a further study on the simplest passive biped robot, with a high modeling accuracy, using the structural viscoelasticity (Kelvin-Voigt) effect as a damping mechanism of the flexible leg has been presented. The updated mathematical model of the walking motion that is described by an impulsive hybrid nonlinear dynamics including nonlinear differential equations and a nonlinear algebraic equation has been obtained by the combined Lagrangian-assumed modes method with specific precise boundary conditions. The main focus has been set on finding the proper initial conditions for stable period-one gait cycles that are the fixed points of the new constructed step function and then simulating the walking pattern as an IVP by defining the resulted stable initial conditions as initial values.
By employing the proposed approach, via numerical simulations, we have shown the period-one gait cycle for the simplest passive walker with the Kelvin-Voigt viscoelastic legs. Moreover, the effects of the system parameters containing the mass ratio and the slope angle have been investigated through bifurcation diagrams. Regarding the material properties of the flexible legs (viscoelastic vs. elastic), it was concluded that the viscoelastic legs increased the stability of the model on shallow slopes, even with high mass ratio or two mode shapes compared to the high damped elastic legs. These results convincingly confirm the important role of the viscoelastic leg as an energy dissipater in the flexible passive walker.
In the long term, like new researches in [50,51], studying more complex walking behavior of the introduced flexible passive biped such as Neimark-Sacker Bifurcation and other types of bifurcations will be undertaken. Besides, we aim at increasing the degree of freedom to achieve more realistic biped robot locomotion and considering the influence of the leg parameters such as leg length asymmetry. Meanwhile, using active control such as OGY-based control method, constructing simple experimental instantiations of the model to verify our theoretical results and investigating the effects of smart materials such as piezoelectric energy harvesting would be relevant for our future work.
Data Availability Data sharing is not applicable to this article as no datasets were generated or analyzed during the current study.