# A Three-Dimensional Mass-Spring Walking Model Could Describe the Ground Reaction Forces

Authors: Huiqi Liang; Wenbo Xie; Zhiqiang Zhang; Peizi Wei; Changhui Cui
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

The oscillatory behavior of the center of mass (CoM) and the ground reaction forces (GRFs) of walking people can be successfully explained by a 2D spring-loaded inverted pendulum (SLIP) model. However, the application of the 2D model is just restricted to a two-dimensional plane as the model fails to take the GRFs in the lateral direction into consideration. In this article, we simulated the gait cycle with a nonlinear dynamic model—a three-dimensional bipedal walking model—that compensated for defects in the 2D model. An experiment was conducted to compare the simulation results with the experimental data, which revealed that the experimental data of the ground reaction forces were in good agreement with the results of numerical simulation. A correlation analysis was also conducted between several initial dynamic parameters of the model. Through an examination of the impact of 3D dynamic parameters on the peaks of GRFs in three directions, we found that the 3D parameters had a major effect on the lateral GRFs. These findings demonstrate that the characteristics of human walking can be interpreted from a simple spring-damper system.


---

## Body

## 1. Introduction

Many mathematical models for dynamic systems in the field of engineering science are formulated as nonlinear differential equations. Many researchers have applied quantity theory or other conservative theories to build differential equations to simulate and analyze nonlinear systems [1–3]. Human walking models derived mathematically are widely used in many engineering fields, such as gait analysis [4–11], human walking assistance [12, 13], motion stability control of robots [14–18], and response prediction of civil structures induced by pedestrians [19–21]. A variety of kinetic models, both passive and active ones, are developed from nonlinear differential equations to study the kinetic phenomenon in human walking. One of the typical passive models that deserves particular attention is the bipedal spring-loaded inverted pendulum (SLIP) model, in which the legs are equivalent to springs, and the GRF is generated through the alternation of the front and rear legs [9, 22–24]. The control strategy to guarantee the stability of walking robots based on the SLIP model has been widely studied and successfully applied [25–31]. Later, a damper was added by Kim et al. to the compliant leg in the classic walking model for a better match with the experimental data, and the result showed that the stiffness and damping ratio exhibited a linear relationship with walking speed [11]. Although the classic bipedal spring-mass model gives a good prediction of the trajectory of the centroid and the complex features of the ground reaction force (GRF), the problem lies in that the previous model overestimates the GRF and underestimates the contact time [32]. Moreover, the SILP model was confirmed to demonstrate the regularity of human walking dynamics concerning age [33]. Lim and Park [34] extended the classic SLIP model to a new one which considered the asymmetry of the ankle joints and proved that the dynamics of multijoint lower limbs during the walking process could be represented by mechanical systems. Li et al. [35] proposed an actuated dissipative spring-mass model that mimicked the actual pattern of human walking. Hao et al. [36] applied a hip torque actuator to the model proposed by Alexander [7], the main advantage of which resides in its reduced energy consumption and a further improvement of the stability. Kuo [37] extended the classic planar motions to a three-dimensional motion which made it possible for the model to tilt side to side, and several strategies were implemented to strengthen the stability of the motion. Roos and Dingwell [38] introduced a dynamic walking model that incorporated a lateral step controller to maintain the lateral stability. They also found that there was a nonlinear relationship between neuromuscular noise and the probability of the model falling. With the undamped inverted pendulum modeled by the Duffing equation as an example, Selvam et al. [39] obtained the existence, Hyers–Ulam stability, and Hyers–Ulam Mittag–Leffler stability of solutions for a discrete fractional Duffing equation with forcing term, and stability conditions for the pendulum were obtained.A major limitation of current studies about the 2D walking model is that all models are restricted to a two-dimensional plane and they fail to include the 3D dynamic parameters [26, 40, 41]. From another perspective, the 2D model is defective as to its inability to calculate the lateral GRF in terms that is the cardinal feature of the oscillation behavior of human walking.dTo compensate for the 2D model’s inability to incorporate lateral GRFs, we present a 3D walking model that embodies GRFs in three directions. The contributions of this paper are as follows: (1) we extended the classic two-dimensional model proposed by Kim and Park [11] to a three-dimensional model that enables a more comprehensive consideration of the GRF of human walking in three directions. (2) We proposed an approach to develop dynamic equations for a 3D walking model through the Lagrangian method. The simulation results showed that our proposed model agreed well with the experimental results. (3) We conducted a correlation analysis between the parameters of the 3D model and confirmed a correlation between some parameters. (4) An investigation was undertaken to study the influence of kinetic parameters on the peak value of GRF.This paper is organized as follows. Section2 gives the kinetic formulation of the 3D walking model and provides a detailed description of the numerical simulation and the experiment protocols. In Section 3, the simulation, experimental results, and numerical analysis are displayed. In Section 4, we discuss the findings and the limitations of our study. Finally, we draw the conclusions in Section 5.


## 2. Materials and Methods

### 2.1. Bipedal 3D Inverted Pendulum Model

To describe the ground reaction forces (GRFs) during the human walking, we proposed a 3D inverted pendulum dynamic model based on the previous 2D model proposed by Kim and Park [11]. In the 3D model, the GRFs exerted by the foot against the ground are obtained from the swing of the center of mass (CoM) in a gait cycle. As shown in Figure 1, we introduce the angles θ1 and θ2 between limbs and the plumb line passing through the center of mass (CoM) to assess the vertical centroid oscillation in a gait cycle. The model consists of two legs with the same stiffness k and damping c, and the lengths of the leading and trailing leg are represented by l1 and l2, respectively. θ11 and θ22 denote the angles between the projection of the two springs on the ground and the forward direction, respectively, to evaluate the lateral motion of the model during the walking process. For simplicity, we refer to the two angles θ11 and θ22 as the 3D angles and θ1 and θ2 as the 2D angles. It is worth noting that the direction of forward movement is self-defined and is fixed over a gait cycle which determines the initial angles and relevant parameters of the model. To construct the 3D model and demonstrate its lateral excursion, the semicircle curve of the classic 2D model is replaced by the hemispheric surface, with GRFs in three coordinate directions taken into account. The hemispheric surface with radius R = 0.3 m is used to simulate the foot (Figure 1). To ensure the transmission of force and moments from the CoM to the foot, the lower end of the spring and the spherical center of the hemispherical surface are rigidly tied. The upper end of the spring is hinged to the CoM, allowing the spring to rotate around the vertical axis freely in the three-dimensional space, with the vertical axis passing through the center of mass.Figure 1
Schematic of the three-dimensional walking model (double-support phase).According to the Lagrange formulation (Appendix A) and space velocity relation (Appendix B), the motion equation of the 3D model for the single-support phase can be obtained as follows:(1)mRsinθ1m0ml12+mR2+2l1mRcosθ1mRsinθ1000msinθ12l12θ..1l..1θ..11=Y1Y2Y3,(2)Y1=−mgcosθ1+kl0−l1+ml1θ˙12+sinθ12θ˙112−cl˙1,(3)Y2=mgsinθ1l1−m2Rcosθ1+l1l˙1θ˙1−sinθ1l1Rθ˙12+cosθ1l1θ˙112,(4)Y3=−2mcosθ1sinθ1l12θ˙1θ˙11−2msinθ12l1θ˙11l˙1.And the motion equation of the 3D model for the double-support phase can be obtained as follows:(5)mRsinθ1m0ml12+mR2+2l1mRcosθ1mRsinθ1000msinθ12l12θ..1l..1θ..11=Y1Y2Y3,(6)Y1=−mgcosθ1+kl0−l1+ml1θ˙12+sinθ12θ˙112−cl˙1+4kl0−l2−cl˙2Rcosθ1+cosθ1cosθ2+cosθ11−θ22sinθ1sinθ2l24Rcosθ2+2+2cos2θ2−cos2θ2−θ22+2cos2θ22−cos2θ2+θ22l2,(7)Y2=mgsinθ1l1−m2Rcosθ1+l1l˙1θ˙1−sinθ1l1Rθ˙12+cosθ1l1θ˙112+4kl0−l2−cl˙2−Rsinθ1l1+Rcosθ11−θ22sinθ2l2−cosθ2sinθ1l1l2+cosθ1cosθ11−θ22sinθ2l1l24Rcosθ2+2+2cos2θ2−cos2θ2−θ22+2cos2θ22−cos2θ2+θ22l2,(8)Y3=−2mcosθ1sinθ1l12θ˙1θ˙11−2msinθ12l1θ˙11l˙1+4kl0−l2−cl˙2sinθ1sinθ2sinθ11−θ22l1l24Rcosθ2+2+2cos2θ2−cos2θ2−θ22+2cos2θ22−cos2θ2+θ22l2,where m,R,θi,li,θii,θ˙i,l˙i,θ˙ii,l0,k, and c represent the model mass, the radius of the foot hemisphere, 2D spring angle, spring length, the 3D angles, the 2D angular velocity, the change rate of spring length, the 3D angular velocity, slack length of spring, spring stiffness, and damping constant, respectively. The subscripts i=1,2 define the trailing and leading leg (Figure 1). It is worth noting that the three mutually independent variables θ1,l1, and θ11 can determine the double-support phase of the model due to the existence of geometric and the dynamic constraint.


### 2.2. Simulation Method

As shown in Figure1, the positive direction of the x-axis is set as the initial direction of walking. With a similar pattern to that of the 2D model in [11, 22], the gait cycle in the 3D model is initiated from a single-support phase in which the 6 initial state parameters need to be determined. The parameters are θ1,θ˙1,θ11,θ˙11,l1, and l˙1, the details of the parameters have been described before. During the simulation for the integration of the equations, the 3D model is initially positioned in a state when the point mass is at its apex during a gait cycle. The forward limb is positioned in front of the trailing limb at a predetermined angle θtd. Before the contact of the leading leg with the ground, the motion equations (equation (1)) of the single-support phase are used to determine the dynamical behaviors of the model. It is worth noting that the angle θ1 of the 3D model is not initially zero, but a small value approaches zero, which is different from that in the 2D model to suit the need for convergence of the differential equation. In the numerical simulation, it is found that a stable result could be obtained when θ1 is in the range of 0.02 rad to 0.04 rad. Meanwhile, the initial angle θ11 and the corresponding angular velocity θ˙11 are selected through the Genetic Algorithm (GA) which will be described below.When the model is transformed from the single-support phase to the double-support phase, the following geometric constraint needs to be satisfied:(9)l1cosθ1=l2cosθ2.Especially when the leading leg is just in contact with the ground, the following condition needs to be satisfied:(10)l1cosθ1=l0cosθtd,where θtd represents the touch-down angle. The calculation of the dynamic states of the model in the double-support phase requires the closed kinematic constraints (equation (B.4) in Appendix B) to determine the parameters related to the limb spring velocity l˙2, 2D angular velocity θ.2, and 3D angular velocity θ˙22. If the trailing leg is lifted, and the point mass is at its apex, the trailing leg will be reset in front of the body with the predetermined impact angle. The model is instantaneously rotated along the plumb line passing through the body center, and the 3D angle and angular velocity after rotation are identical in value to those before the rotation but in the opposite direction. The process corresponds to the phase of “instant switch of legs” in Figure 2.Figure 2
Perspective (a), top view (b), and side view (c) of the walking model in the complete gait cycle.
(a)(b)(c)During a gait cycle, it is assumed that the change in the CoM and gait is symmetrical, which is consistent with the symmetric features of human walking. Based on the symmetric features and the consideration of computational efficiency, half of the gait cycle of the 3D model is calculated which is expressed as the single-double-single process in Figure2. By iterating the initial conditions to be consistent with the end conditions, a stable GRFs curve can be finally obtained. The vertical GRF deserves more attention than the transverse and longitudinal components of the GRF for the vertical GRF is assumed to better reflect the oscillation characteristics of human walking. To optimize the initial parameter values, the Genetic Algorithm (GA) Solver of MATLAB® is utilized to identify the optimal solution for a specified nonlinear system with the fitness function set as the root mean square error between the data of the simulation of vertical GRF and the experimental ones. To prevent blind selection of parameters from causing nonconvergence of the solution, the parameter range is obtained from trial calculations by reference to the results of [22], and the final parameter range is as shown in Table 1. The parameters of GA are given as follows: the population size is set to 100, the mutation rate to 0.01, and the maximum number of generations for each run to 50.Table 1
Normalized model parameters for genetic algorithm calculations.
Global model parameters:M=80kg,l0=0.8m, andR=0.3mLeg stiffness,kN/m17500∼22500Damping ratio,ξ0.02∼0.05Initial leg angle,θ1rad0.02∼0.04Initial leg angular velocity,θ˙1rad/s1.25∼2.00Initial3D leg angle,θ11rad0.3∼0.6Initial3D leg angular velocity,θ˙11rad/s−1.20∼−0.40Touch−down angle,θ2degree18∼24


### 2.3. Experiment

To validate the model proposed for 3D walking, the experimental data were collected to compare with the simulation results. Thirteen young and healthy subjects (body mass of 67.4 ± 12.2 kg and height of 1.67 ± 0.12 m) with no gait disorders participated in the experiment. The subjects walked on a walkway equipped with three force plates with piezoelectric 3-component force sensors. The participants were instructed to keep pace with a metronome placed nearby, which provided a fixed frequency each time. During the experiment, each participant was guided to walk at the frequency of the metronome, three times at each frequency. To get the ground reaction force (GRF) of each participant, three six-dimensional force plates (JP6060, CN) were used to measure the real-time GRF at a sample rate of 100 Hz. To ensure that the participants would cross the force plate at a comfortable and stable pace, the pre- and post-transition parts of the walkway made of wood were placed before and after the force plate (Figure3). An AHRS system (AH100B, CN) was used to measure the acceleration and direction angle of CoM during the walking process at a sample rate of 100 Hz, as shown in Figure 4. The AHRS system could accurately and stably perform real-time 3D dynamic measurements with a Kalman filter. The velocity of the CoM could be calculated by integrating the acceleration data collected by the AHRS system.Figure 3
Test site layout of the walkway.Figure 4
Participant walks on the walkway at a fixed frequency.


## 2.1. Bipedal 3D Inverted Pendulum Model

To describe the ground reaction forces (GRFs) during the human walking, we proposed a 3D inverted pendulum dynamic model based on the previous 2D model proposed by Kim and Park [11]. In the 3D model, the GRFs exerted by the foot against the ground are obtained from the swing of the center of mass (CoM) in a gait cycle. As shown in Figure 1, we introduce the angles θ1 and θ2 between limbs and the plumb line passing through the center of mass (CoM) to assess the vertical centroid oscillation in a gait cycle. The model consists of two legs with the same stiffness k and damping c, and the lengths of the leading and trailing leg are represented by l1 and l2, respectively. θ11 and θ22 denote the angles between the projection of the two springs on the ground and the forward direction, respectively, to evaluate the lateral motion of the model during the walking process. For simplicity, we refer to the two angles θ11 and θ22 as the 3D angles and θ1 and θ2 as the 2D angles. It is worth noting that the direction of forward movement is self-defined and is fixed over a gait cycle which determines the initial angles and relevant parameters of the model. To construct the 3D model and demonstrate its lateral excursion, the semicircle curve of the classic 2D model is replaced by the hemispheric surface, with GRFs in three coordinate directions taken into account. The hemispheric surface with radius R = 0.3 m is used to simulate the foot (Figure 1). To ensure the transmission of force and moments from the CoM to the foot, the lower end of the spring and the spherical center of the hemispherical surface are rigidly tied. The upper end of the spring is hinged to the CoM, allowing the spring to rotate around the vertical axis freely in the three-dimensional space, with the vertical axis passing through the center of mass.Figure 1
Schematic of the three-dimensional walking model (double-support phase).According to the Lagrange formulation (Appendix A) and space velocity relation (Appendix B), the motion equation of the 3D model for the single-support phase can be obtained as follows:(1)mRsinθ1m0ml12+mR2+2l1mRcosθ1mRsinθ1000msinθ12l12θ..1l..1θ..11=Y1Y2Y3,(2)Y1=−mgcosθ1+kl0−l1+ml1θ˙12+sinθ12θ˙112−cl˙1,(3)Y2=mgsinθ1l1−m2Rcosθ1+l1l˙1θ˙1−sinθ1l1Rθ˙12+cosθ1l1θ˙112,(4)Y3=−2mcosθ1sinθ1l12θ˙1θ˙11−2msinθ12l1θ˙11l˙1.And the motion equation of the 3D model for the double-support phase can be obtained as follows:(5)mRsinθ1m0ml12+mR2+2l1mRcosθ1mRsinθ1000msinθ12l12θ..1l..1θ..11=Y1Y2Y3,(6)Y1=−mgcosθ1+kl0−l1+ml1θ˙12+sinθ12θ˙112−cl˙1+4kl0−l2−cl˙2Rcosθ1+cosθ1cosθ2+cosθ11−θ22sinθ1sinθ2l24Rcosθ2+2+2cos2θ2−cos2θ2−θ22+2cos2θ22−cos2θ2+θ22l2,(7)Y2=mgsinθ1l1−m2Rcosθ1+l1l˙1θ˙1−sinθ1l1Rθ˙12+cosθ1l1θ˙112+4kl0−l2−cl˙2−Rsinθ1l1+Rcosθ11−θ22sinθ2l2−cosθ2sinθ1l1l2+cosθ1cosθ11−θ22sinθ2l1l24Rcosθ2+2+2cos2θ2−cos2θ2−θ22+2cos2θ22−cos2θ2+θ22l2,(8)Y3=−2mcosθ1sinθ1l12θ˙1θ˙11−2msinθ12l1θ˙11l˙1+4kl0−l2−cl˙2sinθ1sinθ2sinθ11−θ22l1l24Rcosθ2+2+2cos2θ2−cos2θ2−θ22+2cos2θ22−cos2θ2+θ22l2,where m,R,θi,li,θii,θ˙i,l˙i,θ˙ii,l0,k, and c represent the model mass, the radius of the foot hemisphere, 2D spring angle, spring length, the 3D angles, the 2D angular velocity, the change rate of spring length, the 3D angular velocity, slack length of spring, spring stiffness, and damping constant, respectively. The subscripts i=1,2 define the trailing and leading leg (Figure 1). It is worth noting that the three mutually independent variables θ1,l1, and θ11 can determine the double-support phase of the model due to the existence of geometric and the dynamic constraint.


## 2.2. Simulation Method

As shown in Figure1, the positive direction of the x-axis is set as the initial direction of walking. With a similar pattern to that of the 2D model in [11, 22], the gait cycle in the 3D model is initiated from a single-support phase in which the 6 initial state parameters need to be determined. The parameters are θ1,θ˙1,θ11,θ˙11,l1, and l˙1, the details of the parameters have been described before. During the simulation for the integration of the equations, the 3D model is initially positioned in a state when the point mass is at its apex during a gait cycle. The forward limb is positioned in front of the trailing limb at a predetermined angle θtd. Before the contact of the leading leg with the ground, the motion equations (equation (1)) of the single-support phase are used to determine the dynamical behaviors of the model. It is worth noting that the angle θ1 of the 3D model is not initially zero, but a small value approaches zero, which is different from that in the 2D model to suit the need for convergence of the differential equation. In the numerical simulation, it is found that a stable result could be obtained when θ1 is in the range of 0.02 rad to 0.04 rad. Meanwhile, the initial angle θ11 and the corresponding angular velocity θ˙11 are selected through the Genetic Algorithm (GA) which will be described below.When the model is transformed from the single-support phase to the double-support phase, the following geometric constraint needs to be satisfied:(9)l1cosθ1=l2cosθ2.Especially when the leading leg is just in contact with the ground, the following condition needs to be satisfied:(10)l1cosθ1=l0cosθtd,where θtd represents the touch-down angle. The calculation of the dynamic states of the model in the double-support phase requires the closed kinematic constraints (equation (B.4) in Appendix B) to determine the parameters related to the limb spring velocity l˙2, 2D angular velocity θ.2, and 3D angular velocity θ˙22. If the trailing leg is lifted, and the point mass is at its apex, the trailing leg will be reset in front of the body with the predetermined impact angle. The model is instantaneously rotated along the plumb line passing through the body center, and the 3D angle and angular velocity after rotation are identical in value to those before the rotation but in the opposite direction. The process corresponds to the phase of “instant switch of legs” in Figure 2.Figure 2
Perspective (a), top view (b), and side view (c) of the walking model in the complete gait cycle.
(a)(b)(c)During a gait cycle, it is assumed that the change in the CoM and gait is symmetrical, which is consistent with the symmetric features of human walking. Based on the symmetric features and the consideration of computational efficiency, half of the gait cycle of the 3D model is calculated which is expressed as the single-double-single process in Figure2. By iterating the initial conditions to be consistent with the end conditions, a stable GRFs curve can be finally obtained. The vertical GRF deserves more attention than the transverse and longitudinal components of the GRF for the vertical GRF is assumed to better reflect the oscillation characteristics of human walking. To optimize the initial parameter values, the Genetic Algorithm (GA) Solver of MATLAB® is utilized to identify the optimal solution for a specified nonlinear system with the fitness function set as the root mean square error between the data of the simulation of vertical GRF and the experimental ones. To prevent blind selection of parameters from causing nonconvergence of the solution, the parameter range is obtained from trial calculations by reference to the results of [22], and the final parameter range is as shown in Table 1. The parameters of GA are given as follows: the population size is set to 100, the mutation rate to 0.01, and the maximum number of generations for each run to 50.Table 1
Normalized model parameters for genetic algorithm calculations.
Global model parameters:M=80kg,l0=0.8m, andR=0.3mLeg stiffness,kN/m17500∼22500Damping ratio,ξ0.02∼0.05Initial leg angle,θ1rad0.02∼0.04Initial leg angular velocity,θ˙1rad/s1.25∼2.00Initial3D leg angle,θ11rad0.3∼0.6Initial3D leg angular velocity,θ˙11rad/s−1.20∼−0.40Touch−down angle,θ2degree18∼24


## 2.3. Experiment

To validate the model proposed for 3D walking, the experimental data were collected to compare with the simulation results. Thirteen young and healthy subjects (body mass of 67.4 ± 12.2 kg and height of 1.67 ± 0.12 m) with no gait disorders participated in the experiment. The subjects walked on a walkway equipped with three force plates with piezoelectric 3-component force sensors. The participants were instructed to keep pace with a metronome placed nearby, which provided a fixed frequency each time. During the experiment, each participant was guided to walk at the frequency of the metronome, three times at each frequency. To get the ground reaction force (GRF) of each participant, three six-dimensional force plates (JP6060, CN) were used to measure the real-time GRF at a sample rate of 100 Hz. To ensure that the participants would cross the force plate at a comfortable and stable pace, the pre- and post-transition parts of the walkway made of wood were placed before and after the force plate (Figure3). An AHRS system (AH100B, CN) was used to measure the acceleration and direction angle of CoM during the walking process at a sample rate of 100 Hz, as shown in Figure 4. The AHRS system could accurately and stably perform real-time 3D dynamic measurements with a Kalman filter. The velocity of the CoM could be calculated by integrating the acceleration data collected by the AHRS system.Figure 3
Test site layout of the walkway.Figure 4
Participant walks on the walkway at a fixed frequency.


## 3. Results

### 3.1. The Experimental Results Fit Well with the Numerical Simulation of the 3D Walking Model

To minimize the differentiation of the GRF due to the variation of participants in mass, the GRF is standardized by dividing the ground reaction force by the bodyweight as an evaluation criterion (Figure5). With the initial parameters optimized by the genetic algorithm (GA) method, the simulation results are consistent with the experimental results well, with the root mean square error of 0.065±0.017 for vertical forces. As can be seen from Figures 5(a) and 5(c), for GRF in the vertical and longitudinal direction, the simulation fits the data well in peak and valley values. There is a minor difference between the observed and simulated data in the valley values of the lateral GRF, which can be identified from Figure 5(b).Figure 5
Comparison of the experimental and simulated data on GRF. (a) Vertical forces. (b) Lateral forces. (c) Forward forces. The average and standard deviation of the experimental data are indicated by thin solid lines and gray shade, respectively, and the average of the simulated data is indicated by thick dashed lines.
(a)(b)(c)


### 3.2. Effect of Model Dynamics Parameters on GRF Shape

The effects of spring stiffness and gait frequency on the GRF simulation results are shown in Figure6(a). As demonstrated in Figure 6(a), the greater stiffness and gait frequency bring higher peaks, shorter gait cycles, and a more distinct gap between the peak and valley values but have little impact on the lateral and longitudinal forces. However, there is a significant effect of the 3D dynamics parameters on lateral forces, which is indicated in Figure 6(b). The results in Figure 6(b) reveal that the lateral GRF exhibits an increase in peaks with the increase of 3D angles and a decrease of 3D angular velocity. The effect of other parameters on lateral forces will be discussed later in Figure 7.Figure 6
Simulation of GRF for different oscillation frequencies and initial parameters. (a) GRF in three directions. (b) GRF in the lateral direction. The leg stiffness ranges from 18 kN/m to 20 kN/m and then 22 kN/m, respectively. The vertical, forward, and lateral forces are indicated by thick, thin solid, and dashed lines, respectively (a). By setting the 3D angular velocity to −0.1 rad/s, −0.8 rad/s, and −1.5 rad/s and the 3D angles to 0.022 rad, 0.164 rad, and 0.317 rad respectively, the simulation result of lateral forces can be obtained in (b).
(a)(b)Figure 7
Effect of initial 3D angular velocity and initial 2D angle on the normalized peak values of GRF in the lateral direction: (a) initial 3D angleθ11=0.4rad; (b) initial 3D angular velocity θ˙11=−0.9rad/s.
(a)(b)


### 3.3. Correlations between the Initial Dynamic Parameters of the 3D Model

To explore the potential relationship between various parameters of the model and features of human gait, the kinetic parameterfexp is defined as 1/Tc1, and Tc1 is obtained by experiment (Figure 8). A significantly positive correlation (r = 0.793, p<0.0001, r = correlation coefficient) between walking velocity vexp and gait frequency fexp is observed (Figure 9(a)).Figure 8
Schematic illustration of the definition of the gait frequency.Figure 9
(A) Relationship between walking velocity and gait frequency. (B) Relationship between spring stiffness constants and gait frequency. (C) Relationship between initial angular velocity and gait frequency. (D) Relationship between impact angle and gait frequency.
(a)(b)(c)(d)When the kinetic parameters are identified by the GA to achieve a better fit between the simulations and the experimental results, a significant positive linear correlation can be observed between gait frequency and spring stiffness (r = 0.597, p<0.0001) (Figure 9(b)) and initial 2D angular velocity (r = 0.338, p<0.05) (Figure 9(c)). Meanwhile, there is a significant negative correlation between the impact angle and gait frequency (r = −0.469, p<0.0001), as illustrated in Figure 9(d).Figure10(a) shows the impact angle versus the initial 2D angle, with a correlation coefficient of 0.394 (p<0.05). In contrast, a negative correlation is observed between impact angle and initial 3D angle (r = −0.306, p<0.05) (Figure 10(b)) and impact angle and damping ratio (r = −0.372, p<0.01) (Figure 10(c)).Figure 10
(A) Relationship between initial 2D angle and impact angle. (B) Relationship between initial 3D angle and impact angle. (C) Relationship between the damping ratio and the impact angle.
(a)(b)(c)Figures9 and 10 show the correlation coefficient of the linear fit results for some of the dynamic parameters. As shown in Figure 11, a scatter plot matrix is designed to detect possible correlations among the initial parameters of the 3D model optimized by the genetic algorithm.Figure 11
Scatterplot matrix: Correlation between dynamic and gait characteristics of the 3D mode.A scatter plot matrix is an effective multivariate visualization tool that plots the correlation between all variables. The nondiagonal part of theith row and jth column represents the scatter plot of the ith variable against the jth variable, which is intended to describe the correlation between these two variables. If the two variables are positively correlated, the dependent variable y increases as the independent variable x becomes larger. If the two variables are negatively correlated, the dependent variable y becomes smaller as the independent variable x becomes larger. If the two variables are uncorrelated, the dependent variable y does not change as the independent variable x changes. In Figure 11, the solid line represents a least square fit of the data to a straight line, and the positive or negative sign of its slope determines whether the two variables are positively or negatively correlated. The thick dashed line indicates a smooth nonparametric curve under a smoother loess (the loess span 0.75) and the fine dashed line indicates the boundaries of the 95% confidence interval of the smoothed curve. The box plots of the dynamic parameters of the 3D model are distributed diagonally. From Figure 11, it can be drawn that all parameters are correlated to some extent.To illustrate the level of correlation between the optimized initial variables, the linear fit results of the 10 parameters are used to generate a heat map (Figure12). Correlation analysis in the heat map (Figure 12) shows that walking velocity is positively correlated with gait frequency (r = 0.79, p<0.001), step length (r = 0.63, p<0.001), and leg stiffness (r = 0.33, p<0.5), but has a negative correlation with impact angle (r = −0.40, p<0.01). Among the initial parameters, only the impact angle and initial 3D angular velocity are correlated with the damping ratio (r = −0.37 and r = 0.39, p<0.05). The impact angle is correlated with initial variables, except 2D and 3D angular velocity and step length. It is worth noting that in addition to its positive correlation with velocity, gait frequency is also positively correlated with leg stiffness (r = 0.6, p<0.001) and initial 2D angular velocity (r = 0.34, p<0.05), and there is a relatively significant negative correlation with impact angle (r = −0.47, p<0.001). There is a weak positive correlation between initial 2D angle and initial 3D angle (r = 0.32, p<0.05) and a weak negative correlation between the initial 2D angular velocity and initial 3D angular velocity (r = −0.38, p<0.01).Figure 12
Heat map showing the correlation among the dynamic and gait characteristics of the 3D walking model. The heat map (a) with the level of significance indicated by asterisks (p∗<0.05, p∗∗<0.01, and p∗∗∗<0.001) and insignificant coefficients (p>0.05) marked with the corresponding p value. The heat map (b) with the correlation coefficient indicated by the color code and marked with the r value.
(a)(b)


### 3.4. Effect of Initial Dynamic Parameters on GRFs

#### 3.4.1. Effect of Initial Dynamic Parameters on GRF in the Longitudinal Direction

Kinetics studies are conducted to investigate the effect of various initial conditions on GRF in the longitudinal, vertical, and lateral direction by setting the initial walking velocity asv0=1.5m/s and initial 2D angular velocity as θ˙1=v0/L0+R for estimation.The impacts that initial 3D angular velocity has on GRF in the longitudinal direction are grouped here into two phases (Figure13). Phase I corresponds to the case where initial 3D angular velocity and peak forward forces show a negative correlation. In Phase II, with the increasing 3D angular velocity, the GRF values first increase till reaching a peak value and then decrease, which is exhibited in Figure 13(a) (θ1=0.06), Figure 13(b) (θ1>0.045), and Figure 13(c) (θ1>0.04). It is manifested that the peak values tend to increase with increasing initial 2D angle θ1. The effects of different initial factors on GRF peaks in the longitudinal direction are explored under the following conditions: damping ratio, 0.04; impact angle, 22 deg; leg stiffness, 20000 N/m; human mass, 80 kg; spring resting length, 0.8 m; and radius of foot, 0.3 m.Figure 13
Effect of initial 3D angular velocity on the normalized peak values of GRF in the longitudinal direction: (a) initial 3D angleθ11=0.2rad, (b) initial 3D angle θ11=0.4rad, and (c) initial 3D angle θ11=0.6rad.
(a)(b)(c)


#### 3.4.2. Effect of Initial Dynamic Parameters on GRF in the Vertical Direction

As Figure14 shows, the initial 3D angular velocity has a similar impact on vertical forces. An increase in the 3D angular velocity results in the decrease of vertical forces when the initial angle θ1 is small in Figure 14(a) (θ1<0.04), Figure 14(b) (θ1<0.05), and Figure 14(c) (θ1<0.06). In contrast, the initial 3D angular velocity has little effect on the peak values when the initial angle θ1 is large, as shown in Figure 14(b) (θ1≥0.05) and Figure 14(c) (θ1≥0.06). However, the trend is different when θ11=0.2, and the vertical values first increase and then decrease as θ˙11 grows when θ1=0.06 and remain basically unchanged when 0.04≤θ1<0.06 (Figure 14(a)). As evident from Figure 14, with a fixed initial 3D angle θ11, the larger the θ1, the less influence the initial 3D angular velocity will have on the vertical peak values. The effects of different initial factors on GRF peaks in the vertical direction are explored under the following conditions: damping ratio, 0.04; impact angle, 22 deg; leg stiffness, 20000 N/m; human mass, 80 kg; spring resting length, 0.8 m; and radius of foot, 0.3 m.Figure 14
Effect of initial 3D angular velocity on the normalized peak values of GRF in the vertical direction: (a) initial 3D angleθ11=0.2rad, (b) initial 3D angle θ11=0.4rad, and (c) initial 3D angle θ11=0.6rad.
(a)(b)(c)


#### 3.4.3. Effect of Initial Dynamic Parameters on GRF in the Lateral Direction

Figure7(a) reveals that an increase in the initial 3D angular velocity results in a decrease of peak values of GRF in the lateral direction while the decreasing initial θ1 shows an opposite effect. The initial 3D angle θ11 is expected to have a negligible impact on the GRF values, which is illustrated in Figure 7(b). The effects of different initial factors on GRF peaks in lateral direction are explored under the following conditions: damping ratio, 0.04; impact angle, 22 deg; leg stiffness, 20000 N/m; human mass, 80 kg; spring resting length, 0.8 m; and radius of foot, 0.3 m.


## 3.1. The Experimental Results Fit Well with the Numerical Simulation of the 3D Walking Model

To minimize the differentiation of the GRF due to the variation of participants in mass, the GRF is standardized by dividing the ground reaction force by the bodyweight as an evaluation criterion (Figure5). With the initial parameters optimized by the genetic algorithm (GA) method, the simulation results are consistent with the experimental results well, with the root mean square error of 0.065±0.017 for vertical forces. As can be seen from Figures 5(a) and 5(c), for GRF in the vertical and longitudinal direction, the simulation fits the data well in peak and valley values. There is a minor difference between the observed and simulated data in the valley values of the lateral GRF, which can be identified from Figure 5(b).Figure 5
Comparison of the experimental and simulated data on GRF. (a) Vertical forces. (b) Lateral forces. (c) Forward forces. The average and standard deviation of the experimental data are indicated by thin solid lines and gray shade, respectively, and the average of the simulated data is indicated by thick dashed lines.
(a)(b)(c)


## 3.2. Effect of Model Dynamics Parameters on GRF Shape

The effects of spring stiffness and gait frequency on the GRF simulation results are shown in Figure6(a). As demonstrated in Figure 6(a), the greater stiffness and gait frequency bring higher peaks, shorter gait cycles, and a more distinct gap between the peak and valley values but have little impact on the lateral and longitudinal forces. However, there is a significant effect of the 3D dynamics parameters on lateral forces, which is indicated in Figure 6(b). The results in Figure 6(b) reveal that the lateral GRF exhibits an increase in peaks with the increase of 3D angles and a decrease of 3D angular velocity. The effect of other parameters on lateral forces will be discussed later in Figure 7.Figure 6
Simulation of GRF for different oscillation frequencies and initial parameters. (a) GRF in three directions. (b) GRF in the lateral direction. The leg stiffness ranges from 18 kN/m to 20 kN/m and then 22 kN/m, respectively. The vertical, forward, and lateral forces are indicated by thick, thin solid, and dashed lines, respectively (a). By setting the 3D angular velocity to −0.1 rad/s, −0.8 rad/s, and −1.5 rad/s and the 3D angles to 0.022 rad, 0.164 rad, and 0.317 rad respectively, the simulation result of lateral forces can be obtained in (b).
(a)(b)Figure 7
Effect of initial 3D angular velocity and initial 2D angle on the normalized peak values of GRF in the lateral direction: (a) initial 3D angleθ11=0.4rad; (b) initial 3D angular velocity θ˙11=−0.9rad/s.
(a)(b)


## 3.3. Correlations between the Initial Dynamic Parameters of the 3D Model

To explore the potential relationship between various parameters of the model and features of human gait, the kinetic parameterfexp is defined as 1/Tc1, and Tc1 is obtained by experiment (Figure 8). A significantly positive correlation (r = 0.793, p<0.0001, r = correlation coefficient) between walking velocity vexp and gait frequency fexp is observed (Figure 9(a)).Figure 8
Schematic illustration of the definition of the gait frequency.Figure 9
(A) Relationship between walking velocity and gait frequency. (B) Relationship between spring stiffness constants and gait frequency. (C) Relationship between initial angular velocity and gait frequency. (D) Relationship between impact angle and gait frequency.
(a)(b)(c)(d)When the kinetic parameters are identified by the GA to achieve a better fit between the simulations and the experimental results, a significant positive linear correlation can be observed between gait frequency and spring stiffness (r = 0.597, p<0.0001) (Figure 9(b)) and initial 2D angular velocity (r = 0.338, p<0.05) (Figure 9(c)). Meanwhile, there is a significant negative correlation between the impact angle and gait frequency (r = −0.469, p<0.0001), as illustrated in Figure 9(d).Figure10(a) shows the impact angle versus the initial 2D angle, with a correlation coefficient of 0.394 (p<0.05). In contrast, a negative correlation is observed between impact angle and initial 3D angle (r = −0.306, p<0.05) (Figure 10(b)) and impact angle and damping ratio (r = −0.372, p<0.01) (Figure 10(c)).Figure 10
(A) Relationship between initial 2D angle and impact angle. (B) Relationship between initial 3D angle and impact angle. (C) Relationship between the damping ratio and the impact angle.
(a)(b)(c)Figures9 and 10 show the correlation coefficient of the linear fit results for some of the dynamic parameters. As shown in Figure 11, a scatter plot matrix is designed to detect possible correlations among the initial parameters of the 3D model optimized by the genetic algorithm.Figure 11
Scatterplot matrix: Correlation between dynamic and gait characteristics of the 3D mode.A scatter plot matrix is an effective multivariate visualization tool that plots the correlation between all variables. The nondiagonal part of theith row and jth column represents the scatter plot of the ith variable against the jth variable, which is intended to describe the correlation between these two variables. If the two variables are positively correlated, the dependent variable y increases as the independent variable x becomes larger. If the two variables are negatively correlated, the dependent variable y becomes smaller as the independent variable x becomes larger. If the two variables are uncorrelated, the dependent variable y does not change as the independent variable x changes. In Figure 11, the solid line represents a least square fit of the data to a straight line, and the positive or negative sign of its slope determines whether the two variables are positively or negatively correlated. The thick dashed line indicates a smooth nonparametric curve under a smoother loess (the loess span 0.75) and the fine dashed line indicates the boundaries of the 95% confidence interval of the smoothed curve. The box plots of the dynamic parameters of the 3D model are distributed diagonally. From Figure 11, it can be drawn that all parameters are correlated to some extent.To illustrate the level of correlation between the optimized initial variables, the linear fit results of the 10 parameters are used to generate a heat map (Figure12). Correlation analysis in the heat map (Figure 12) shows that walking velocity is positively correlated with gait frequency (r = 0.79, p<0.001), step length (r = 0.63, p<0.001), and leg stiffness (r = 0.33, p<0.5), but has a negative correlation with impact angle (r = −0.40, p<0.01). Among the initial parameters, only the impact angle and initial 3D angular velocity are correlated with the damping ratio (r = −0.37 and r = 0.39, p<0.05). The impact angle is correlated with initial variables, except 2D and 3D angular velocity and step length. It is worth noting that in addition to its positive correlation with velocity, gait frequency is also positively correlated with leg stiffness (r = 0.6, p<0.001) and initial 2D angular velocity (r = 0.34, p<0.05), and there is a relatively significant negative correlation with impact angle (r = −0.47, p<0.001). There is a weak positive correlation between initial 2D angle and initial 3D angle (r = 0.32, p<0.05) and a weak negative correlation between the initial 2D angular velocity and initial 3D angular velocity (r = −0.38, p<0.01).Figure 12
Heat map showing the correlation among the dynamic and gait characteristics of the 3D walking model. The heat map (a) with the level of significance indicated by asterisks (p∗<0.05, p∗∗<0.01, and p∗∗∗<0.001) and insignificant coefficients (p>0.05) marked with the corresponding p value. The heat map (b) with the correlation coefficient indicated by the color code and marked with the r value.
(a)(b)


## 3.4. Effect of Initial Dynamic Parameters on GRFs

### 3.4.1. Effect of Initial Dynamic Parameters on GRF in the Longitudinal Direction

Kinetics studies are conducted to investigate the effect of various initial conditions on GRF in the longitudinal, vertical, and lateral direction by setting the initial walking velocity asv0=1.5m/s and initial 2D angular velocity as θ˙1=v0/L0+R for estimation.The impacts that initial 3D angular velocity has on GRF in the longitudinal direction are grouped here into two phases (Figure13). Phase I corresponds to the case where initial 3D angular velocity and peak forward forces show a negative correlation. In Phase II, with the increasing 3D angular velocity, the GRF values first increase till reaching a peak value and then decrease, which is exhibited in Figure 13(a) (θ1=0.06), Figure 13(b) (θ1>0.045), and Figure 13(c) (θ1>0.04). It is manifested that the peak values tend to increase with increasing initial 2D angle θ1. The effects of different initial factors on GRF peaks in the longitudinal direction are explored under the following conditions: damping ratio, 0.04; impact angle, 22 deg; leg stiffness, 20000 N/m; human mass, 80 kg; spring resting length, 0.8 m; and radius of foot, 0.3 m.Figure 13
Effect of initial 3D angular velocity on the normalized peak values of GRF in the longitudinal direction: (a) initial 3D angleθ11=0.2rad, (b) initial 3D angle θ11=0.4rad, and (c) initial 3D angle θ11=0.6rad.
(a)(b)(c)


### 3.4.2. Effect of Initial Dynamic Parameters on GRF in the Vertical Direction

As Figure14 shows, the initial 3D angular velocity has a similar impact on vertical forces. An increase in the 3D angular velocity results in the decrease of vertical forces when the initial angle θ1 is small in Figure 14(a) (θ1<0.04), Figure 14(b) (θ1<0.05), and Figure 14(c) (θ1<0.06). In contrast, the initial 3D angular velocity has little effect on the peak values when the initial angle θ1 is large, as shown in Figure 14(b) (θ1≥0.05) and Figure 14(c) (θ1≥0.06). However, the trend is different when θ11=0.2, and the vertical values first increase and then decrease as θ˙11 grows when θ1=0.06 and remain basically unchanged when 0.04≤θ1<0.06 (Figure 14(a)). As evident from Figure 14, with a fixed initial 3D angle θ11, the larger the θ1, the less influence the initial 3D angular velocity will have on the vertical peak values. The effects of different initial factors on GRF peaks in the vertical direction are explored under the following conditions: damping ratio, 0.04; impact angle, 22 deg; leg stiffness, 20000 N/m; human mass, 80 kg; spring resting length, 0.8 m; and radius of foot, 0.3 m.Figure 14
Effect of initial 3D angular velocity on the normalized peak values of GRF in the vertical direction: (a) initial 3D angleθ11=0.2rad, (b) initial 3D angle θ11=0.4rad, and (c) initial 3D angle θ11=0.6rad.
(a)(b)(c)


### 3.4.3. Effect of Initial Dynamic Parameters on GRF in the Lateral Direction

Figure7(a) reveals that an increase in the initial 3D angular velocity results in a decrease of peak values of GRF in the lateral direction while the decreasing initial θ1 shows an opposite effect. The initial 3D angle θ11 is expected to have a negligible impact on the GRF values, which is illustrated in Figure 7(b). The effects of different initial factors on GRF peaks in lateral direction are explored under the following conditions: damping ratio, 0.04; impact angle, 22 deg; leg stiffness, 20000 N/m; human mass, 80 kg; spring resting length, 0.8 m; and radius of foot, 0.3 m.


## 3.4.1. Effect of Initial Dynamic Parameters on GRF in the Longitudinal Direction

Kinetics studies are conducted to investigate the effect of various initial conditions on GRF in the longitudinal, vertical, and lateral direction by setting the initial walking velocity asv0=1.5m/s and initial 2D angular velocity as θ˙1=v0/L0+R for estimation.The impacts that initial 3D angular velocity has on GRF in the longitudinal direction are grouped here into two phases (Figure13). Phase I corresponds to the case where initial 3D angular velocity and peak forward forces show a negative correlation. In Phase II, with the increasing 3D angular velocity, the GRF values first increase till reaching a peak value and then decrease, which is exhibited in Figure 13(a) (θ1=0.06), Figure 13(b) (θ1>0.045), and Figure 13(c) (θ1>0.04). It is manifested that the peak values tend to increase with increasing initial 2D angle θ1. The effects of different initial factors on GRF peaks in the longitudinal direction are explored under the following conditions: damping ratio, 0.04; impact angle, 22 deg; leg stiffness, 20000 N/m; human mass, 80 kg; spring resting length, 0.8 m; and radius of foot, 0.3 m.Figure 13
Effect of initial 3D angular velocity on the normalized peak values of GRF in the longitudinal direction: (a) initial 3D angleθ11=0.2rad, (b) initial 3D angle θ11=0.4rad, and (c) initial 3D angle θ11=0.6rad.
(a)(b)(c)


## 3.4.2. Effect of Initial Dynamic Parameters on GRF in the Vertical Direction

As Figure14 shows, the initial 3D angular velocity has a similar impact on vertical forces. An increase in the 3D angular velocity results in the decrease of vertical forces when the initial angle θ1 is small in Figure 14(a) (θ1<0.04), Figure 14(b) (θ1<0.05), and Figure 14(c) (θ1<0.06). In contrast, the initial 3D angular velocity has little effect on the peak values when the initial angle θ1 is large, as shown in Figure 14(b) (θ1≥0.05) and Figure 14(c) (θ1≥0.06). However, the trend is different when θ11=0.2, and the vertical values first increase and then decrease as θ˙11 grows when θ1=0.06 and remain basically unchanged when 0.04≤θ1<0.06 (Figure 14(a)). As evident from Figure 14, with a fixed initial 3D angle θ11, the larger the θ1, the less influence the initial 3D angular velocity will have on the vertical peak values. The effects of different initial factors on GRF peaks in the vertical direction are explored under the following conditions: damping ratio, 0.04; impact angle, 22 deg; leg stiffness, 20000 N/m; human mass, 80 kg; spring resting length, 0.8 m; and radius of foot, 0.3 m.Figure 14
Effect of initial 3D angular velocity on the normalized peak values of GRF in the vertical direction: (a) initial 3D angleθ11=0.2rad, (b) initial 3D angle θ11=0.4rad, and (c) initial 3D angle θ11=0.6rad.
(a)(b)(c)


## 3.4.3. Effect of Initial Dynamic Parameters on GRF in the Lateral Direction

Figure7(a) reveals that an increase in the initial 3D angular velocity results in a decrease of peak values of GRF in the lateral direction while the decreasing initial θ1 shows an opposite effect. The initial 3D angle θ11 is expected to have a negligible impact on the GRF values, which is illustrated in Figure 7(b). The effects of different initial factors on GRF peaks in lateral direction are explored under the following conditions: damping ratio, 0.04; impact angle, 22 deg; leg stiffness, 20000 N/m; human mass, 80 kg; spring resting length, 0.8 m; and radius of foot, 0.3 m.


## 4. Discussion

The application of the conventional SLIP model is restricted to a two-dimensional plane, and the model can only predict the trajectory of CoM and the corresponding values of GRF in the vertical and longitudinal directions [11, 22] in the human walking process. In the conventional model, a curved foot is employed to mimic the pressure center excursions [8, 11, 22, 42]. In this article, we introduce the spatial angle and the corresponding angular velocity to construct a 3D model, with the foot treated as a hemispheric surface (Figure 1). Similar to the traditional SLIP model, the proposed model simulates bipedal walking through a simple mass, spring, and damper system, with legs equivalent to springs and dampers. The difference between the original and the proposed model lies in the fact that the proposed one can simulate the gait process and take the 3D dynamic parameters into account.With the initial parameters optimized by the genetic algorithm (GA) method, the simulation results are consistent with the experimental results of the GRF well, which validates the 3D model. Due to the introduction of 3D kinematic parameters, the model proposed in this paper can generate GRF forces in three directions, rather than a calculation of the forces just in two directions conducted in previous models [8–10, 22, 24, 32, 43–45]. The initial dynamics parameters of the 3D model have a significant effect on the shape and peak of the GRF, which is in line with those results presented in previous studies [22, 34, 46].To thoroughly explore the relationship between several dynamic parameters, the correlation analysis is carried out with a linear correlation (Pearson’sr) test. The results show that there is a correlation between some of the initial parameters (Figure 12). A positive significant correlation is found between walking velocity and gait frequency and walking velocity and step length, a pattern similar to previous studies [33, 46, 47]. The results reveal a positive association between the frequency and stiffness and initial 2D angular velocity and a negative association with the impact angle. There also exists a negative and significant correlation between impact angle and frequency, velocity, stiffness, and damping ratio, in line with [22]. The initial 2D angle is significantly positively correlated with the impact angle and initial 3D angle; 2D angular velocity is negatively correlated with 3D angular velocity. A moderate positive correlation is found between damping and 3D angular velocity; a moderate negative correlation is found between impact angle and initial 3D angle. Regarding leg stiffness k, our research suggests that there is a positive correlation between walking speed and stiffness, which is consistent with the results of [33, 48]. However, the finding that the damping ratio is not significantly associated with walking speed does not support the previous research [33]. The possible reason is that the spring constant k taken by [33] is in the range of 15–30 kN/m, but we take it in the range of 17.5–22.5 kN/m. Contrary to the previous study [11], our study shows that there is a correlation between the initial angle, initial angular velocity, and damping with walking speed. The existence of 3D angles and 3D angular velocities may explain such an inconsistency.A detailed parametric study is conducted to analyze the impact of initial dynamic parameter values on the GRF peaks. The results show that both 2D and 3D parameters have a nonnegligible effect on the peak of the GRF, but relative to the GRF in the vertical and longitudinal direction, the 3D parameters have the most significant effect on the GRF in the lateral direction. The possible explanation for the sensitivity of the lateral GRF to the 3D angular velocity may be as follows. The increase in absolute values of 3D angular velocity means that under the same impact angle and time for leg switching, the 3D angle will increase when the foot against the ground. Since the 3D angle is small and the lateral force is the product of the sine of the 3D angle and the spring force in the horizontal projection, the lateral force shows a positive correlation with the 3D angle and a negative correlation with 3D angular velocity (note: the 3D angular velocity is a negative value).It must be acknowledged that our study has some limitations. First, the 3D model results in computational complexity for its adding of 3D variables entails the restriction of variables range. Second, the 3D model cannot perform a complete gait cycle continuously and spontaneously if compared to the model proposed by Lim and Park [48]. Third, the objective function of the genetic algorithm is set to be the root mean square error between the simulation results and the experimental data of the vertical GRF. Despite the advantages of simplified calculations, it must be conceded that the optimization approach may omit some correlation features between the parameters.


## 5. Conclusions

In this study, we proposed a three-dimensional bipedal walking model that could generate the repeated gait cycle with the simulated results of GRF in three directions fitting with the experimental data well. The correlation analysis between the initial dynamic parameters of the model was conducted. We examined the impact of 3D dynamics parameters on the peaks of GRF in three directions and found that the 3D parameters had a major effect on the lateral GRF. The results of our study implied that the dynamic characteristics presented by the bipedal walking pedestrian could be simulated by a spring-damping system.



---

*Source: 6651715-2021-07-15.xml*