Comparing system identification techniques for identifying human-like walking controllers.

While human walking has been well studied, the exact controller is unknown. This paper used human experimental walking data and system identification techniques to infer a human-like controller for a spring-loaded inverted pendulum (SLIP) model. Because the best system identification technique is unknown, three methods were used and compared. First, a linear system was found using ordinary least squares. A second linear system was found that both encoded the linearized SLIP model and matched the first linear system as closely as possible. A third nonlinear system used sparse identification of nonlinear dynamics (SINDY). When directly mapping states from the start to the end of a step, all three methods were accurate, with errors below 10% of the mean experimental values in most cases. When using the controllers in simulation, the errors were significantly higher but remained below 10% for all but one state. Thus, all three system identification methods generated accurate system models. Somewhat surprisingly, the linearized system was the most accurate, followed closely by SINDY. This suggests that nonlinear system identification techniques are not needed when finding a discrete human gait controller, at least for unperturbed walking. It may also suggest that human control of normal, unperturbed walking is approximately linear.

Introduction

While most people can easily walk without conscious thought, it requires appropriate coordination of leg joints, which arises from the central nervous system [1,2] and/or the natural dynamics of the system [3]. Given the redundancy in the human body, gait control is often assumed to follow a hierarchical structure such as the template and anchor paradigm [4]. Even with a simple template model, there are still many possible gaits, each of which can be described by a set of control parameters. The process that dictates how these control parameters are chosen for an imminent step is the controller. Given that the central nervous system is largely a black box, its exact controller is unknown. However, its controller can be inferred by examining the resulting gait. Thus, one objective of this work is to infer a human-like controller given a set of human experimental data.

To accomplish this, the first task is to choose an appropriate template. For walking, there are two common template models—the inverted pendulum model [5] and the spring-loaded inverted pendulum (SLIP) model [6-8]. This work focuses on the SLIP model, partly because previous work has inferred a human-like controller for the inverted pendulum model [9,10]. The SLIP model has two massless, springy legs connected at the hip (figure 1) [6,7]. The hip contains all of the model’s mass. For walking, the model alternates between a single support period when only one leg is on the ground and a double support period when both legs are on the ground. This model has two control parameters, leg stiffness and touchdown angle, both of which need to be controlled appropriately to achieve a desired motion. In an ideal case, a single choice of leg stiffness and touchdown angle produces locally stable, periodic walking. However, human gait is only approximately periodic [11], with small perturbations and variations occurring every step. Thus, the controller must adjust the control parameters every step [12].

Figure 1

Typical SLIP model, shown during the double support period. The plus symbol represents the properties of the new stance leg and the minus symbols represent properties of the previous stance leg.

While a human-like SLIP controller has not yet been identified for walking, many studies have investigated the relationships between the control parameters and walking performance. For humans, leg stiffness tends to increase as walking speed increases, until the transition to running [8,13]. In modelling studies, changing leg stiffness [14] or touchdown angle [7] drastically alters the shape of the ground reaction force (GRF) profile during unperturbed, periodic walking. For perturbed walking, step-to-step leg stiffness adjustment is critical for a stable and efficient gait; for a given biped state, there is a limited stiffness range that allows the next step to be completed [12,15,16]. The range of permitted stiffnesses is primarily determined by walking speed and touchdown angle [7,16]. It is also possible to develop controllers for the SLIP model without considering how human-like the resulting gait is. These studies typically solve the dynamics of the SLIP model, then derive a controller to create a closed system model [17-20]. Beyond this, work with both the inverted pendulum and SLIP models have shown that different controllers can drastically change important gait properties [14,21]. Thus, it is clear that developing an appropriate controller is critical to producing a desired gait.

One way to infer the human controller is to use system identification methods. To do so, both system states (controller input) and observable control parameters (controller output) are measured. These data are then used to find the transformation from the input to the output that best represents the black-box transformation occurring in the central nervous system. While traditionally used for engineering systems, system identification methods have been used to infer neural controllers for both animals [22,23] and humans [9,10,24,25]. Most relevant to this work, system identification has been used to infer a human-like linear controller for a three-dimensional inverted pendulum model [9,10]. This model had two rigid legs connected by a hip joint and transitioned between steps using an instantaneous double-support period. The continuous-time system was converted into a discrete system by taking a Poincaré section at mid-stance. To find the controller, the system states were hip position and velocity at mid-stance, and the control parameters were foot placement and push off impulse. Using an ordinary least-squares method, the human data were converted into a linear system model that directly transformed the mid-stance state to the next step’s mid-stance state. This model was then used to find a linear controller that transformed the mid-stance state into the control parameters. Walking data both with [10] and without [9] external perturbations were used; the results were similar in both cases.

The work in [9,10] treated biped walking as a discrete linear system, which may underestimate nonlinear effects. There has been considerable work showing that gait has nonlinear properties, even when analysed as a discrete system [11,26,27]. Thus, nonlinear system identification techniques may be more appropriate. One method for nonlinear system identification is sparse identification of nonlinear dynamics (SINDY) [28,29]. This method uses a library of nonlinear functions of the system states. It finds a linear combination of nonlinear functions that best models the nonlinear system using sparse identification and experimental data. While SINDY has accurately modelled and controlled many systems, it has not yet been applied to walking. Thus, the second objective of this work was to apply SINDY to walking data and compare the accuracy of SINDY with the simpler and more standard linear system identification techniques used by [9,10].

Therefore, to both (i) infer a human-like controller given a set of human experimental data and (ii) compare different system identification techniques, human walking data were used to identify SLIP-based controllers using three methods. For the first method, a linear Poincaré return map and controller were found using ordinary least squares. A second linear Poincaré return map and controller were found that both encoded the linearized SLIP model and matched the first linear Poincaré return map and controller as closely as possible. A third nonlinear Poincaré return map and controller were found by applying SINDY to the same data. For convenience, the Poincaré return map will be called a model for the remainder of the paper. The results were evaluated using a separate set of human data.

Methods

As just discussed, three methods were used to extract a system model, J, and a controller, K, from human data. The first method performed ordinary least-squares analysis on human walking and leg stiffness data to identify a system model and controller, referred to as the linearized system J and controller K [9]. The second method linearized the SLIP model dynamics and then used the linearized dynamics to optimize a new controller, referred to as the optimized system J and controller K [10]. The third method used SINDY with control and thus used a library of nonlinear functions to identify a linear combination of nonlinear functions that formed the system model and controller, referred to as the SINDY system J and controller K [28,29]. The three methods were then compared with each other.

Human data

Walking data from 12 adult subjects (three male, nine female, 18–63 years of age, mass 69 ± 12 kg, leg length 0.928 ± 0.053 m) were recorded using a VICON motion capture system. Prior to data collection, institutional review board approval from Penn State and informed consent from subjects was obtained. During the experiment, subjects walked on a split-belt instrumented treadmill at 1.1 m s−1 for 6 min, resulting in about 600 steps per subject. Ten per cent of the steps from each subject (approx. 60 steps each or 707 steps total) were held out for evaluation purposes. The remaining steps were used for system identification.

SLIP model

The SLIP model contains a point mass at the hip supported by two legs of different stiffnesses (figure 1). A brief overview of the model is given here; the reader is referred to [6,30] for detailed mathematical derivations. The hip had mass m, which was equal to the mass of the subject. Uncompressed leg length was l0, which was equal to the subject’s leg length. Hip position (x, y) was measured from the hip to the centre of pressure under the foot, both of which were known directly from the measured experimental data. The GRF for each leg was Fleg, which was also measured directly. To convert the measured human data into the model’s leg stiffness k, the effective stiffness at each instant was calculated usingwhere Δl was the stance leg displacement, which was the difference between the subject’s leg length l0 and the distance from the hip to the centre of pressure under the foot. The stiffness was then averaged over the step to find a step stiffness k. The stiffness of each leg was held constant throughout a step and changed during that leg’s swing phase. The touchdown angle αTD was found by calculating the angle between the ground and the line from hip to centre of pressure at heel contact. A step began at touchdown and ended when the foot lifted off the ground. This model can be mathematically described using a set of continuous-time differential equations [30].

To convert the continuous-time SLIP model into a discrete model, a Poincaré section was taken at vertical leg orientation (VLO). VLO occurs when the centre of mass passes over the stance foot. The discrete system’s states were defined aswhere Q is the state vector, the superscript − represents VLO of the current step, and the superscript + represents VLO of the next step. Since x is measured relative to the stance foot, x is always zero at the Poincaré section since it is taken at VLO. Therefore, x is not considered a system state. The control inputs for a step werewhere U− is the control input vector at the current step. To account for subjects of different sizes, Q and U were normalized using subject leg length l0 and mass m. Specifically, hip position was normalized by l0, hip velocity was normalized by , and stiffness was normalized by l0/(mg), where g was the acceleration of gravity. The system states and control inputs were calculated for each step of each subject.

Linearized system identification

Once the system states and inputs were known, the next task was to perform system identification to find both a model and a controller. In all cases, the goal was to find a mapping such that given a hip state Q− and current control input U−, the hip state at the next step Q+ could be predicted. The first method was simple linear system identification. This has been done previously using a stiff-legged three-dimensional model [9]. That model used the same form for Q, but used foot placement and push off impulse for U, instead of stiffness and touchdown angle. Thus, the desired system for this work waswhere and . For each subject, and were the mean of the states and control input, averaged over all their steps [9] (table 1). Thus, J represented the Jacobian of the linearized system. Additionally, a controller was desired such that, given hip state Q− and current control input U−, the control inputs for the next step could be calculatedwhere K was the controller for the linearized system.

Table 1

Unit-normalized state and control input means and standard deviations, averaged across subjects. αTD is in radians. The hip position is (x, y), the leg stiffness is k, and the leg angle at heel contact is αTD.

state/input	mean	std. dev.
x˙	0.321	0.030
y	0.973	0.077
y˙	0.027	0.011
k	28.32	8.45
αTD	1.245	0.104

In order to use the data to find J and K, the Poincaré sections were arranged into matrices such that each column in and represented one step, and the columns of and were the corresponding next step. Using ordinary least squares, J and K were found usingandwhere indicates the Moore–Penrose pseudo-inverse [31]. This calculation used all the training steps for one subject to find J and K such that equations (2.4) and (2.6) held (approximately) true. A different J and K were found for each subject based on their step data. Thus, the linearized model and controller were given by equations (2.4) and (2.6). This was the simplest system identification method, but it did not explicitly include the SLIP model.

Optimized system identification

By contrast, the optimized system identification method explicitly included the SLIP model when finding J and K while still resulting in a final linear system. Similar to the previous method, the Poincaré section was taken at VLO, and the discrete model and controller had the same form, using states at the previous step to anticipate and control the next stepandEquation (2.9) was further broken down into where J = (A + B K). A and B were Jacobians that determined how the states at the next step were affected by small perturbations to previous states and control inputs. Here, all terms were found for each subject separately.

Finding A and B first required finding a periodic walking gait for the SLIP model using optimization. The optimization was performed using the Matlab optimization function fmincon. For each optimization function evaluation, the continuous-time SLIP model was simulated starting at VLO and ending at the next VLO; the objective function was the difference between the beginning and ending states. Both the state vector (giving initial conditions for the SLIP simulation) and the control vector were optimized together to find a periodic gait. The initial conditions for each optimization were and for that subject, resulting in periodic states closest to the subject’s average state. The optimizations were constrained to ensure that the hip did not begin higher than the length of the subject’s leg, and to ensure the control inputs did not fall outside of the ranges used by the subject. Once a periodic gait was found, each state was perturbed and simulated for one step, and the differences in states at the next step were used to calculate the Jacobians A and B.

Once A and B were known, the next task was to find a controller K such thatwas an accurate linear system model [10]. To do this, another optimization was performed to choose K so that the optimized system matched the linear system as closely as possible. Because K appeared in both the model (equation (2.9)) and controller (equation (2.10)), the error in both needed to be considered. The error in the controller was defined as the vector e = K − K for i = 1 : 10. Similarly, the error in the model was defined as the vector e = (A + B K) − J for j = 1 : 15. These vectors were combined into a single error vector e = [e; e]. Then, the Matlab function fmincon was used to solve the optimization problemλerr, was a vector of scaling parameters for each entry of J and K that represented uncertainty, calculated as the inverse of the confidence interval of each entry of K and J. These confidence intervals were found via bootstrap statistics, in which , as well as the corresponding entries in , , and , were sampled randomly with replacement to find equal length vectors of steps, then K and J were calculated again with the sampled steps using equations (2.7) and (2.8). This process was repeated 1000 times, and the confidence interval for each entry of the bootstrap K and J was calculated based on this set of 1000 [32-34]. The solution to equation (2.13) gave K while equation (2.12) gave J. Thus, the optimized model and controller were given by equations (2.9) and (2.10) with the terms calculated using the solutions to equations (2.12) and (2.13).

SINDY system identification

The previous methods linearized the system and the controller. However, human walking is a very complex system, and may not be linear outside of very narrow margins. Thus, finding a nonlinear system and controller may better capture human behaviour. To do so SINDY and SINDY with control (SINDYc) were used [28,29], which required a system of formandwhere and were libraries of nonlinear functions. Thus, the system model (equation (2.14)) depended on the current states and control parameters, as well as the control parameters for the next step. By contrast, the controller (equation (2.15)) only depended on the current states and control parameters. J and K were sparse matrices that specified which functions in and were relevant. System identification required choosing the functions in the library and finding J and K.

The library of functions contained a wide range of nonlinear functions of the states ΔQ− and control parameters ΔU− such as and . A variety of functions were used because it was not clear a priori what all of the important relationships were. From previous work, it is known that horizontal velocity is important [7,8,13,16], so that was included. But since other relationships may also be important, functions not previously related to stiffness or touchdown angle were included as well; for example, horizontal velocity times vertical height, which is related to angular momentum. Because SINDY is well able to balance model accuracy and complexity even with a large library of functions, the functions were chosen systematically knowing that many would not be included in the final model. Thus, the library included all first- and second-order polynomial combinations of the states and control inputs. Additionally, these polynomials were combined with sine and cosine functions of the states and inputs. The trigonometric functions were included mainly because the touchdown angle was a control input and was likely to have a trigonometric relationship with the system. It seemed possible that this may also be true for the other states, so trigonometric functions of all the states and control inputs were included. In total, contained 274 nonlinear functions and contained 131, with the complete libraries given in the electronic supplementary material.

Once the libraries were chosen, the next task was to determine the appropriate weighting for each function. Given a set of data (, , , and ), matrices and were constructed such that every row represented a step and every column represented a different nonlinear function, turning the nonlinear problem into a quasi-linear problem. Thus,andIn many cases, including here, the columns of and must be normalized to ensure that the restricted isometry property holds [35]; this was done by dividing each column of and by the norm of that column. Additionally, after J and K were calculated, each row was divided by the corresponding column norm of and . To solve for J and K, a sequential thresholded least-squares algorithm was used. This algorithm began by solving equations (2.16) and (2.17) using the mldivide function in Matlabandwhich gives the least-squares solution using a QR solver. Then, a sparsity threshold λSINDY was used to find the indices of J and K that were too small, meaning that the corresponding columns of and and the nonlinear functions those columns represent were not important to the system. Those entries were set to zero and the corresponding columns of and were ignored. The process was repeated iteratively with the remaining indices, comparing the values in J and K to λSINDY each time until J or K converged. This generally took less than 10 iterations. Similar to the other methods, a separate J and K were found for each subject. Thus, the SINDY model and controller were given by equations (2.14) and (2.15) with J and K found using the iterative procedure just described.

Each state and control input had a different λSINDY associated with it because the order of magnitude for each state or control input used could be different. Therefore, λSINDY for each column of J and K was found by sweeping through values from 10−4 to 1015, evenly spaced in log coordinates. To determine the appropriate value, the resulting J and K matrices from this sweep were used to predict steps for the held-out data, and the error between the actual and predicted states was calculated for each subject. Section 2.6 describes the procedure used in more detail. In most cases, there was a clear inflection point where higher λSINDY values caused the error to increase substantially (figure 2). The values just before this point were recorded for each subject, and the logarithmic mean over all of the subjects was used as the final λSINDY value (table 2).

Figure 2

Average error in predicting states as a function of the sparsity threshold λSINDY for a representative subject. Increasing λSINDY too much causes errors to increase.

Table 2

Sparsity thresholds used.

associated variable	λSINDY
x˙	7.77 × 103
y	30.6
y˙	9.14 × 104
k	9.58 × 104
αTD	1.70 × 102

Evaluation

In order to evaluate how well the different methods identified a human-like system, the models and controllers were evaluated using the held-out data (§2.1). Several different tests were performed. First, model accuracy was evaluated. For each subject and every test step, the experimental states and control input were used with the linearized model (equation (2.4)), optimized model (equation (2.9)) and SINDY model (equation (2.14)) to predict the states at the next step. For the SINDY model, the output from the controller (equation (2.15)) was used as well because explicitly includes the control inputs at the next step. (J and J incorporate their corresponding controllers implicitly.) The predicted states were then compared with experimental data.

The same test data were used with equations (2.6), (2.10) and (2.15) to evaluate the performance of the identified controllers. The calculated control outputs from each step were then compared with experimental data to determine how well the identified controllers matched the human controller.

Finally, a one-step, continuous-time SLIP simulation was performed for each test step. The experimental states and current control parameters were used as the initial conditions for the simulation. The touchdown angle and new leg stiffness were determined using equation (2.6), (2.10) or (2.15) depending on the model being tested. In some cases, the walker fell over during the simulation; the number of these falls was recorded. For steps that were completed, the simulated and experimental hip states were compared.

Before calculating summary statistics, all of the test results from all subjects were grouped together and outliers were removed using the Matlab function rmoutliers. Outliers were defined using the default threshold of three scaled median absolute deviations. For the linearized system, less than 5% of the test data were removed. For the optimized system, approximately 18% of the test data were removed when testing for accuracy, and approximately 8% were removed when evaluating the simulation results. For SINDY, approximately 8% of the test data were removed when testing for accuracy, and approximately 2% were removed when evaluating the simulation results. The mean (signed) error, the median absolute error, the RMSE and the standard deviation of the signed error were computed for each of the three evaluation tests for each of the three system identification methods for each of the five extended states (, y, , k and αTD). The signed error was approximately normal, motivating the use of a mean while the absolute error did not follow a normal distribution, motivating the use of a median. To determine if there was a consistent bias in the results, Student t-tests were used to check if the average (signed) error of each state was statistically different from zero. To compare system identification methods, paired Wilcoxon signed rank tests were used with the absolute error to determine if one method produced errors that were statistically lower and thus better captured the experimental data. In all cases, α = 0.05.

Results

System identification results

All three identified system models had mean errors of the order of 10−3 in normalized units for all three states although approximately half of the cases had statistically significant non-zero bias (table 3 and figure 3). The median absolute error and RMSE were also of the order of 10−3–10−2. (For comparison, the average and across subjects is given in table 1). This is less than a 10% error for horizontal hip velocity and hip height and a less than 25% error for vertical hip velocity. The large per cent error in vertical hip velocity is primarily due to the fact that experimental vertical hip velocity is very close to zero at VLO. While the optimized system model predicted with the lowest bias as indicated by the lowest signed error, it had the highest values for the other measures. The linearized system model had the lowest median absolute error and RMSE for and was significantly better than the other two systems (p < 0.001). The linearized and SINDY models predicted Δy with similar low median absolute error and RMSE (p = 0.75); both were significantly better than the optimized system model (p ≪ 0.001). The linearized system model predicted significantly better than the other two models (p ≪ 0.001). Therefore, the linearized model had the most accurate results for all three states in general.

Table 3

Mean signed error (μ), median absolute error (|μ|), standard deviation (σ) and RMSE for identified system models when calculating the next state ΔQ+ based on ΔQ− and ΔU−. All values are multiplied by 103.

	Δx˙^+				Δy+				Δy˙^+
system	μ	|μ|	σ	RMSE	μ	|μ|	σ	RMSE	μ	|μ|	σ	RMSE
JL	−0.93	5.58	8.60	8.60	−0.08	1.53	2.40	2.40	−0.52	3.65	5.50	5.50
JO	−0.56	11.4	18.0	18.0	0.35	2.52	4.20	4.20	0.01	4.30	6.40	6.40
JS	−1.40	6.12	9.30	9.40	−0.10	1.50	2.40	2.40	−1.10	4.31	6.50	6.60

Figure 3

Results from testing identified models. In all cases, all three system identification methods were well able to predict the human response. Statistically significant differences between methods are indicated with bars; *** indicates p < 0.001.

The terms in J and J for all subjects are visualized in the electronic supplementary material. The electronic supplementary material also gives J for a representative subject. For the two linear models J and J, high-valued matrix entries correspond to variables that the methods identified as important, while entries that had a consistent sign across subjects indicated that the variable had a consistent effect on gait. For J, non-zero entries correspond to functions that were important; states or control parameters that have many non-zero functions within the library and across subjects were assumed to be important. For all three system identification methods, entries related to were typically larger in magnitude, which was expected since speed is known to be correlated with stiffness. Particularly for the linearized model, the signs of the terms across subjects were consistent. Interestingly, large entries related to Δy were also prevalent in all three models. Touchdown angle did not have large magnitudes, but the signs of the terms tended to be consistent for the linear models. The matrices calculated with SINDY also tended to have larger entries related to the sine and cosine of αTD.

Controller ID results

All of the identified controllers predicted leg stiffness with mean errors below 1 in normalized stiffness units, indicating minimal bias (p > 0.2). The median absolute error and RMSE ranged between 5% and 15% of mean leg stiffness (table 4 and figure 4). The linearized system controller mimicked the human stiffness response the best (p ≪ 0.001) followed by the SINDY controller (p ≪ 0.001). For touchdown angle, the controllers all had mean signed errors of the order of 10−3 radians. Despite the small error magnitude, there was consistent bias for all three controllers (p ≪ 0.001). When considering accuracy, the median absolute error and RMSE were less than 3% of the mean experimental values. Similar to stiffness, the linearized controller predicted the human response best (p < 0.03) followed closely by the SINDY controller. The optimized controller was by far the least accurate (p ≪ 0.001). Therefore, the linearized model had the most accurate controller in general.

Table 4

Mean error (μ), median absolute error (|μ|), standard deviation (σ), and RMSE for identified controllers when calculating ΔU+ based on ΔQ− and ΔU−.

	Δk				ΔαTD (rad)
system	μ	|μ|	σ	RMSE	μ (10−3)	|μ| (10−3)	σ (10−3)	RMSE (10−3)
KL	−0.11	1.40	2.37	2.37	−2.50	7.76	11.6	11.9
KO	0.30	2.53	4.51	4.54	−7.30	13.8	26.0	26.6
KS	−0.08	1.48	2.67	2.67	−2.20	8.01	12.5	12.6

Figure 4

Results from testing identified controllers. Similar to the states, the errors were generally low, indicating good agreement between the experimental and predicted control parameters. Statistically significant differences between methods are indicated with bars; * indicates p < 0.05 and *** indicates p < 0.001.

The terms in K and K for all subjects are visualized in the electronic supplementary material. The electronic supplementary material also gives K for a representative subject. Similar to the model results, entries related to and Δy were typically larger in magnitude for all three models, although the sign of the terms varied by subject. For the linearized controller, the new touchdown angle was also strongly influenced by the previous touchdown angle. For SINDY, αTD had the most non-zero functions although all three states had large numbers of non-zero functions as well.

Simulation results

When the controllers (equations (2.6), (2.10) and (2.15)) were used in conjunction with a continuous-time SLIP model simulation, all three states for all three controllers had significant bias (p ⋘ 0.001, table 5 and figure 5). Horizontal hip velocity at VLO was consistently overestimated by 5–10% on average. Vertical hip velocity at VLO was consistently underestimated; the magnitude of the error was similar to the average experimental value. Regardless of the system identification method, vertical hip velocity had the largest absolute error and largest per cent error of the three states. Vertical hip height at VLO was consistently underestimated by approximately 2% on average. The optimized controller produced steps that matched the experimental states significantly better than both the linearized controller and the SINDY controller (p < 0.001); it had the lowest mean signed error, lowest median absolute error and lowest RMSE for all three states. The linearized controller and SINDY controller performed similarly in simulation (p = 0.01 for Δy and p > 0.45 for and ). However, the SLIP model fell during simulation: 0 times for the linearized controller, 74 times for the optimized controller (including 59 times for one subject) and five times for the SINDY controller (table 6). Thus, the linearized and SINDY controllers appear to be more robust but less accurate than the optimized controller when performing continuous-time simulations.

Table 5

Mean error (μ), median absolute error (|μ|), standard deviation (σ), and RMSE for identified controllers when calculating ΔU+ based on ΔQ− and ΔU−, then simulating for one step using SLIP model to find ΔQ+. All values are multiplied by 103.

	Δx˙^+				Δy+				Δy˙^+
system	μ	|μ|	σ	RMSE	μ	|μ|	σ	RMSE	μ	|μ|	σ	RMSE
JL	21.7	21.1	14.9	26.3	−20.4	19.7	9.4	22.5	−109.2	105	55.5	123
JO	17.4	18.0	15.7	23.4	−16.8	18.4	11.0	20.1	−99.4	90.8	59.1	116
JS	21.5	21.4	15.2	26.3	−20.6	20.1	10.0	22.9	−110.5	105	56.5	124

Figure 5

Results from using identified controllers in simulation. There was significant bias in the predicted states regardless of which system identification method was used. Nevertheless, the errors in and y were acceptable. Statistically significant differences between methods are indicated with bars; ** indicates p < 0.01, *** indicates p < 0.001.

Table 6

Falls recorded during simulations.

	subject number
controller	1	2	3	4	5	6	7	8	9	10	11	12
KL	0	0	0	0	0	0	0	0	0	0	0	0
KO	0	10	0	0	0	0	0	0	59	5	0	0
KS	0	1	1	0	0	0	0	0	3	0	0	0

Discussion

All three system identification methods performed well, with errors below 10% in most cases. Somewhat surprisingly, the linearized system model and controller were the most accurate overall. The SINDY system model and controller were almost, but not quite, as accurate. Given that the linearized model performed best, it appears that nonlinear system identification techniques are not needed when finding a discrete human gait model, at least for unperturbed walking. The optimized controller was generally the least accurate model, although it produced the most accurate results when the controller was used with a continuous-time SLIP model simulation. However, the optimized controller was significantly more likely to choose control parameters that caused a fall in simulation. For the set of data used here, the optimized controller failed completely for one subject and produced significantly more falls for the remaining subjects. Data from these falls were not included in the summary values; if they were, the values for optimized controller would have been much worse.

Regardless of the system identification method, the simulations were far less accurate than directly finding a mapping from one step to the next. When simulated, all three states had significant errors and bias, but the error in was particularly high. In some cases, the hips ended the simulation step moving in the wrong vertical direction. Because the control parameters were generally similar to the human values, but the simulations still performed poorly, this suggests that the SLIP model did not properly capture all human walking dynamics [12]. This may partially explain why the optimized controller was the least accurate method. Because the optimized model directly fit the data to the SLIP model, it was limited to SLIP dynamics. By contrast, the other two methods were informed by the SLIP model but did not directly include it and were therefore able to identify important aspects not present in the SLIP model.

All three models found that horizontal hip velocity was an important factor in predicting at least some of the states and/or control parameters. Previous work using the three-dimensional inverted pendulum models also found that horizontal hip velocity was important [9,10]. Beyond horizontal hip velocity, the system identification methods also seemed to identify Δy as an important factor in walking. This may be because y is directly related to how ‘compressed’ the leg is, which in turn alters the GRF. Δy was not an important factor for the three-dimensional inverted pendulum model in [9]. This difference is not particularly surprising because hip height is primarily a measure of how tilted the leg is for an inverted pendulum model. In other words, hip height captures different dynamical information between the two models, so the impact on the system was expected to be different. Hip vertical velocity did not appear to be important for either the linearized or optimized models, and appeared to be slightly less important than horizontal hip velocity for the SINDY model. This is not entirely surprising because the biped was analysed at VLO when hip height should be close to a maximum for the step, so there should be little vertical velocity. Leg stiffness appeared to be an important factor for the linearized and optimized systems. This was expected because many previous works have shown that leg stiffness alters gait [12,15,16]. Surprisingly, there were very few terms containing stiffness in the SINDY model and controller. Touchdown angle was a consistently important factor for all three system identification methods. While touchdown angle has not been studied as much as stiffness, previous simulation studies have also shown that it is important [7].

Although humans do not completely follow simple continuous time models, their response to normal walking variability appears to be approximately linear. This seems to hold for both the inverted pendulum model [9] and the SLIP model studied in this work. Switching from the very simple linearized system identification method to the nonlinear system identification method SINDY slightly reduced accuracy. This may be because the library of functions used for SINDY did not include some important, but currently unknown, function or measurement, such as angular velocity or angular momentum. It is also possible that SINDY needed more data to find a more accurate model. However, because the linear controller performed well in the first place, it is possible that steady state walking is actually well approximated by a linear approximation of a Poincaré section. Because this work only used unperturbed walking, it is unclear if the linear approximation holds even with noticeable perturbations to walking. It may, since linear system identification found very similar systems for both unperturbed and perturbed walking [9,10]. However, neither of those works compared the linear system with a nonlinear system, and Joshi & Srinivasan [10] noted that there was some indication of nonlinearity. Future work using data that include highly perturbed steps may shed more light on this question. It may also lead to a more robust solution, since the behaviours to correct the gait will be exaggerated and therefore easier for the algorithms to identify.

Nonlinear models may also be required if the control parameters are allowed to vary during a step. The system identification methods used in this work assumed that gait was quasi-periodic and that it could be analysed by taking a Poincaré section. The methods further assumed that stiffness only changed at leg touchdown. This is not entirely realistic because humans can and do modulate their muscles in the middle of a step, which in turn alters the stiffness during the step [36,37]. Thus, a method that identifies a continuous time controller or model may be more appropriate. In particular, such a controller may improve the continuous-time simulation results because it could capture dynamics beyond the basic SLIP dynamics. Finding a continuous time controller for a quasi-periodic system is likely to be prohibitively difficult using linearization techniques, but may be possible through SINDY, which has a continuous time form not used here [29].

In summary, both linear and nonlinear system identification methods can be used to find human-like, discrete-time walking models and controllers based on the SLIP model. At least for the data used here (unperturbed walking from 12 subjects with approx. 600 steps each), a basic linear system identification method worked best, followed closely by the nonlinear system identification method SINDY. Explicitly incorporating the SLIP model into the linear system identification method did not work as well. These results may indicate that the central nervous system uses an approximately linear controller when reacting to small perturbations caused by noise in the neuromuscular system, as opposed to the nonlinear controller that is often assumed.

Click here for additional data file.