Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
IEEE Trans Robot. Author manuscript; available in PMC 2017 August 1.
Published in final edited form as:
PMCID: PMC5222530

Incorporating Human-like Walking Variability in an HZD-Based Bipedal Model


Predictive simulations of human walking could be used to investigate a wide range of questions. Promising moderately complex models have been developed using the robotics control technique hybrid zero dynamics (HZD). Existing simulations of human walking only consider the mean motion, so they cannot be used to investigate fall risk, which is correlated with variability. This work determines how to incorporate human-like variability into an HZD-based healthy human model to generate a more realistic gait. The key challenge is determining how to combine the existing mathematical description of variability with the dynamic model so that the biped is still able to walk without falling. To do so, the commanded motion is augmented with a sinusoidal variability function and a polynomial correction function. The variability function captures the variation in joint angles while the correction function prevents the variability function from growing uncontrollably. The necessity of the correction function and the improvements with a reduction of stance ankle variability are demonstrated via simulations. The variability in temporal measures is shown to be similar to experimental values.

Keywords: Biomimetics, human locomotion, legged robots, motion control, variability

I. Introduction

The ability to model and predict human gait can be used to answer scientific questions that are difficult or impossible to answer with human subject studies. For example, investigations of periodic gait have provided insights into muscle contributions [1]. Further, predictive simulations could allow for preliminary testing of assistive devices, such as validating controllers for a powered leg prosthesis prior to hardware implementation with human amputee subjects [2].

Despite the common assumption that human gait is periodic, significant variability exists. This variability has long-range correlations between steps [3], [4] and structure at the joint level [5]. The joint-level variability has been mathematically described but no attempt to incorporate it into a dynamic model has been made [5]. Increased variability is correlated with increased fall risk [3], but somewhat surprisingly, variability and dynamic stability are not equivalent [6].

Dynamic biped models typically focus on periodic gait, regardless of model complexity (e.g., [1], [7], [8]). The few “human” models that include variability have been relatively simple, with point masses and no knees [9]–[12]. The variability arises from the chaotic nature of the system or from control noise. These simple models have consistently predicted the long-range correlations observed in human gait, although the region of stable model gaits is often much smaller than for human walking. Because fall risk is the most clinically relevant stability metric, in this paper, stability refers to a biped’s ability to avoid falling. Aperiodic walking has received more attention in the field of dynamic bipedal robots, although the focus has typically been on the effects of external perturbations such as ground slope on stability. Various methods of measuring stability and optimizing for robustness have been proposed, including analytic [13] and stochastic [14] measures, BMI optimization [15], and heuristic methods [16]. In all cases, the ideal gait is periodic so the methods cannot be directly used when the desired gait itself is aperiodic. Very little work has been devoted to designing robot gaits in which variability is an essential feature. A notable exception is [17], which developed a systematic method of designing a finite set of both periodic and aperiodic gaits with provable orbital stability.

Investigating the root cause of variability in human walking will likely require extremely complex models that integrate neural control models with muscle models. These individual models are still being developed and combining them for a challenging task such as walking is likely to be very difficult. However, it may be possible to answer many of the more applied scientific questions about variability using a moderately complex model that mimics human variability but does not attempt to explain why it occurs. Unfortunately, such a model does not yet exist. Perhaps the closest existing model [8] is based on the bipedal robot control technique hybrid zero dynamics (HZD), which uses feedback linearization to track a commanded trajectory [18]. This model is currently capable of predicting the average, periodic motion of healthy human walking, but the resulting model simulation has no variability [8]. It is unknown if an HZD-based controller can tolerate human-like variability. Because of both the high feedback gains required and the fact that the instantaneous step-to-step transition tends to amplify errors, introducing even relatively small amounts of variability into the controller may destabilize the system. Nevertheless, physical systems with large deviations from the design model have successfully walked [19], suggesting that the controller may be able to accommodate variability.

This paper adds human-like variability to an existing, HZD-based, moderately-complex, predictive human model [8]. The updated model’s variability is generated by (noise in) the controller, and each step is unique. The mathematical description of the variability is taken from [5], although additional changes to the controller were required to maintain adequate stability. The key theoretical contribution is combining the existing method of predicting mean motion [8] with the existing description of variability [5] in such a way that the biped is consistently able to walk. Sec. II provides an overview of the model and control method. Sec. III describes how the commanded trajectory was modified to generate stable walking with human-like variability, the main contribution. Sec. IV compares the simulation results to human experimental data.

II. Modeling and Control

This section provides a brief overview of HZD-based modeling and control and is included for completeness. The original formulation for point-foot robots was developed in [18] and then modified to capture periodic human walking in [8], [19].

A. Model Overview

The healthy human model developed in [8] was used. The planar six-link model had a point-mass at the hip representing the upper body and two legs with knee and ankle joints (Fig. 1). The parameters were anthropomorphic [20]. To keep model complexity relatively low, the function of the foot and ankle was modeled using a circular foot plus an ankle joint to capture both the center of pressure (CoP) movement [21] and the positive work performed at the stance ankle [22]. Circular feet are commonly used in simple to moderately-complex human biped models (e.g., [7], [8], [10], [23]), and for this model, the ankle joint was needed to capture the other joint-level behaviors and step-level characteristics of human gait [8]. Because the circular model foot captured some of the physiological ankle movement, the motion of the model ankle was expected to be somewhat different than the motion of the physiological ankle. The model ankle variability may also be different than the physiological ankle variability. The biped rolls without slip, so the foot-ground interface was unactuated. This unactuated degree of freedom (DoF) was captured with q1. The remaining joints had ideal actuators that generated torque u2u6.

Fig. 1
Schematic of the model. The model joint angles are indicated with the qi’s. The unactuated angle is q1; the actuated angles are q2q6. The results are analyzed (Fig. 2) using the biomechanics convention indicated with the θi ...

Each step was modeled with a finite-time single support period and an instantaneous double support period during which the stance leg switched. The single support period was modeled with continuous, second-order differential equations:


where x=[qTq.T]TR12 contains the generalized coordinates and their derivatives, and f and g are twelve-dimensional vector-valued functions that can be found from the equations of motion [18]. The instantaneous, impulsive double support period was modeled using an algebraic mapping that related the state of the biped at the instant before impact to the state at the instant after impact [19]:



where SR6 is the resetting map, A is the impact map, and ‘’ and ‘+’ refer to the instants before and after impact.

B. Control Overview

As in the human gait simulations of [8], feedback linearization was used to control the biped during the single support period [18]. The commanded motion of the actuated variables was encoded in output functions to be zeroed [18]:


where h is a five-dimensional vector-valued function to be zeroed, H0R5×6 is a matrix that maps the generalized coordinates to the actuated angles, hc is a five-dimensional vector-valued function of the commanded joint angles, and s is the phase variable. The phase variable parameterizes the commanded motion as a function of hip progression, rather than by time explicitly. Different gaits can be achieved by varying hc.

To drive the output (Eq. 4) to zero, an appropriate controller must be chosen. Differentiating Eq. 4 twice and substituting in the equations of motion (Eq. 1) gives the output dynamics:


To cancel the nonlinearities in the output dynamics, set y¨=v and solve for the input torques:


where v is a stabilizing controller. If the output starts at zero at the beginning of the step, the control law (Eq. 6) ensures that the output remains zero during the single-support period (i.e., continuously invariant). If the gait is perturbed, the stabilizing controller v drives the output to zero. If the output is zero just before impact and is still zero just after impact, then the gait is hybrid invariant. While hybrid invariance is certainly not necessary for stable walking, it does allow the greatest level of control over the biped’s trajectory.

III. Defining the Commanded Motion

Since the goal of this work is to create a model that can be used to study the effects of variability, it is not important to introduce the variability in a physiological manner. Thus, the variability is added to the mean commanded motion. The challenge is in choosing hc such that the simulation matches both the mean motion and the variability in human walking while simultaneously not falling. This section describes the variability refinements of hc and evaluates each refinement’s effect on fall frequency.

The total commanded motion can be parameterized as


where hcM defines the mean motion, hcV defines the variability, and hcC defines a correction polynomial. Using the optimization method from [8], periodic simulation gaits were designed for three of the healthy adult subjects in [5] (a slower subject, a faster subject, and an approximately average subject). Variability was then introduced into the commanded motion, and ten simulations of 100 continuous steps at all speeds were run for each version of hc, although some simulations fell before completing 100 steps. Because it is unlikely humans deliberately corrupt their control signal, the proposed method of adding (known) physiological variability to the optimized mean motion may retain the predictive nature of the original model. The verification of this is left for future work.

A. Defining the Mean Motion

The mean motion is defined using fifth order Bézier polynomials as in [8]. For a periodic gait without variability, the phase variable is normalized to lie between s+ = 0 (impact) and s = 1 (lift-off). When variability is included, the phase variable is unlikely to begin at 0 and end at 1. Unfortunately, the Bézier polynomials quickly diverge to undesirable values outside of the design range [23]. This is particularly problematic at the end of a step when the commanded joint angles frequently throw the swing leg over the hip when s > 1, preventing foot contact and causing the biped to fall forward. To prevent this, the mean motion is switched to a constant velocity profile when s > 1. Thus, the mean motion is parameterized as


where αk are the polynomial coefficients.

B. Defining the Variability

To define the joint angle variability, the method described in [5] and summarized here was used. The variability was defined using a second-order Fourier series for the stance joints and a first-order Fourier series for the swing joints:


where K is the order of the Fourier series, ω is the dominant frequency of the variability, and ak and bk are magnitude coefficients. Based on human subject studies [5], the joint-level variability is approximately sinusoidal, making Fourier series an obvious choice of basis function. Other basis functions could also be used, which would change the details of how variability is represented but would probably not change the overall results because correction polynomials would likely still be needed for this model. Because the only controlled hip joint is the angle between the legs (q2 in Fig. 1), it was treated as a swing joint and stance hip variability was not directly encoded. The magnitude coefficients belong to one of three normal (Gaussian) distributions and the frequencies belong to one of four normal distributions. In addition, there are linear relationships between the magnitude coefficients for the stance hip and knee joints, the magnitude coefficients for the swing hip and knee joints, and the magnitude coefficients for the stance ankle. Finally, there is a nonlinear relationship between the magnitude coefficients for the swing knee and the frequency at the swing knee. The exact method of generating coefficients varied by joint within a general framework [5]. For all joints, the frequency (ω) was chosen randomly from the appropriate normal distribution. For most joints, the bk coefficients were chosen next and used to enforce hybrid invariance in velocity:


Eq. 10 is found by differentiating the output (Eq. 4) with respect to time, substituting in the commanded motion (Eq. 7), setting the output equal to zero, and rearranging. The weighting between randomly choosing bk , satisfying a between-coefficient relationship, and ensuring velocity hybrid invariance was tuned in [5] to ensure human-like variability when only considering kinematics; no additional tuning was done here even though a dynamic model was now being used. Once ω and the bk coefficients were known, the ak coefficients were chosen to enforce hybrid invariance in position:


Similar to the bk coefficients, the weighting found in [5] between randomly choosing ak, satisfying a between-coefficient relationship, and ensuring position hybrid invariance was used.

Initially, hc did not have the correction polynomial hcC so the last term in the hybrid invariance conditions (Eqs. 10 and 11) did not exist. This resulted in a positive feedback loop and uncontrolled growth of the magnitude coefficients. To satisfy Eqs. 10 and 11, at the start of a step, the variation function and its derivative must equal the difference between the actual joint motion and the mean commanded motion. When the difference is large, the magnitude coefficients must also be large, which results in large oscillations and significant differences from the mean commanded motion at the next impact. The impact phase tends to amplify velocity differences, so the next step requires even larger magnitude coefficients to satisfy Eqs. 10 and 11. As a result, defining the commanded motion using just the mean trajectory (hcM) and the variability (hcV) always resulted in the biped simulation falling within ten steps (Table I).

Ability of Simulation to Walk 100 Steps without Falling

C. Adding a Correction Polynomial

To stabilize the variability magnitude, a correction polynomial was added to hc. As is typical in robotic motion planning, a polynomial was used because it is straightforward to ensure smooth control both during the motion and across transitions. The correction polynomial eliminates some of the start-of-step difference between the actual (q) and mean trajectory (hcM) [24]:


The polynomial coefficients ck are found by solving a square linear system to satisfy five continuity constraints:

  • Reduce the start-of-step differences:
    where 0 ≤ λP, λV ≤ 1 are weighting factors defining what percentage of the difference the correction polynomial handles, and
  • Smoothly join the correction polynomial to the mean commanded motion at the mid-point of the step:

λP and λV control how much of the start-of-step difference between the actual motion and the mean motion is zeroed using the correction polynomial vs. the variability function. (Note that Eqs. 13 and 14 are rearranged versions of Eqs. 11 and 10 without the variability terms.) A value of zero means that the variability function must handle all of the difference (recovering version 1 of hc) and a value of one means that the correction polynomial alone ensures hybrid invariance. The choice of where to join the correction polynomial to the mean commanded trajectory only has a small effect on the gait. The choices of λP and λV have a much larger effect on the gait. If they are too small, the variability still grows uncontrollably. If they are too large, there is not enough variability. λP and λV were hand-tuned to keep the magnitude coefficients within their experimental ranges while still allowing appropriate variability (Table II). The tuning was done using two gaits at about the self-selected speed originally developed in [8]. This is version 2 of hc.

Correction Polynomial Constants

The addition of the correction polynomial greatly improved the biped’s robustness, but the biped still fell backwards frequently (Table I). When a step failed, the previous step was slow and ended sooner than expected (s < 1). To a large extent, human gait involves an exchange of potential and kinetic energy with transition periods to redirect the center of mass (CoM) velocity from down to up [7], [25]. The joints perform positive work to replace the energy lost during the transition. In addition, the stance ankle appears to have a large influence on the control of the CoM velocity [25]. The stance ankle generates a significant percentage of its work at the end of the step, so if a step ends too soon, the following step begins with much less kinetic energy than usual. In general, the slower a step, the less kinetic energy it has. If the start-of-step kinetic energy is less than a step’s required increase in potential energy, the biped falls over backwards. This is the failure mechanism for version 2.

D. Reducing Stance Ankle Variability

Because the ankle motion plays a large role in determining the start-of-step kinetic energy and because the model ankle does not correspond to the physiological ankle, the allowable range of the stance ankle magnitude coefficients for the variability (Eq. 9) was reduced by a factor of 12. To ensure that it was actually possible to find coefficients in the reduced stance ankle coefficient range, the stance ankle λP and λV were increased to 0.88 and 0.95, respectively. This final version 3 of the commanded motion function consistently resulted in stable walking (Table I).

Since the model stance ankle variability is reduced compared to the physiological ankle variability, this suggests that humans regulate the foot and ankle function captured by the model ankle fairly precisely. In other words, humans may regulate the control of the CoM velocity closely. Presumably, the CoP movement is regulated less closely and therefore accounts for much of the physiological ankle variability. Alternatively, for this model, the stance ankle is much more sensitive to start-of-step perturbations than the human ankle.

IV. Comparison to Experimental Data

Using version 3 of hc, the 30 simulations were compared to the experimental results. The average speed for the simulations with variability were slower than both the experimental data and the periodic simulation (which matched the experimental speed, Table III). Shorter step lengths tend to correspond to an early impact, which results in a reduction of stance ankle work and tends to slow down the gait. Although the increased impact losses associated with longer step lengths have a slowing effect, longer steps usually correspond to a late impact, which tends to increase stance ankle work and speed up the gait. The net result is a slower than expected average speed. The variability in speed was similar between simulation and experiment. The average simulated stride period was somewhat longer than the experimental stride period, and the variability was higher (Table III). Overall, the simulation did an acceptable job capturing the variability in the step-level temporal parameters.

Temporal Data (Mean ± Standard Deviation)

The mean simulated and experimental hip and knee angles were similar while the ankle angles were significantly different due to the non-physiological ankle-foot model (Fig. 2). Even with the reduction to the variability described in Sec. III-D, the simulated stance ankle had more variability than the experimental data. At the other joints, the simulated joint-level variability was similar to the experimental levels, indicating that the simulation is capable of capturing human-like variability.

Fig. 2
The (a) hip, (b) knee, and (c) ankle angles at a periodic speed of 1.30 m/s from experiment and version 3 (Table I) of the simulation. The results are very similar for the other gaits. The biomechanics angle convention (θi’s in Fig. 1 ...

The variability in the temporal parameters arises from the variability in the joint-level control. Consistent with human experimental results, increased variability tends to result in increased fall risk [3]. The need for the correction polynomials indicates that there is a stabilizing factor in human gait that is not present in the model dynamics. The most likely feature of human walking is the double support period. In human gait, the double support period is both over-actuated and finite time. As a result, the double support period may be used to easily reject destabilizing disturbances since there are redundant actuators. (Some variability is optimal for many biological systems, which may be why humans do not regulate their gait as tightly as possible [3], [4].) In contrast, for the model, the double support period is both unactuated and instantaneous, which tends to magnify differences between the mean and actual motion. As a result, the correction polynomials are needed to artificially reduce the error and prevent failure. These results are consistent with the results for the simple variability models. The models with instantaneous transfers of support [9], [10] were limited to much less variability than the model with a finite-time double support period [12].

V. Conclusion

This paper incorporates human-like joint-level variability into a moderately complex human model by augmenting the mean desired motion with two additional functions – a sinusoidal variability function that randomly injects variability into the system, and a polynomial correction function to help stabilize the system. While this is straightforward from a mathematical perspective, it ignores the fact that human variability likely arises from noise in the neuromuscular control signal. However, even using the formulation presented in this paper (Eq. 7), the input torques (Eq. 6) can be written as the desired mean motion plus a perturbation torque:


where uP is a perturbation term given by


hM is the output function (Eq. 4) containing just the desired mean motion (Eq. 8), and hP is the output function (Eq. 4) containing just the variability functions (Eqs. 9 and 12). Obviously, if the formulation in Eq. 16 is used, a more straightforward equation for uP would be needed, but Eq. 17 should approximate that equation. Eq. 17 clearly shows that the perturbation torque uP is a nonlinear function of both the original control signal and random variability. This is consistent with experimental work showing that there is a correlation between the commanded signal and the amount of noise [26].

This improved human model will be used in preliminary testing of a powered above-knee prosthesis to check that the proposed prosthesis controller does not destabilize the human-prosthesis system. The results may also have implications for optimal bipedal robot control. Since some variability is optimal for human gait [3], [4], incorporating human-like variability into robot gait might improve robustness.


This work was supported by the National Institute of Child Health & Human Development of the NIH under Award Number DP2HD080349. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH. R. D. Gregg holds a Career Award at the Scientific Interface from the Burroughs Wellcome Fund.

Contributor Information

Anne E. Martin, Department of Mechanical and Nuclear Engineering, The Pennsylvania State University, University Park, PA, 16802 USA, ude.usp@43mea.

Robert D. Gregg, Departments of Bioengineering and Mechanical Engineering, University of Texas at Dallas, Dallas, TX, 75080 USA, ude.salladtu@ggergr.


[1] Liu MQ, Anderson FC, Schwartz MH, Delp SL. Muscle contributions to support and progression over a range of walking speeds. J. Biomech. 2008;41(15):3243–52. [PMC free article] [PubMed]
[2] Gregg RD, Lenzi T, Hargrove LJ, Sensinger JW. Virtual constraint control of a powered prosthetic leg: From simulation to experiments with transfemoral amputees. IEEE Trans. Robot. 2014;30(6):1455–71. [PMC free article] [PubMed]
[3] Hausdorff JM. Gait dynamics, fractals and falls: Finding meaning in the stride-to-stride fluctuations of human walking. Hum. Movement Sci. 2007;26(4):555–89. [PMC free article] [PubMed]
[4] Dingwell JB, Cusumano JP. Identifying stride-to-stride control strategies in human treadmill walking. PloS ONE. 2015;10(4):e0124879. [PMC free article] [PubMed]
[5] Martin AE, Villarreal DJ, Gregg RD. Characterizing and modeling the joint-level variability in human walking. Journal of Biomechanics. Accepted. [PMC free article] [PubMed]
[6] Dingwell JB, Cusumano JP, Cavanagh PR. Local dynamic stability versus kinematic variability of continuous overground and treadmill walking. J. Biomech. Eng. 2001;123(1):27–32. [PubMed]
[7] Kuo AD. The six determinants of gait and the inverted pendulum analogy: A dynamic walking perspective. Hum. Movement Sci. 2007;26(4):617–56. [PubMed]
[8] Martin AE, Schmiedeler JP. Predicting human walking gaits with a simple planar model. J. Biomech. 2014;47(6):1416–21. [PubMed]
[9] Garcia M, Chatterjee A, Ruina A, Coleman M. The simplest walking model: Stability, complexity, and scaling. J. Biomech. Eng. 1998;120(2):281–8. [PubMed]
[10] Gates DH, Su JL, Dingwell JB. Possible biomechanical origins of the long-range correlations in stride intervals of walking. Physica A. 2007;380:259–70. [PMC free article] [PubMed]
[11] Roos PE, Dingwell JB. Influence of simulated neuromuscular noise on movement variability and fall risk in a 3D dynamic walking model. J. Biomech. 2010;43(15):2929–35. [PMC free article] [PubMed]
[12] Ahn J, Hogan N. Long-range correlations in stride intervals may emerge from non-chaotic walking dynamics. PLoS ONE. 2013;8(9):e73239. [PMC free article] [PubMed]
[13] Hobbelen DGE, Wisse M. A disturbance rejection measure for limit cycle walkers: The gait sensitivity norm. IEEE Trans. Robot. 2007;23(6):1213–24.
[14] Byl K, Tedrake R. Metastable walking machines. Int. J. Robot. Res. 2009;28(8):1040–64.
[15] Hamed KA, Buss BG, Grizzle JW. Exponentially stabilizing continuous-time controllers for periodic orbits of hybrid systems: Application to bipedal locomotion with ground height variations. Int. J. Robot. Res. To appear.
[16] Post DC, Schmiedeler JP. IEEE/RSJ Int. C. Intell. Robot. Syst. Chicago, IL: 2014. Velocity disturbance rejection for planar bipeds walking with HZD-based control; pp. 4882–7.
[17] Yang T, Westervelt ER, Serrani A, Schmiedeler JP. A framework for the control of stable aperiodic walking in underactuated planar bipeds. Auto. Robot. 2009;27(3):277–90.
[18] Westervelt ER, Grizzle JW, Chevallereau C, Choi JH, Morris B. Feedback Control of Dynamic Bipedal Robot Locomotion. CRC Press; 2007.
[19] Martin AE, Post DC, Schmiedeler JP. Design and experimental implementation of a hybrid zero dynamics-based controller for planar bipeds with curved feet. Int. J. Robot. Res. 2014;33(7):988–1005.
[20] Winter DA. Biomechanics and Motor Control of Human Movement. 4th John Wiley and Sons, Inc.; 2009.
[21] Hansen AH, Childress DS, Knox EH. Roll-over shapes of human locomotor systems: Effects of walking speed. Clin. Biomech. 2004;19(4):407–14. [PubMed]
[22] Anderson FC, Pandy MG. Individual muscle contributions to support in normal walking. Gait Posture. 2003;17(2):159–69. [PubMed]
[23] Quintero D, Martin AE, Gregg RD. IEEE/RAS-EMBS Int. C. Rehabil. Robot. Singapore: 2015. Unifying the gait cycle in the control of a powered prosthetic leg. [PMC free article] [PubMed]
[24] Sreenath K, Park HW, Poulakakis I, Grizzle JW. A compliant hybrid zero dynamics controller for stable, efficient and fast bipedal walking on MABEL. Int. J. Robot. Res. 2011;30(9):1170–93.
[25] Donelan JM, Kram R, Kuo AD. Simultaneous positive and negative external mechanical work in human walking. J. Biomech. 2002;35(1):117–24. [PubMed]
[26] Harris CM, Wolpert DM. Signal-depended noise determines motor planning. Nature. 1998;394:780–4. [PubMed]