Skip to main content
NIHPA Author Manuscripts logoLink to NIHPA Author Manuscripts
. Author manuscript; available in PMC: 2019 Dec 1.
Published in final edited form as: IEEE Trans Neural Syst Rehabil Eng. 2018 Nov 5;26(12):2342–2350. doi: 10.1109/TNSRE.2018.2879570

Modeling the Kinematics of Human Locomotion over Continuously Varying Speeds and Inclines

Kyle R Embry 1, Dario J Villarreal 2, Rebecca L Macaluso 3, Robert D Gregg 4
PMCID: PMC6358210  NIHMSID: NIHMS996046  PMID: 30403633

Abstract

Powered knee and ankle prostheses can perform a limited number of discrete ambulation tasks. This is largely due to their control architecture, which uses a finite-state machine to select among a set of task-specific controllers. A non-switching controller that supports a continuum of tasks is expected to better facilitate normative biomechanics. This paper introduces a predictive model that represents gait kinematics as a continuous function of gait cycle percentage, speed, and incline. The basis model consists of two parts: basis functions that produce kinematic trajectories over the gait cycle, and task functions that smoothly alter the weight of basis functions in response to task. Kinematic data from ten able-bodied subjects walking at twenty-seven combinations of speed and incline generate training and validation data for this data-driven model. Convex optimization accurately fits the model to experimental data. Automated model order reduction improves predictive abilities by capturing only the most important kinematic changes due to walking tasks. Constraints on range of motion and jerk ensure the safety and comfort of the user. This model produces a smooth continuum of trajectories over task, an impossibility for finite-state control algorithms. Random sub-sampling validation indicates basis modeling predicts untrained kinematics more accurately than linear interpolation.

Keywords: Human locomotion, optimization, predictive models, prosthetic limbs, robot control

I. INTRODUCTION

Biological biped locomotion is capable of a diverse array of tasks and feats of agility unmatched by modern prosthetic legs. The goal of powered knee and ankle prostheses is to improve transfemoral amputee quality of life by restoring normative biomechanics for as many locomotion tasks as possible. One example is the Vanderbilt Generation 3 powered knee and ankle prosthesis. This device has been shown to allow amputees to walk with gait kinematics that are similar to the normal gait kinematics of subjects with healthy limbs, and is showing promising results in terms of metabolic performance and minimizing back muscle activation related to the risk of low back pain [1].

Modern powered knee and ankle prostheses require complex control algorithms to provide amputees with intuitive control over their prostheses. A common control strategy, borrowed from gait analysis [2], is to create separate controllers for discrete locomotion activities [3]–[6] and phases of the gait cycle [7]–[9]. A high level controller selects among these controllers using classifier algorithms, which generally depend on input from inertial measurement units (IMU) or electromyography (EMG) sensors [10]–[13]. This finite-state machine (FSM) paradigm requires many controllers to approximate the continuum of all tasks and phases. As the number of controllers increases, so does the time and expertise needed to tune the system to a subject [7]. Therefore a new method capable of supporting a wide range of tasks and gait phases with a single controller is expected to reduce tuning time and promote normative biomechanics.

A different approach to restoring normative biomechanics without the drawbacks inherent to the FSM method is using normative kinematic trajectories as a virtual constraint [14]. Virtual constraints enslave the progression of actuated joints to the progression of a mechanical signal called a phase variable [15]–[18]. Our lab has used healthy human subject kinematic trajectories as the desired trajectories of prosthetic joints and a phase variable derived from the motion of the intact hip joint. This method was shown to improve amputee subjects’ walking performance at discrete speed and incline combinations from 0.67 to 1.21 m/s and −2.5 to 9.0 deg, while reducing the necessary number of tuning parameters and therefore the configuration time [19]. We have also shown promising results that tracking normative trajectories can reduce compensations associated with passive prostheses, including vaulting and hip circumduction [20]. However, this approach is still limited by the use of an FSM to switch between locomotion tasks.

Other researchers in the fields of biped robotics and prosthetics have demonstrated the potential of replacing FSM algorithms with continuous parameterizations. Holgate et al. [21] proposed a mapping from percent gait and stride length to desired actuator position for an ankle prosthesis. This model represents a continuum of phase and walking speed, as stride length is related to walking speed. A more recent approach for controlling powered ankle prostheses modeled a surface representing the time evolution of selected locomotion variables, and regressed that surface over a velocity range to create a manifold [22]. This method creates realistic gait at many speeds in simulation, and provided more power at higher speeds in experiments with an ankle prosthesis. Some biped robots depend on a finite set of optimal gaits, each designed for a specific task, and either define non-periodic transitional gaits to guide the robot back to a preprogrammed periodic gait [23] or interpolate the optimal kinematics for all other tasks in between [24]–[26]. This results in robust gait and the ability to handle a variety of terrains, but requires a large and well-structured set of optimal gaits to interpolate between. The required density of training data is infeasible for human subject experiments.

This paper introduces a method to unify phase and task via a reduced-order data-driven model which will be referred to as the basis model (BM). The joint kinematics over many combinations of speed and incline are parameterized as a function of a phase variable and a novel addition, task variables. This model uses normative data from healthy human subjects performing a variety of walking speeds and inclines to determine optimal model parameters. The goal is to predict kinematic gait patterns for a range of walking speeds and inclines given relatively little training data. Our contributions over prior work [27] include an automated method to select a minimal set of bases (Section II-D3). This step helps to create a simpler model that provides greater model generalization by avoiding overfitting to noisy data [28]. Other new contributions are range of motion constraints on the output, (Section II-D4), minimum jerk regularization of the kinematic surface (Section II-D5), a comparison of this method to another interpolative modeling method, and a more thorough cross-validation process (Section III-B).

The data used to train this model was collected via the procedure given in Section II-A. The model format and training method are described in detail in Section II-C. The predictive accuracy of the basis method, compared to linear interpolation, is discussed in Section II-F, and the numerical results are shown in Section III. Finally, the significance of this work and future plans to utilize the basis model to generate reference trajectories for a transfemoral powered prosthetic leg are discussed in Section V.

II. METHODS

A. Experimental Protocol

The experimental protocol was approved by the Institutional Review Board at the University of Texas at Dallas (UTD). All 10 subjects (5 female) were able-bodied and provided written informed consent. The subject were mean age 23 years (σ = 2.8 years), mean height 170 cm (σ = 8.2 cm), and mean weight 64 kg (σ = 7.7 kg). This selection of subjects represents a population of young, active, and healthy walkers. This should produce the desired kinematics of otherwise healthy amputee community ambulators who can fully utilize the potential capabilities of an advanced powered knee and ankle prosthesis like the 2nd generation UTD LEG [29]. However, it should be noted that the emphasis of this paper is on the solution method, not the specific kinematic data. The only assumptions made about the training data kinematics is that they are periodic and continuously adapt to changes in speed and incline (within the bounds tested). Researchers with different patient populations of interest are encouraged to replicate this optimization method with different training data.

All tests consisted of subjects walking at a steady speed and incline on a Bertec instrumented split-belt treadmill for one minute. For each test the subject walked at a constant speed of 0.8 m/s, 1.0 m/s, or 1.2 m/s and a constant ground slope ranging from −10 degrees to +10 degrees at 2.5 degree increments. All subjects walked at every combination of speed and incline, resulting in 27 different tasks with unique identifiers, χj with j = 1, 2, …, 27. A 10-camera Vicon motion capture system recorded the subjects’ kinematics at 100 Hz. The order of trials was randomized and subjects took breaks to prevent fatigue. A supplementary data set from these experiments is available for download [30].

B. Data Processing

Kinematic data were filtered using Vicon’s Plug-In Gait pipeline [31] which fills gaps and smooths trajectories with a Woltring filter (smoothing factor of ten). The Dynamic Plug-in Gait Model was applied to calculate joint angles [31]. Positive joint angles represent hip flexion, knee flexion, and ankle dorsiflexion throughout this paper.

Kinematic data for each subject were separated into gait cycles and interpolated to contain 150 points per stride. The number of data points was chosen to slightly exceed the number of recorded frames for any stride. All gait cycles begin at heel strike (1) and end just before the next heel strike (150). Heel strike was defined as when the force plate measurement crossed 10% of the subject’s weight, or as the time of maximal anterior position of the heel marker if the subject stepped on the wrong force plate. For each speed and incline combination, the intra-subject mean of each joint angle was calculated across the phase dimension using MATLAB. Outliers, defined as a trajectory that is more than three standard deviations away from the mean trajectory at any point in the stride, were removed. After outlier rejection, the inter-subject mean and standard deviation were calculated across the phase dimension.

C. Gait Model

The basis model represents joint kinematics as a function of phase and task. Phase is measured by a phase variable, φ{|0φ<1,φ˙>0}, which is a cyclic and monotonic scalar that increases from 0 to 1 once per stride. Each independent dimension of the ambulation task is measured by a task variable, which are concatenated into a task vector χ. In our case, χ = (v,α), where v is the subject’s forward speed linearly mapped from a range of 0.6 m/s to 1.4 m/s to a range of 0 to 1, and α is the ground slope, linearly mapped from −10 degrees to 10 degrees to a range of 0 to 1. The methods in this paper can be expanded to include more dimensions of χ, like a sit/stand variable or a flat/stairs variable, but those transitions are out of the scope of this paper.

Gait kinematics are modeled as the weighted summation of N basis functions of phase, bk(φ). The weight of each basis function changes for each unique task, as determined by the task functions ck(χ). This yields the following separable expression for the joint angle q of the hip, knee, or ankle:

q(φ,χ)=k=1Nbk(φ)ck(χ), (1)

where the total number of basis functions is N, indexed by k.

The basis functions bk are finite Fourier series of degree F:

bk(φ)=β00k+i=1F(β1ikcos(iφ)+β2iksin(iφ)), (2)

where the βtik ∈ R terms are coefficients that will be solved for with convex optimization. The t index denotes the type of function (constant, cosine, or sine) that the coefficient pairs with, i indexes the order of the terms, and the k index shows what basis function bk the terms belong to. This general form is widely used in the interpolation of periodic functions. The order F = 10 is selected to ensure our model can match the significant frequency content of human kinematics, as in [32].

The scalar task functions ck describe how the basis function’s weight continuously changes with respect to changing tasks variables. We use Bernstein basis polynomials as task functions because of their frequent use to parameterize geometric shapes on finite intervals [33]. This also means that our resulting model is a Bezier curve, which will produce closed-form solutions for derivatives and guaranteed boundedness. All task functions have the format:

ck(χ)=(γκ)f(χ)κ(1f(χ))γκ, (3)

where (κγ) is a binomial coefficient. The f(χ), γ and κ terms will vary depending on whether a given task function is meant to model changes due to speed, incline, or constant components, as shown by the following sets:

K0={1},K1={2,3,4},K2={5,,8}, (4)

where all task functions ck with k ∈ K0 are constant functions, k ∈ K1 are speed functions, and k ∈ K2 are incline functions. There are three speed terms because we elected to use a 2nd order Bernstein function for speed, the highest order appropriate for only three recorded speeds. A third order Bernstein basis (four terms) was selected for inclines due to the relative complexity of kinematic changes due to incline. The other parameters of (3) take the following values:

kK0f(χ)=0,γ=0,κ=0,
kK1f(χ)=ν,γ=2,κ=0,1,2,
kK2f(χ)=α,γ=3,κ=0,1,2,3.

D. Objectives and Constraints of the Basis Model

Convex optimization is used to solve for the parameters β that determine the shape of the basis functions. The basis model has several objectives and constraints, and each will be discussed separately.

1). Fitting to Training Data:

Our main goal is that the basis model (1) can approximate average human kinematics for a range of tasks. To do this we solve for β in bk that satisfy:

d¯φiχjk=1Nbk(φi)ck(χj), (5)

where d¯φiχj is the inter-subject mean of recorded kinematics of a given joint at discrete phase point φi and discrete task vector χj.

To write a standard form convex optimization problem that will satisfy (5), we define the column vector x:

xk=[β00kβ11kβ21kβ1Fkβ2Fk]T1+2F,
x   =[x1x2xN]TN(1+2F).

A matrix Ai stores the phase dependent terms from (2):

ai=[1cos(φi)sin(φi)cos(Fφi)sin(Fφi)]1×(1+2F),Ai=[ai000ai000ai]N×N(1+2F),

which yields Aixk = bk(φi) from (2). Similarly the task functions from (3) are stored in a row vector:

yj=[c1(χj)c2(χj)cN(χj)]1×N,

which gives the equivalent statements:

yjAix=k=1Nbk(φi)ck(χj).

This matrix format makes it explicit that the basis model is linear with respect to x. Minimizing any norm of d¯φiχjyjAix will constitute a convex objective function, guaranteeing the existence of a global minimum. However, one degree of error has different effects at different points in the gait cycle. To represent this, our objective function will divide the norm of d¯φiχjyjAix by the standard error of the data at that phase point and task, SE(dφiχj). This follows the assumption that points in the gait cycle with high variance do not need to be followed as accurately.

2). Testing against Validation Data:

To avoid overfitting to the specific tasks recorded in our experiment, and to test the feasibility of recording fewer trials in future experiments, we employed a form of repeated random sub-sampling validation called Monte Carlo cross-validation [34]. Cross validation, in general, consists of breaking a data set into two groups, a training set and a validation set. The model parameters are solved in order to approximate the training set, and the model’s performance is tested on the remaining data in the validation set. Monte Carlo cross-validation starts with the data being randomly divided into training and validation sets in z unique ways. Each random instance of cross-validation will be called a test split. In our case, the model parameters are optimized with respect to the training data from tasks χj with j ∈ Hh, where Hh is a subset of the tasks recorded in experiments, and h = 1, 2,… , z indexes through the unique randomized sets of all test splits. The performance of the model is tested against the validation set HhC.

A risk with data-driven models is the potential need for exorbitant amounts of data to produce accurate results. For example, if we added other task variables to our task vector (stairs, running, transitions between) it would become a combinatorial problem to measure a fine grid of every combination of tasks as was done for this experiment. Because of this, it is important that our model can be trained using sparse data in the task space. To test this ability, we created the random training sets Hh using Latin Hypercube Sampling (LHS), a common technique for the design of sparse experiments [35]. LHS discretizes the independent variables of an experiment into n bins, and chooses n points to sample such that each bin of each independent variable is sampled exactly once. To simulate LHS for our current experiment, we take nine unique samples from the available 27 as training data for each test split, such that each incline is represented once and each speed is represented exactly three times. We will generate z = 250 unique test splits with the LHS method to test if our modeling method can make predictions from sparse data.

3). Automated Model Order Reduction:

The basis model has a total of N = 8 basis functions and task functions (4). It would be advantageous to reduce this number if possible. Low rank models which select only the most important model features to explain the response are more likely to provide predictive accuracy by avoiding fitting to model noise, and are simpler to interpret [36]. To facilitate the use of sparsity inducing norms to create lower order models, we introduce a selector function g(K) which returns a vector similar to yj where only task functions with k ∈ K are nonzero:

gk(K)={ck(χj)ifkK,0otherwise,k=1,2,,N. (6)

It is important to note that making ck = 0 also eliminates the contribution of bk towards the full model (1). During the optimization process, it is safe to assume bk = 0 if kK. The g function is used to calculate the constant terms of lower-order basis models, denoted as:

ΛijK=g(K)AiN(1+2F),

which leads to the lower order model:

ΛijKx=kKbk(φi)ck(χj). (7)

We will use convex, sparsity-inducing norms on ΛijKx to reduce the model and use only the most important bases, a set that is defined as K* ⊆ {K0K1K2}.

4). Range of Motion Limits:

The basis model will calculate desired joint trajectories for a prosthetic leg, so bounding the model’s output is a matter of safety for the user and machine. The magnitude of each task function is known before solving for β, and we can rely on this fact to bound the range of motion of each type of task variable (constant, speed, or incline) over a grid of input values indexed by m:

Rimax=Λmmaxim1x+kK1*mmaxΛimkx+kK2*mmaxΛimkx,Rimin=Λmminim1x+kK1*mminΛimkx+kK2*mminΛimkx,i=1,2,,150and m=1,2,,100

where Ki*=KiK* is the set of each type of task function that was determined to be significant, and i indexes over all discrete points of phase. The grid of checkpoints are combinations of 100 speeds and 100 inclines uniformly spanning the domain. The goal is to build a dense enough grid that it is unlikely for the function to leave the range of motion bounds between grid points. Rimax gives the maximum position of the model at discrete phase φi while Rimin gives the minimum. Both can be constrained over the range of i to limit the range of motion at all phase points.

5). Jerk Minimization:

It is important that the output of the basis model resembles natural human motion, and minimizing the jerk of the kinematic trajectories is an effective way to do so [37], [38]. Jerk-minimization will smooth the model output, and is very simple to implement as a convex objective. The jerk at each point in the model can be written as:

JijK=3φ3ΛijK,vec(JKx)=[J11KxJL1KxJ1MKxJLMKx].

Minimizing a norm of vec(JKx) will reduce joint jerk.

E. Solving for Optimal Model Parameters

Our optimization is a two-step process. Step 1 uses all basis functions and reduces the model dimensionality using group L1 regularization [39], as follows:

minimizex,ρ     ρ+λΩ(x), (8)
such thatρSE(dφiχj)d¯φiχjyjAixρSE(dφiχj),Ω(x):=k=1Nmaxi,j|Λijkx|,i=1,2,,150 and j=1,2,,27,

where ρ is a bound on the absolute error between the experimental data d¯φiχj and the model’s prediction yjAix, divided by SE(dφiχj). The function Ω penalizes the number of nonzero basis functions used. The L norm (max of absolute values) is taken of each low order model Λijkx to determine its size and the L1 norm (summation of strictly positive elements inside the Ω function) is taken across these groups. Minimizing the L1 norm encourages sparsity in the residual [40]. In our case, this will encourage a sparse set of basis functions/task functions. The scalar λ is a weighting term, which weights the relative importance of our two objectives: fitting the data accurately (ρ), and using a minimal set of bases (Ω).

The optimization problem (8) has a guaranteed global minimum because it is a convex optimization problem. Each element of d¯φiχjyjAix is affine with respect to vector x. This means the upper bound ρSE(dφiχj) and lower bound ρSE(dφiχj) form affine constraints. The regularizing function Ω is convex, because Λijkx is affine with respect to x, the L norm of an affine expression is always convex, and the sum of convex functions is convex [40]. Similarly, the objective function ρ+γΩ(x) is the sum of an affine expression and a convex function, making it a convex function. As (8) has a convex objective and affine constraints, it is a standard convex optimization problem.

To determine the desired order of the model, we perform a Principle Component Analysis (PCA) of the complete data set. We define the desired order of the data to be the number of eigenvalues greater than three (a threshold that separated the very large and very small eigenvalues) plus one to account for a constant offset (PCA is based on zero mean data, in contrast to our method). This results in a desired order of four for the hip and five for the knee and ankle. We then iteratively weight λ with the bisection method until only the desired number of terms are sufficiently large (> 1).

The second optimization step is smoothly fitting the model to data, using only the basis functions deemed significant by the previous step, denoted K*. This problem is formulated as follows:

minimizex     ρ+δvec(JK*x)2, (9)
such thatρSE(dφiχj)d¯φiχjΛijK*xρSE(dφiχj),R1<Rimin,R2>Rimax,i=1,2,,150 and jHh,

where the term δ is the relative weight between our two objectives, fitting the data accurately (ρ) and smoothing the model (vec(JK*x)2). After some iterations we settled on a value of δ = 1e − 5, but the solution is not highly sensitive to differences in this value. The lower and upper range of motion bounds of each joint, R1 and R2 respectively, are based on values from [41].

The convexity analysis of (9) is similar to (8). We drop the basis regularization objective in favor of a smoothing regularization term vec(JK*x)2. The L2 norm of linear JK*x is convex, so the objective function ρ+δvec(JK*x)2 is convex. The minimization across multiple linear functions is convex, and the positive sum of convex functions is convex, so the lower bound constraint R1<Rimin is valid. Similarly, the maximization of multiple linear functions is concave, and the positive sum of concave functions is concave, making the upper bound R2>Rimax likewise a valid constraint.

F. Determining Predictive Accuracy

It is important to test if the basis model is capable of predicting kinematics for completely new tasks. To quantify predictive accuracy, we calculate the error between the predicted kinematic value, q(φi, χj), and the real value found via experiments, d¯φi,χj, defined as:

Ghj=maxi{d¯φiχjq(φi,χj)SE(dφiχj)},i=1,2,,150,jHhC, and h=1,2,z. (10)

where taking the maximum over i gives the value of the largest error encountered in a gait cycle. This evaluation is repeated for all validation test splits jHhC and stored in a vector of mean errors, eμ, and a vector of maximum errors em:

(eμ)h=meanjGhj,                  (em)h=maxjGhj.

We will compare the mean, max, and distribution of eμ and em between the two methods compared in this paper: basis modeling (BM) as described in (1) and linear interpolation (LN) of kinematics with respect to phase, speed, and incline. Due to the extreme inaccuracy of linear extrapolation of gaits, our linear interpolation model uses linear interpolation for task vectors in the interior of the convex hull of the training data, and an FSM for task vectors on the boundary or outside the convex hull. The FSM uses the nearest training data to predict extrapolated data, as defined by Delaunay Triangulation [42].

The distributions of each test split vector, eμ and em, are compared between the basis modeling method and linear interpolation method using the Wilcoxon signed-rank test. A non-parametric test was used because the Anderson-Darling test rejected the null hypotheses that eμ or em came from normal distributions [43].

III. RESULTS

A. Parameterizing Experimental Data

In this section, we will model hip, knee, and ankle kinematics and evaluate the ability of the basis model to predict untrained task kinematics. To solve this optimization problem, we used CVX, a package for specifying and solving convex programs [44], [45]. Percent gait acts as the phase variable of (1), which will be defined:

φ=tT,φi=tiT,

where t and ti are continuous and discrete time, respectively, from the beginning of a gait cycle, and T is the total time of a gait cycle.

After order reduction, the most important task functions are plotted in Fig. 1. The order of the legend gives the relative importance of each function. Each of these task functions are paired with a basis function, see Fig. 2.

Fig. 1.

Fig. 1.

The most important task functions for each joint, determined by group L1 regularization. The legend gives the order of importance.

Fig. 2.

Fig. 2.

The matching basis functions for each task function in Fig. 1.

B. Evaluating the Model of Experimental Data

Next, we compare the error in our prediction of all untrained tasks to the errors that would result from the use of linear interpolation. Table I shows the mean of the mean summary statistic vector (eμ¯) for all test splits, and the max of the max summary statistic vector (∨em) for all test splits for both the basis modeling method (1) and linear interpolation. Table I shows that the basis model (BM) and the linear interpolation (LN) method have similar hip error means, but the basis model has lower knee and ankle fitting in terms of eμ¯ and ∨em. A one-sided Wilcoxon signed-rank test shows the eμ and em distributions of test split results have lower medians for the basis model on all three joints with the p-values given in the third row of Table I. This shows that the basis model is repeatedly outperforming linear interpolation in terms of accurately predicting hip, knee, and ankle kinematics.

TABLE I.

Summary statistic distribution means, maximums, and Wilcoxon signed rank results among test splits for all joints

Hip Knee Ankle
eμ¯ em eμ¯ em eμ¯ em
BM 0.492 2.10 0.881 3.41 0.856 2.62
LN 0.518 1.59 1.09 5.83 0.994 3.14
p ≪0.05 ≪0.05 ≪0.05 ≪0.05 ≪0.05 ≪0.05

p is the p-value test result of a one sided Wilcoxon signed rank test that the median of a given summary statistic distribution resulting from the BM method is less than the median of the same distribution resulting from the LN method.

Fig. 3 and Fig. 4 show example kinematic surfaces made from the basis method or linear interpolation, respectively. Both use the same random test split for training. In both figures, model outputs are represented as a blue-green surface, and the data used to train the model is shown in dashed blue lines. Solid red lines show the untrained data, which is what was used to calculate validation error. Note that the basis model is C smooth, while the linear interpolation surface has derivative discontinuities at training points and when extrapolating using an FSM.

Fig. 3.

Fig. 3.

Basis model joint surfaces generated for a random test split. The blue dashed mean kinematic trajectories were used for training, solid red is for validation. The true surface is four dimensional, so this is the 3D projection onto one speed, 1.0 m/s.

Fig. 4.

Fig. 4.

Linear interpolation model joint surfaces generated for a the same test split as Fig. 3, projected onto one speed, 1.0 m/s. The dashed blue mean kinematic trajectories were used for training, solid red is for validation.

IV. DISCUSSION

A. Advantages of Basis Modeling

Any powered prosthesis with a position tracking controller will require reference trajectories to follow. There are two basic data-driven approaches to providing those trajectories: an FSM that provides one trajectory for each supported task, or continuous prediction of gait based on the proximity of the current task to known data. The latter has several advantages that we believe are beneficial for emulating natural gait with a prosthesis. Continuous prediction provides a continuum of possible gaits, providing the user with fine control of their gait. Previous work showed continuous prediction over the gait cycle reduces configuration time of a powered leg to individual users [19], and we expect doing the same for both phase and task will have even greater benefits. We have also shown that continuous prediction of gait is more accurate than an FSM at predicting untrained kinematics, even when the FSM is emulated with perfect classification accuracy [27].

There are many options for how to construct a continuous model of human kinematics. The simplest conceptually is linear interpolation, and this has been used in some cases for biped robot gait [46]. The basis modeling method has several advantages over linear interpolation. Through convex optimization, we can force the basis model to fit experimental data while simultaneously minimizing jerk and constraining range of motion. Jerk-minimized kinematic trajectories have been shown to closely mimic natural human movements [37], [38]. The form of the basis model is also C smooth, regardless of solved parameters, see Fig. 3. The basis model also benefits from automated order reduction and a simple closed form expression. Our statistical analysis has shown that basis modeling more accurately predicts untrained kinematics given sparse irregular training data. Performance on sparse data is important because it allows us to design future experiments that will make better models with fewer gait experiments. This aspect will be vital when recording data on impaired populations that cannot be expected to participate in lengthy gait trials. Irregular data is also important as it allows us to utilize experimental design techniques like LHS, or possibly to combine several publicly available data sets of kinematics sampled on different (or irregular) intervals. This second possibility is of interest but is outside the scope of this paper.

B. Interpretation of Trends

Basis and task functions lend themselves well to an intuitive understanding of how task affects gait. The basis function bk(φ) shows what kinematic trend will become more apparent as ck(χ) increases. For example, take the ankle basis b8(φ) in Fig. 2, which contributes additional dorsiflexion around heel strike. As c8(χ) = α3 increases (i.e., as the incline of walking increases from −10 to +10 degrees), the contribution of b8(φ) will increase. Keeping these two terms in mind, it can be deduced that as incline increases, so will ankle dorsiflexion at heel strike. The intuitive behavior of b8(φ)c8(χ) is a well documented feature of inclined walking, also shown in [47, Fig. 2]. A similar example can be seen on the knee joint b5(φ), which will shrink as incline increases because it is paired with knee c5(χ) = −(α − 1)3, see Fig. 1 and 2. This would imply that the knee angle during late stance will be shallower for high inclines, matching again a discrete trend highlighted in [47, Fig. 2]. Part of the power of basis modeling is how these well-established trends emerge naturally from the solution process. The researcher does not have to explicitly define what trends in gait data are important and add terms intended to capture those phenomena.

V. CONCLUSION

This paper develops a method to model human locomotion as a function of phase, speed, and incline. This model will be used as the reference trajectory generator of a virtually constrained powered knee and ankle prosthesis to achieve a continuum of possible gaits over a range of tasks. Various convex optimization techniques are used to improve the accuracy of the model and deliver several useful properties. Model parameters are optimized to match a set of able-bodied human kinematic data recorded at UTD. Structured sparsity inducing optimization methods ensure that the order of the model is minimized, which helps the model to match important trends in the data while maintaining simplicity. Kinematic jerk minimization helps the model output to mimic natural human motion, and preserves the operational life of actuators. Range of motion constraints respect anatomical constraints and will promote safety in future prosthetics applications.

The most important conclusion of this paper is that it is possible to accurately parameterize the human gait cycle as a continuous function of only phase, speed, and incline. This function is represented well by a separable model with basis functions that parameterize kinematics as a function of phase, and task functions that change the relative weight of basis functions depending on the current task. The basis modeling technique was shown to be more accurate than linear interpolation on sparse training data by Monte Carlo cross validation. This paper has shown that, when using the basis modeling method, very fine grids of task data (as collected for this paper) are unnecessary for producing satisfactory interpolation accuracy. This is an important conclusion, because when scaling the model up to a larger task space (e.g., stair climbing) it will be infeasible to sample on the fine grid necessary for accurate linear interpolation. The basis model also provides direct intuition into how human gait is affected by changing tasks.

The next step in this research is developing sensing algorithms for powered prostheses to measure task in real time. This will rely heavily on existing techniques for IMU estimation of walking speed [48]–[50] and incline [51]–[55]. Given measurable tasks and phase variables [56], the following step is to implement this algorithm on the second-generation UTD leg [29].

ACKNOWLEDGMENT

The authors would like to thank Dr. Stephen Boyd for enlightening conversations about the potential of convex optimization to solve for gait basis functions, and Sanyukta Bihari for her contributions to the collection of this experimental data.

This work was supported by the National Institute of Child Health & Human Development of the NIH under Award Numbers DP2HD080349 and R01HD094772. This work was also supported by NSF Award CMMI-1652514. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH or NSF. Robert D. Gregg, IV, Ph.D., holds a Career Award at the Scientific Interface from the Burroughs Wellcome Fund.

Contributor Information

Kyle R. Embry, Department of Mechanical Engineering, University of Texas at Dallas, Richardson, TX 75080, USA..

Dario J. Villarreal, Department of Electrical Engineering at Southern Methodist University, Dallas, TX 75205, USA..

Rebecca L. Macaluso, Department of Bioengineering, University of Texas at Dallas, Richardson, TX 75080, USA.

Robert D. Gregg, Departments of Bioengineering and Mechanical Engineering, University of Texas at Dallas, Richardson, TX 75080, USA..

REFERENCES

  • [1].Jayaraman C, Hoppe-Ludwig S, Deems-Dluhy S, McGuire M, Mummidisetty C, Siegal R, Naef A, Lawson BE, Goldfarb M, Gordon KE et al. , “Impact of powered knee-ankle prosthesis on low back muscle mechanics in transfemoral amputees: A case series,” Frontiers in neuroscience, vol. 12, p. 134, 2018. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • [2].Perry J, Burnfield J, and Cabico L, Gait Analysis: Normal and Pathological Function, 2nd ed. Thorofare, NJ: Slack, 2010. [Google Scholar]
  • [3].Young AJ, Simon AM, Fey NP, and Hargrove LJ, “Intent recognition in a powered lower limb prosthesis using time history information,” Ann. Biomed. Eng, vol. 42, no. 3, pp. 631–641, 2014. [DOI] [PubMed] [Google Scholar]
  • [4].Jiménez-Fabián R and Verlinden O, “Review of control algorithms for robotic ankle systems in lower-limb orthoses, prostheses, and exoskeletons,” Med. Eng. Phys, vol. 34, no. 4, pp. 397–408, 2012. [DOI] [PubMed] [Google Scholar]
  • [5].Tucker MR, Olivier J, Pagel A, Bleuler H, Bouri M, Lambercy O, del R Millán J, Riener R, Vallery H, and Gassert R, “Control strategies for active lower extremity prosthetics and orthotics: a review,” J. Neuroeng. Rehabil, vol. 12, no. 1, p. 1, 2015. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • [6].Lawson BE, Mitchell J, Truex D, Shultz A, Ledoux E, and Gold-farb M, “A robotic leg prosthesis: design, control, and implementation,” IEEE Robotics & Automation Magazine, vol. 21, no. 4, pp. 70–81, 2014. [Google Scholar]
  • [7].Simon AM, Ingraham KA, Fey NP, Finucane SB, Lipschutz RD, Young AJ, and Hargrove LJ, “Configuring a powered knee and ankle prosthesis for transfemoral amputees within five specific ambulation modes,” PloS one, vol. 9, no. 6, p. e99387, 2014. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • [8].Eilenberg MF, Geyer H, and Herr H, “Control of a powered ankle–foot prosthesis based on a neuromuscular model,” IEEE Trans. Neural Syst. Rehabil. Eng, vol. 18, no. 2, pp. 164–173, 2010. [DOI] [PubMed] [Google Scholar]
  • [9].Shultz AH, Lawson BE, and Goldfarb M, “Running with a powered knee and ankle prosthesis,” IEEE Trans. Neural Syst. Rehabil. Eng, vol. 23, no. 3, pp. 403–412, 2015. [DOI] [PubMed] [Google Scholar]
  • [10].Bartlett HL and Goldfarb M, “A phase variable approach for imu-based locomotion activity recognition,” IEEE Trans. Biomed. Eng, 2017. [DOI] [PubMed] [Google Scholar]
  • [11].Stolyarov RM, Burnett G, and Herr H, “Translational motion tracking of leg joints for enhanced prediction of walking tasks,” IEEE Trans. Biomed. Eng, 2017. [DOI] [PubMed] [Google Scholar]
  • [12].Huang H, Zhang F, Hargrove LJ, Dou Z, Rogers DR, and Englehart KB, “Continuous locomotion-mode identification for prosthetic legs based on neuromuscular–mechanical fusion,” IEEE Trans. Biomed. Eng, vol. 58, no. 10, pp. 2867–2875, 2011. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • [13].Young A, Kuiken T, and Hargrove L, “Analysis of using emg and mechanical sensors to enhance intent recognition in powered lower limb prostheses,” Journal of neural engineering, vol. 11, no. 5, p. 056021, 2014. [DOI] [PubMed] [Google Scholar]
  • [14].Gregg RD, Lenzi T, Hargrove LJ, and Sensinger JW, “Virtual constraint control of a powered prosthetic leg: From simulation to experiments with transfemoral amputees,” IEEE Trans. Robot, vol. 30, no. 6, pp. 1455–1471, 2014. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • [15].Westervelt ER, Grizzle JW, Chevallereau C, Choi JH, and Morris B, Feedback Control of Dynamic Bipedal Robot Locomotion. New York, NY: CRC Press, 2007. [Google Scholar]
  • [16].Sreenath K, Park H-W, Poulakakis I, and Grizzle JW, “A compliant hybrid zero dynamics controller for stable, efficient and fast bipedal walking on MABEL,” The International Journal of Robotics Research, vol. 30, no. 9, pp. 1170–1193, 2010. [Google Scholar]
  • [17].Ramezani A, Hurst JW, Hamed KA, and Grizzle JW, “Performance analysis and feedback control of atrias, a three-dimensional bipedal robot,” Journal of Dynamic Systems, Measurement, and Control, vol. 136, no. 2, p. 021012, 2013. [Google Scholar]
  • [18].Powell MJ, Zhao H, and Ames AD, “Motion primitives for human-inspired bipedal robotic locomotion: walking and stair climbing,” in IEEE Int. Conf. Robot. Autom. IEEE, 2012, pp. 543–549. [Google Scholar]
  • [19].Quintero D, Villarreal D, Lambert D, Kapp S, and Gregg R, “Continuous-phase control of a powered knee-ankle prosthesis: Amputee experiments across speeds and inclines,” vol. 34, no. 3, pp. 686–701, 2018. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • [20].Rezazadeh Siavash & Quintero David & Divekar Nikhil & Reznic Emma & Gray Leslie, and D. Gregg Robert. (2018). A Phase Variable Approach for Improved Volitional and Rhythmic Control of a Powered Knee-Ankle Prosthesis. ArXiv e-prints. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • [21].Holgate MA, Sugar TG, and Böhler AW, “A novel control algorithm for wearable robotics using phase plane invariants,” in IEEE Int. Conf. Robot. Autom. IEEE, 2009, pp. 3845–3850. [Google Scholar]
  • [22].Dhir N, Dallali H, Ficanha EM, Ribeiro GA, and Rastgaar M, “Locomotion envelopes for adaptive control of powered ankle prostheses,” in IEEE Int. Conf. Robot. Autom. IEEE, 2018, pp. 1488–1495. [Google Scholar]
  • [23].Liu C, Atkeson CG, and Su J, “Biped walking control using a trajectory library,” Robotica, vol. 31, no. 2, pp. 311–322, 2013. [Google Scholar]
  • [24].Hartley R, Da X, and Grizzle JW, “Stabilization of 3d underactuated biped robots: Using posture adjustment and gait libraries to reject velocity disturbances,” in IEEE Conf. Control Tech. and Applications. IEEE, 2017, pp. 1262–1269. [Google Scholar]
  • [25].Da X and Grizzle J, “Combining trajectory optimization, supervised machine learning, and model structure for mitigating the curse of dimensionality in the control of bipedal robots,” arXiv preprint arXiv:1711.02223, 2017.
  • [26].Nguyen Q, Agrawal A, Da X, Martin WC, Geyer H, Grizzle JW, and Sreenath K, “Dynamic walking on randomly-varying discrete terrain with one-step preview,” in Robotics: Science and Systems (RSS), 2017. [Google Scholar]
  • [27].Embry KR, Villarreal DJ, and Gregg RD, “A unified parameterization of human gait across ambulation modes,” in IEEE Int. Conf. Eng. in Medicine and Biology Society. IEEE, 2016, pp. 2179–2183. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • [28].Bermingham ML, Pong-Wong R, Spiliopoulou A, Hayward C, Rudan I, Campbell H, Wright AF, Wilson JF, Agakov F, Navarro P et al. , “Application of high-dimensional feature selection: evaluation for genomic prediction in man,” Scientific reports, vol. 5, p. 10312, 2015. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • [29].Elery T, Rezazadeh S, Nesler C, Doan J, Zhu H, and Gregg R, “Design and benchtop validation of a powered knee-ankle prosthesis with high- torque, low-impedance actuators,” in IEEE Int. Conf. Robot. Autom. IEEE, 2018, pp. 2788–2795. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • [30].Embry KR, Villarreal DJ, Macaluso RL, and Gregg RD, “The effect of walking incline and speed on human leg kinematics, kinetics, and EMG,” 2018. [Online]. Available: 10.21227/gk32-e868 [DOI]
  • [31].Vicon Motion Systems Ltd. (2018) Vicon nexus reference guide. [Online]. Available: https://docs.vicon.com/display/Nexus26/Vicon+Nexus+Reference+Guide
  • [32].Quintero D, Martin AE, and Gregg RD, “Toward unified control of a powered prosthetic leg: A simulation study,” IEEE Trans. Control Syst. Technol, vol. 26, no. 1, pp. 305–312, 2018. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • [33].Farouki RT, “The Bernstein polynomial basis: A centennial retrospective,” Computer Aided Geometric Design, vol. 29, no. 6, pp. 379–419, 2012. [Google Scholar]
  • [34].Kuhn M and Johnson K, Applied predictive modeling. Springer, 2013, vol. 26. [Google Scholar]
  • [35].Iman RL, Encyclopedia of Quantitative Risk Analysis and Assessment, 2008.
  • [36].James G, Witten D, Hastie T, and Tibshirani R, An introduction to statistical learning. Springer, 2013, vol. 112. [Google Scholar]
  • [37].Lenzi T, Hargrove L, and Sensinger J, “Speed-adaptation mechanism: Robotic prostheses can actively regulate joint torque,” IEEE Robotics & Automation Magazine, vol. 21, no. 4, pp. 94–107, 2014. [Google Scholar]
  • [38].Flash T and Hogan N, “The coordination of arm movements: an experimentally confirmed mathematical model,” Journal of neuroscience, vol. 5, no. 7, pp. 1688–1703, 1985. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • [39].Bach F, Jenatton R, Mairal J, and Obozinski G, “Structured sparsity through convex optimization,” Statistical Science, pp. 450–468, 2012. [Google Scholar]
  • [40].Boyd SP and Vandenberghe L, Convex Optimization. New York, NY: Cambridge University Press, 2004. [Google Scholar]
  • [41].“Range of joint motion evaluation chart,” Washington State Department of Social & Health Services, Report DSHS 13–585A, 2014.
  • [42].MathWorks. (2017) scatteredinterpolant. [Online]. Available: https://www.mathworks.com/help/matlab/ref/scatteredinterpolant.html
  • [43].Whitley E and Ball J, “Statistics review 6: Nonparametric methods,” Critical care, vol. 6, no. 6, p. 509, 2002. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • [44].Grant M and Boyd S, “CVX: Matlab software for disciplined convex programming, version 2.1,” http://cvxr.com/cvx, Mar. 2014.
  • [45].Grant M and Boyd S, “Graph implementations for nonsmooth convex programs,” in Recent Advances in Learning and Control, ser. Lecture Notes in Control and Information Sciences, Blondel V, Boyd S, and Kimura H, Eds. Springer-Verlag Limited, 2008, pp. 95–110, http://stanford.edu/~boyd/graphdcp.html. [Google Scholar]
  • [46].Da X, Harib O, Hartley R, Griffin B, and Grizzle JW, “From 2d design of underactuated bipedal gaits to 3d implementation: Walking with speed tracking,” IEEE Access, vol. 4, pp. 3469–3478, 2016. [Google Scholar]
  • [47].Leroux A, Fung J, and Barbeau H, “Postural adaptation to walking on inclined surfaces: I. normal strategies,” Gait & posture, vol. 15, no. 1, pp. 64–74, 2002. [DOI] [PubMed] [Google Scholar]
  • [48].Yang S and Li Q, “Inertial sensor-based methods in walking speed estimation: A systematic review,” Sensors, vol. 12, no. 5, pp. 6102–6116, 2012. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • [49].Rueterbories J, Spaich EG, Larsen B, and Andersen OK, “Methods for gait event detection and analysis in ambulatory systems,” Medical engineering & physics, vol. 32, no. 6, pp. 545–552, 2010. [DOI] [PubMed] [Google Scholar]
  • [50].Li Q, Young M, Naing V, and Donelan J, “Walking speed estimation using a shank-mounted inertial measurement unit,” Journal of biomechanics, vol. 43, no. 8, pp. 1640–1643, 2010. [DOI] [PubMed] [Google Scholar]
  • [51].Chen B and Wang Q, “Combining human volitional control with intrinsic controller on robotic prosthesis: A case study on adaptive slope walking,” Proc. Annu. Int. Conf. IEEE Eng. Med. Biol. Soc. EMBS, vol. 2015-Novem, pp. 4777–4780, 2015. [DOI] [PubMed] [Google Scholar]
  • [52].Sabatini AM, Martelloni C, Scapellato S, and Cavallo F, “Assessment of walking features from foot inertial sensing,” IEEE Trans. Biomed. Eng, vol. 52, no. 3, pp. 486–494, 2005. [DOI] [PubMed] [Google Scholar]
  • [53].Svensson W and Holmberg U, “An autonomous control system for a prosthetic foot ankle,” in 4th IFAC Symposium on Mechatronic Systems. International Federation of Automatic Control (IFAC), 2006, pp. 856–861. [Google Scholar]
  • [54].Sup F, Varol HA, and Goldfarb M, “Upslope walking with a powered knee and ankle prosthesis: initial results with an amputee subject,” IEEE Transactions on Neural Systems and Rehabilitation Engineering, vol. 19, no. 1, pp. 71–78, 2011. [DOI] [PubMed] [Google Scholar]
  • [55].Rehbinder H and Hu X, “Drift-free attitude estimation for accelerated rigid bodies,” Automatica, vol. 40, no. 4, pp. 653–659, 2004. [Google Scholar]
  • [56].Villarreal DJ, Poonawala HA, and Gregg RD, “A robust parameterization of human gait patterns across phase-shifting perturbations,” IEEE Trans. Neural Syst. Rehabil. Eng, vol. 25, no. 3, pp. 265–278, 2017. [DOI] [PMC free article] [PubMed] [Google Scholar]

RESOURCES