Arbitrary Symmetric Running Gait Generation for an Underactuated Biped Model

This paper investigates generating symmetric trajectories for an underactuated biped during the stance phase of running. We use a point mass biped (PMB) model for gait analysis that consists of a prismatic force actuator on a massless leg. The significance of this model is its ability to generate more general and versatile running gaits than the spring-loaded inverted pendulum (SLIP) model, making it more suitable as a template for real robots. The algorithm plans the necessary leg actuator force to cause the robot center of mass to undergo arbitrary trajectories in stance with any arbitrary attack angle and velocity angle. The necessary actuator forces follow from the inverse kinematics and dynamics. Then these calculated forces become the control input to the dynamic model. We compare various center-of-mass trajectories, including a circular arc and polynomials of the degrees 2, 4 and 6. The cost of transport and maximum leg force are calculated for various attack angles and velocity angles. The results show that choosing the velocity angle as small as possible is beneficial, but the angle of attack has an optimum value. We also find a new result: there exist biped running gaits with double-hump ground reaction force profiles which result in less maximum leg force than single-hump profiles.


Introduction
The field of legged locomotion has proposed some fundamental models, including the Spring Loaded Inverted Pendulum (SLIP), the active SLIP, and the Point Mass Biped (PMB). These models allow one to investigate both walking and running gaits. The SLIP model remains a popular tool to investigate bipedal walking and running gaits [1]. This passive model, consisting of a point mass and a massless spring leg, can generate a center-of-mass (COM) trajectory and a ground-reaction-force (GRF) profile similar to that of human running [2]. This has encouraged researchers to investigate SLIP to try to understand the fundamental dynamics of bipedal running [3], and also to use it as a template for controlling real biped robots [4]. For any forward velocity, SLIP can generate periodic running gaits that are passively stable in a narrow region of initial condition parameters, but that are unstable for parameters outside it [5,6]. Therefore, some researchers have proposed stabilizing flight-phase and stance-phase follow SLIP, tested successfully in simulations. The results showed that motor torque saturation remains an obstacle to building practical biped running robots.
One aim of our work is to develop techniques that require less peak force for the stance leg. The success of the template-anchor methods described above motivates the approach in this paper where we use a more general and versatile template than the SLIP model. We propose a novel control law to make an underactuated robot's COM follow any smooth trajectory with any arbitrary initial and final stance conditions (using the PMB model); this result should make it easier to control real biped robots subject to real-world constraints like motor torque saturation. (Note that biped robots with point feet are underactuated in stance phase making them uncontrollable using classic feedback control methods.) Here we compare a circular arc and polynomials of various degrees as COM trajectories in stance phase, but our methods can also be applied to other smooth analytical paths. The variation of cost-of-transport (COT), the maximum GRF with respect to angle of attack, and the COM velocity angle serve as criteria to evaluate the proposed methods.

Biped Models
The benefits of SLIP include its simplicity and its ability to capture some important dynamics of human running [25]. However, it is totally passive. We propose to use PMB in this work in order to capture the advantages of SLIP while avoiding its passivity disadvantage, although we do compare our results to SLIP. PMB consists of a modified SLIP model where a force actuator replaces the spring (Fig 1). The stance leg in both models endures only compressive forces to keep the foot in contact with the ground. Unlike SLIP, PMB can generate arbitrary running gaits.
One step of a running gait includes a stance phase, a take-off, a flight phase and a touchdown. Stance and flight phases are continuous time phases, while take-off and touch-down are instantaneous events. In the stance phase the foot is assumed to be an ideal frictionless passive pivot to the ground. Take-off is an instantaneous transition from stance to flight phase with no discontinuities, occurring when the vertical component of the GRF becomes zero. In flight the robot behaves as a ballistic point mass. Due to the massless legs, there is no impact at touchdown and so no changes occur in the robot's configuration and velocities. The touch-down state becomes the next stance phase initial condition, with no discontinuities.
In stance phase both SLIP and PMB models have 2 DOFs that can be modeled using Cartesian or polar coordinate systems. The body mass is denoted by m, the leg length by r, the leg angle with respect to vertical by θ, the leg axial force in PMB by F r , the spring stiffness in SLIP by K s , the Cartesian unit vectors by i and j, and the polar unit vectors by e r and e θ . We assume stance state vector x s ¼ ½r; y; _ r; _ y. So x s 0 ¼ ½r 0 ; y 0 ; _ r 0 ; _ y 0 constitutes the initial condition of stance. The stance phase equations of motion for SLIP using Cartesian coordinates are ( which constitutes a passive system that utilizes elastic potential energy to boost the mass. On the other hand, PMB ( is an active model without elastic elements and uses only control input F r to boost the mass.
In flight, the ballistic motion of the robot can be described in a straightforward manner using Cartesian coordinates. The massless leg is assumed to be able to move to any desired angle instantaneously with no energy consumption. The flight state vector is x f ¼ ½x; y; _ x; _ y. The dynamic model for running can be expressed as in the equation below, where + andsuperscripts show the instances just after and just before events, respectively.

Running Gait Generation
For the stance phase trajectories we try a circular arc and polynomials of different degrees. Our method generates running gaits and calculates the necessary leg force profiles in stance phase. The proposed method will work for any desired stance phase path profile consistent with the initial condition.

Circular-arc trajectory in stance phase
We desire a symmetric circular-arc trajectory to be generated using the single force actuator of the leg as shown in Fig 2. The center and the radius of the circle are defined using the touchdown angle, leg length and velocity vector direction. The circle is tangent to the initial velocity vector and its center is located at the intersection of the vertical line passing through the foot contact point and the radius line perpendicular to the initial velocity vector. So the stance phase trajectory becomes symmetric with initial and final positions mirrored with respect to the vertical, and initial and final speeds equal to each other. Using a polar coordinate system with the origin at foot contact point, the velocity of COM isṽ with r and θ defined in Fig 1. This equation is stated in Cartesian coordinates as v ¼ ðÀ _ r sin y À r _ y cos yÞĩ þ ð_ r cos y À r _ y sin yÞj: ð5Þ The length of the vertical line OC is The angle between the radial line CH and the vertical line OC is denoted by θ 1 and the radius of the circle is obtained as To calculate the necessary leg force we use a new polar coordinate system with the origin at C, radius of R and angle of θ 1 as shown in Fig 3. The unit vectors of this coordinate system are denoted by e r 1 and e y 1 . The dynamic equations of PMB in the new coordinate system are F r sin ðy 1 þ yÞ À mg sin ( where F r is the leg force and mg is the weight of the robot. Because the radius R of the circle is constant during stance phase, its time derivatives will be zero in Eq (9). The necessary leg force is calculated using the second equation of Eq (9) as According to Fig 3 the angle θ is Substituting Eqs (10) and (11) into the first equation of Eq (9) results in an angular acceleration on the circular trajectory as This nonlinear second order differential equation with respect to θ 1 is solved numerically with the initial condition of to obtain θ 1 and _ y 1 versus time. By substituting this solution into Eq (11) the necessary leg force profile is calculated numerically from Eq (10) to generate the circular-arc trajectory.

Polynomial trajectory in stance phase
As a general framework we define the trajectory of the robot COM to be a polynomial of degree 2n. Since periodic running gaits of SLIP have symmetric stance trajectory, we consider symmetric polynomials with even degrees. The trajectory equation and its time derivatives are written as in which parameters a k can be chosen or designed using COM initial position, velocity and acceleration for stance phase. The Newtonian equations of motion in this phase are written using Eq (2) as Substituting F r from Eq (16), y from Eq (14) and € y from Eq (16) into the first equation of Eq (17), the horizontal displacement differential equation is obtained as The second order nonlinear differential Eq (18) should be solved numerically with the initial condition of ( to obtain x, _ x, € x versus time. Then y versus time is calculated using Eq (14) and substituting those into the first equation of Eq (17) results in the necessary leg force profile to generate the desired polynomial trajectory by COM in stance phase as Verification of the proposed control law To verify the proposed control laws to generate arbitrary trajectories for PMB in stance phase, we apply the calculated force profile F r and the weight to the point mass, with the initial conditions, to generate a trajectory. This trajectory should be identical to the primarily desired circular-arc or polynomial trajectory. To do so, we use equations of motion in polar coordinate system shown in Fig 1b as ( Eq (21) can be solved with known F r and stance initial condition r 0 , θ 0 , _ r 0 , _ y 0 to obtain generated trajectories by the proposed control laws in stance phase.

Cost of transport
The cost of transport (COT) is an energy expenditure index that is used to evaluate biped waking and running efficiency. COT is defined as the energy exhausted by the motors per unit weight of the robot per unit distance traveled. To calculate the exhausted energy of motors, we use the absolute value of the product of the motor force and its displacement.
where W is the exhausted energy of the motors, mg is the total weight of the robot, L is the range of one running step on the ground, t CS is the time duration of one complete step, F r is the leg actuator force, and _ r is the leg length variation rate.

SLIP model
A SLIP with a mass of 61.9 kg, a spring stiffness of 16 kN/m and a spring free length of 0.95 m is considered here to generate the running gait. These values are for the equivalent SLIP model for a typical biped robot named ATRIAS [24]. Starting from an appropriate stance phase initial condition, the SLIP running dynamic Eq (3) can be solved to calculate stance and flight trajectories for one step, as well as the initial condition for stance phase of the next step. A SLIP running Poincare map is composed of dynamic equations of one complete step of running and maps the stance phase initial state vector of one step to the initial stance sate vector of the next step. For a given running velocity, a root finding problem can be solved to find the fixed point of this map which shows a stance phase initial condition that generates a periodic running gait. We solved this problem for the above-mentioned SLIP model with touch-down horizontal and vertical velocity components of V

Circular-arc trajectory for PMB in stance phase
Using the desired initial condition we perform the procedure described in section 3.1 step by step to calculate the necessary leg force to generate a circular-arc trajectory in stance phase. The desired trajectory for stance phase is found using L, R and θ 1 (0) from Eqs (7), (8) and (13).
The necessary force profile for the leg actuator is calculated from Eq (9) and then is substituted into the differential equations in Eq (21). These equations are solved with the desired initial condition to calculate the actual trajectory in stance phase. The desired and resulting COM trajectory are exactly the same (Fig 5). This verifies our gait planning method.
The COM trajectory for one complete running step is shown in Fig 6 in which the thick and thin solid lines stand for stance and flight trajectory respectively, and dash-dot lines show leg potion in touch-down and take-off events. It can be seen that the overall trajectory is similar to a SLIP trajectory. The resulting necessary force profile for the leg actuator is depicted in Fig 7 by a thin solid line and its parameters are shown in Table 1. The leg force starts at its maximum value of 1038 N and reaches its minimum value of 872 N in mid-stance. Note it has a different trend than the SLIP force profile which starts at zero and has maximum value of 1463 N in mid-stance. The maximum required leg force for a circular-arc trajectory is 29% less than with SILP, however it starts from its maximum value. Providing such a large force instantaneously just after touch-down is very difficult in practice.  Polynomial trajectories for PMB in stance phase Polynomials with two parameters. In this section various polynomial trajectories of form Eq (14) are generated for the stance phase of PMB steady running. The first challenge is how to define the polynomial coefficients using the stance phase initial conditions. The stance phase initial state vector x s 0 defines the COM initial position [x 0 , y 0 ] T and velocity ½ _ Substituting these points into Eqs (14) and (15) results in two algebraic linear equations with two unknown polynomial coefficients. So we consider even functions with only two coefficients: y = ax 2 + b, y = ax 4 + b, and y = ax 6 + b. These paths with coefficients calculated from the stance phase initial state vector are shown in Fig 8. Using the procedure of Section 3.2 the required leg actuator force profile for each of these trajectories are calculated from Eq (20). Fig  7 and Table 1 show the resulting force profiles. All of them start and end at their maximum value, and reach their minimum value in mid-stance, similar to a circular-arc path. As the number of the degree of the polynomial increases, the maximum value of force grows and its minimum value decreases. The force profile corresponding to a degree-2 polynomial is very close to a circular trajectory force profile. The maximum force increases with polynomial degree. Because practical robots have motor torque restrictions, providing such large forces at the beginning of stance phase is not feasible. Tuned degree-4 polynomial. To make the leg force profile start and end at zero, one additional constraint is required in the stance phase initial condition. The acceleration of the body just before touch-down is ½€ x 0 ; € y 0 T ¼ ½0; À g T , where g is the acceleration of gravity. To have the GRF start at zero, this acceleration constraint should be added to the stance phase initial condition. So in this case we have three points (position, velocity and acceleration) to be substituted into Eqs (14), (15) and (16). Therefore three coefficients of the path polynomial are  required by these initial conditions, and we choose an even function of degree 4 for the polynomial y = ax 4 + bx 2 + c. The required leg force to make the robot undergo this trajectory during stance is calculated from Eq (20). Then the resulting force profile is substituted into differential Eq (21) to derive the actual stance trajectory. The desired and actual trajectories are verified to be coincident, similar to Fig 5. Fig 9 shows the PMB trajectory for one step of running using this tuned degree-4 polynomial. This trajectory is very similar to SLIP with the same parameters. The stance phase GRF components for both SLIP and PMB with a tuned degree-4 polynomial are shown in Fig 10. It is observed that these two models have almost the same horizontal GRF components but the vertical component of SLIP is sharper, i.e. the PMB model requires less maximum motor torque and less ground friction coefficient than SLIP. These are the main advantages of PMB with a tuned degree-4 polynomial compared to SLIP running. Note also that physical implementation SLIP-based robots, such as ATRIAS [20], have motors in series to the spring. The COM trajectories for one step of running for the all previously generated gaits are shown in Fig 11. All of these gaits have the same touch-down and take-off configuration but different stance trajectories. The flight phases of all these gaits are coincident. In stance phase, it is noticeable that the circular arc and degree-2 polynomial trajectories are very close to each other and the tuned degree-4 polynomial is relatively close to the SLIP trajectory. Horizontal components of COM velocities for these trajectories during stance phase are shown in Fig 7. It is noticeable that again the velocity profiles and the time duration of the circular-arc and degree-2 polynomial trajectories are almost coincident, and the velocity profiles and time duration of SLIP model and degree-4 polynomial trajectories are also close to each other. All of these trajectories start and end with the same velocity of V x = 3 m/s, but their velocities differ slightly during stance phase. Because the start and end points are fixed, this causes different stance times for these trajectories. The SLIP model has a horizontal velocity of 2.83 m/s in . Note all of the generated gaits in this paper use the same initial and final position/velocity for stance phase, and different trajectories differ less than 2% in average velocity.
The corresponding leg force profiles for the various trajectories are shown in Fig 12. It is desirable for the leg actuator force in stance phase to have an initial and final value of zero, and a peak value as small as possible. By this criteria, the tuned degree-4 polynomial has the best force profile specifications among these trajectories, with 11% less maximum force than SLIP. The smallest maximum leg force is found with the degree-2 polynomial, and it is 30% less than SLIP. Although this gait needs less motor force, it has the disadvantage of starting and ending at its maximum value. The COTs of the generated gaits are compared in Table 1. Notice that SLIP is a passive model with a COT of zero. In order to gauge whether the PMB gaits are energy efficient, we compare them to a SLIP model where the leg spring has been replaced by a force actuator that needs to generate the same force profile and displacement i.e. a force-actuated SLIP. The resulting force-actuated SLIP has the best COT with a value of 0.2422, but the PMB gaits have very similar COTs ( Table 1). The worst PMB COTs are found with the degree-2 polynomial and the circular-arc trajectories, both 1.7% higher than force-actuated SLIP. The tuned degree-4 polynomial PMB has a COT only 1.3% higher.
Optimized degree 6 polynomial. We conclude that the PMB tuned degree-4 polynomial trajectory meets the main requirements of biped running. In an effort to discover more about optimal gaits, we next use an even function of degree 6 polynomial y = ax 6 + bx 4 + cx 2 + d, where the coefficient a is a free parameter that satisfies the three initial conditions (similar to the tuned degree-4 polynomial). The relevant trajectories with coefficient a from -30 to 10 are shown in Fig 13, which illustrates that the leg length shortens with increasing a. The corresponding leg force profiles are shown in Fig 14, which depicts single-hump profiles for a > −7.22 and double-hump profiles for a < −7.22. (Note the case a = 0 in Fig 14 is identical to the tuned degree-4 polynomial in the previous section). This is an interesting result because a single-hump force profile is known as a characteristics of bipedal running gaits, while a double hump is characteristic of walking [15,16,25]. The double-hump force profile for a running gait is a new result to the best of our knowledge. It is likely that one could capture running gaits with triple and quadruple-hump force profiles by choosing higher-order polynomial trajectories. We claim that a single-hump force profile is not a characteristics of bipedal running gaits (at least in the field of robotics). Although all observed biologic running gaits demonstrate single-hump force profiles, either single or double-hump profiles can be chosen for the running gaits of robots according to engineering requirements. One important restriction for the running gaits of robots is the maximum torque of the motors. To overcome this restriction we need gaits with less maximum leg force. It can be observed in Fig 14 that there is an optimum value for the coefficient a that minimizes maximum leg force, resulting in a double-hump profile. The variation of maximum leg force with respect to a is shown in Fig 15, which shows that the optimum value of the free parameter is a = −12.24. The leg force profile corresponding to the optimal value of coefficient a is shown in Fig 14 by a dash-dot line.
The important specifications of the force profiles of all previously mentioned stance trajectories are summarized in Table 1. The row for the degree-6 polynomial in this table uses the optimal value a = −12.24 and for mid-stance leg force we use the maximum leg force in the double-hump profile. Notice that stance times differ a bit for the various trajectories because we constrained only stance starting velocities, ending velocities, and its path profile (not average running speed). These differences are not significant and do not cause any problems to our investigation.

Effects of Stance Initial Condition on Running Efficiency
We defined the stance phase state vector as x s ¼ ½r; y; _ r; _ y. Equivalently the velocity can be defined by V x , V y , where the vertical velocity can be written as V y = V x tan β according to Fig  16, β is the angle of velocity vector V with respect to the horizontal, and θ is the leg angle with respect to the vertical. So instead we can use [r, θ, V x , β] as the stance phase state vector and [r 0 , θ 0 , V x0 , β 0 ] as the stance initial condition. In this section we examine the effects of attack angle θ 0 and velocity angle β 0 on the efficiency of the running of a biped robot with free leg length of r 0 and touch-down speed of V x0 . Thus r 0 , V x0 are assumed to be constant with values of 0.95m and 3m/s respectively.
Because the tuned degree-4 polynomial trajectory includes the necessary biped running dynamics and since we aim to investigate trends of gait specifications, we choose the tuned degree-4 polynomial to generate gaits in this section. For each initial condition [r 0 , θ 0 , V x0 , β 0 ] there is a unique tuned degree-4 polynomial. Parameters θ 0 , β 0 are varied from 1˚to 34˚to generate feasible running gaits. Using the proposed control strategy, the stance trajectory is uniquely defined for each initial condition. The COT and the maximum leg force are calculated in these cases, shown in Figs 17 and 18. It turns out that the COT increases with both θ 0 and β 0 , and it approaches zero when both θ 0 and β 0 approach zero. Also it can be seen that maximum leg force decreases with θ 0 and increases with β 0 .
Optimal values of stance initial condition within the intervals are summarized in Table 2. The best COT, with value 0.0176, is obtained by the smallest attack angle and velocity angle (θ 0 = 1˚, β 0 = 1˚), but this produces a relatively large leg force F r,max = 1789.8 N. The smallest maximum leg force is obtained by the smallest velocity angle β 0 = 1˚and the largest attack angle θ 0 = 34˚, but it results in a COT = 0.346, which is a large value relative to the other gaits shown in Table 1. To find an overall optimum gait considering the trade-off between COT and Arbitrary Bipedal Running Gait Generation maximum leg force, we define a normalized index as which is normalized relative to the COT and F r,max values of the tuned degree-4 polynomial   Fig 19. It is worth noticing that all these gaits have the same horizontal velocity but considerably different time duration and step length. This implies that shorter gaits will need a higher frequency of running steps to travel the same horizontal distance in the same amount of time.
Here we put only the normalized sum of COT and F r,max as the objective function to obtain an optimal gait. However, humans, animals and robots have some more constraints which should be taken into account when choosing an appropriate gait. For example choosing a relatively small attack angle θ 0 = 4˚causes an optimal COT and F r,max , but this requires a high frequency of feet switching that would cause muscle fatigue in humans/animals and challenges in robots due to motor response-time restrictions. Thus choosing a larger attack angle makes more physical sense.
Choosing velocity angle β 0 as small as possible is beneficial in minimizing both COT and F r,max . This is convenient because it has been observed that COM travels on an (almost) horizontal line during high-speed, high-efficiency running in both humans and animals.

Conclusions
A point mass biped (PMB) model provided the basis for generating arbitrary symmetric trajectories in the stance phase of steady running. To show the generality of the method, a circulararc trajectory and polynomials of various degrees (all satisfying the initial stance conditions) constituted desired trajectories in stance. The algorithm calculated the required force profile for the leg actuator for each trajectory. We found that paths from polynomials with only two parameters (satisfying initial position and velocity) do not start and end with zero force, which is not desirable. To overcome this problem we added the stance phase initial acceleration to the initial condition, then considered paths with three and four parameters. A degree-4 polynomial with three parameters resulted in a path and force profile very similar to those from a SLIP model, with nearly-equal COT and a maximum leg force 11% less. For a degree-6 polynomial with four parameters, optimizing the maximum leg force resulted in a gait with maximum leg force of 25% less than SLIP. The degree-6 polynomial trajectory can generate biped running gaits with either single-hump or double-hump GRF (previously double-hump profiles was counted as a characteristics of walking [25]). It is worth mentioning that for a defined touch-down velocity, a SLIP model needs a unique attack angle to generate a periodic running gait. However, PMB can generate periodic running gaits with arbitrary touch-down velocities and attack angles (within reasonable limits). We then investigated the effects of attack angle and velocity angle on running efficiency. Choosing the velocity angle as small as possible reduced both COT and maximum leg force. Larger attack angles increased COT and reduced maximum leg force, so there is an optimum value. Using PMB instead of SLIP as a template for multibody biped robots offers some practical advantages: a PMB trajectory can be optimized to overcome real-world limitations such as maximum motor torque, ground coefficient of friction, etc.
We generated symmetric trajectories for the stance phase in this work. Because this was analytic, exact trajectory planning did not deal with stability issues. Implementing these openloop gaits on consecutive running steps of an underactuated robot will generate periodic unstable gaits. Designing controllers to stabilize these planned underactuated gaits (for example using Poincare map control [26] or energy level control [24]) defines future research work. Also, in future work asymmetric trajectories can be investigated in order to reject disturbances in the stance initial condition for steady running, transient running gaits, or running on stairs. Furthermore, this method of gait planning can be applied to 3D PMB and biped robot models with higher degrees of freedom.