Abstract

In this work, the lower extremity physiological parameters are recorded during normal walking gait, and the dynamical systems theory is applied to determine a stability analysis. The human walking gait pattern of kinematic and dynamical data is approximated to periodic behavior. The embedding dimension analysis of the kinematic variable's time trace and use of Taken's theorem allows us to compute a reduced-order time series that retains the essential dynamics. In conjunction with Floquet theory, this approach can help determine the system's stability characteristics. The Lyapunov–Floquet (L-F) transformation application results in constructing an invariant manifold resembling the form of a simple oscillator system. It is also demonstrated that the simple oscillator system, when re-mapped back to the original domain, reproduces the original system's time evolution (hip angle or knee angle, for example). A reinitialization procedure is suggested that improves the accuracy between the processed data and actual data. The theoretical framework proposed in this work is validated with the experiments using a motion capture system.

1 Introduction

The study of human walking gait is an active area of research. For the applications of designing exoskeletons, rehabilitation devices, bionics, prostheses, and robots that imitate human walking, there is a need to understand the inherent dynamical characteristics and behavior of human walking. The dynamical analysis of gait patterns enables us to distinguish stable walking from unstable walking. It also aids in building rehabilitation devices for patients with an unsteady gait, with ailments, such as stroke, spinal cord injury (SCI) conditions, and Parkinson’s disease [1]. These dynamical stability measures could also be applied in designing efficient humanoid robots and mechanisms that could imitate human gait flawlessly. Hence, there is a need to understand human walking gait behavior to facilitate impeccable human–robot interaction.

The earliest dynamical analyses on human gait patterns were performed on the biped models based on the kinematic data [2,3]. The kinematic behavior of a biped was considered as a periodic time-varying system, and Floquet theory was used. A Poincare section was taken, from the phase plots of joint angle variations, at various gait events (such as heel strike and toe-off), and stability of the gait was determined based on the Floquet Multipliers. The divergence of neighboring trajectories in the phase plot was measured using the Lyapunov exponents [4]. This technique was also utilized in determining the local stability of human kinematic gait data variation, along with the Floquet Multipliers, in Refs. [5,6]. Researchers also used the variation of maximum Lyapunov exponents from the kinematic data and correlated it to fall detection in the biped simulation models [7]. The dynamical systems theory’s application for evaluating pathological coordination in movement disorder was demonstrated in Ref. [8]. Researchers applied a robust control design, dependent on sagittal plane hip angles, using the Lyapunov approach for bounded nonlinearity in a phase-based hip exoskeleton control in Ref. [9]. Hence, various stability analysis techniques have been used to investigate human gait patterns based on kinematic data. However, in specific applications, both the kinematic and dynamic features are required for comprehensive stability characterization.

Researchers, in the past, have worked toward imitating the human gait patterns in prosthesis and humanoid robots, using multiple approaches, such as optimization methods, dynamic simulations, predictive models with the aid of dynamic parameters [10]. In Ref. [11], a ground reaction force (GRF) estimation technique was proposed to control a lower-limb prosthesis, and the stability characteristics were verified using the Lyapunov analysis. A humanoid robot was balanced based on the real-time normal and tangential GRF data obtained from wireless instrumented shoes [12]. In Ref. [13], the GRF data collected from a legged robot was utilized in conjunction with a support vector machine (SVM) algorithm in terrain classification. The GRF data were also considered as the parameter for evaluating asymmetry, in the unilateral and bilateral movements for anterior cruciate ligament (ACL) injured patients, during rehabilitation in Refs. [14,15]. The chances of re-occurrence of the same injury after ACL reconstruction was predicted by evaluating the largest Lyapunov exponent from GRF data in Ref. [16]. Thus, it can be inferred that the GRF plays a vital role in the stability and control of human gait and posture’s dynamical behavior. The GRF is also closely correlated with other dynamical features, such as joint torques and moments.

The kinematic features of the human walking gait are collected using laboratory equipment or wearable devices. These features are recorded using these devices as time series. Rosenstein et al., in their article [4], emphasized the use of time-delayed embedding of time series data for proper reconstruction of phase space for the original system variables. In Ref. [17], the application of time-delayed function was considered to assess the dynamical stability in assistive robots. Researchers have also used the time-delay embedding technique for activity recognition in human gait [18]. Researchers have also worked on generating a model-based control system for gait applications [19]. In Ref. [20], a 5-link biped walker was modeled as a periodic system, with an impulse effect that lay within a hybrid invariant manifold, which was later analyzed using dynamical systems theory. In Ref. [21], both the kinematic and dynamic parameters were utilized to design a humanoid biped robot.

An equivalent time-varying system can model the discrete-time series data to understand the human walking gait from a dynamical systems perspective. In this work, the human walking gait pattern is studied using the dynamical system theory, with the following objectives:

  • Consider the human walking behavior as a time-periodic system and perform a conformal mapping to convert it to a time-invariant dynamical system.

  • Evaluate the dynamical characteristics of the human gait pattern using the time-invariant system and validate the periodic behavior.

  • Back-transform/re-map the time-invariant dynamical system to the original coordinates to generate the system state's closed-form solution and validate with the time series data.

  • Demonstrate the proposed approach's application using the lower body joint angles as the kinematic features and ground reaction forces as the dynamic features on the data collected from healthy human subjects.

The application of dynamical system tools to study human gait would better explain the inherent biomechanical behavior. This dynamical systems theory approach would provide insights into efficient exoskeleton design and bionics for rehabilitation applications. Since this work is based on the time series of physiological parameters, a time-delayed embedding method is used to ensure that the conformal mapping onto an invariant manifold captures all the original system features. The dynamical characteristics are evaluated using Floquet theory. A brief description of various concepts, methodology, test results, and conclusive remarks for this work, are provided in the Secs. 24.

2 Mathematical Background

In this section, multiple concepts applicable to this work are briefly explained.

2.1 Dynamical System Theory.

The human gait system is a dynamically complicated system to analyze. Researchers have used different ways to approach the problem while approximating the system’s behavior [8,2231]. In Ref. [29], the human walking gait model was compared with a forced Van der Pol oscillator equation. Similarly, other researchers have expressed human walking behavior as a time-varying dynamical system. In this work, the human walking gait is expressed in the generic state-space form, considering a time-periodic linear part and the nonlinear and forcing terms separately, as indicated below:
x˙=A(t)x+E(x,t)+G(t),xRn
(1)
where A(t) is T periodic, i.e., A(t)=A(t+T). E(x,t) represents the nonlinear part of the system and G(t) represents the forcing term [31]. The nonlinearity in human walking behavior is assumed to correspond to the irregularity and variation in the walking surface or speed. In the case of varying gait speeds or irregular surfaces, the nonlinear part becomes prominent and needs to be accounted for in the human walking gait pattern analysis. The authors propose to use a Lyapunov-based approach [32] to quantify and verify the bounding conditions of the nonlinearity. In the case of walking on an instrumented treadmill, such effects can be assumed to be negligible. The external forcing component is assumed to correspond to the gait assistance provided by assistive devices or robotic systems, such as exoskeletons. The linear form of Eq. (1) without any forcing and nonlinear terms could be represented by
x˙(t)=A(t)x(t)
(2)

Further analysis performed in this article assumes human walking gait as a time-periodic linear system with no external force and no nonlinear terms. In this work, the human walking gait’s local stability on an instrumented treadmill is analyzed for the dynamical characteristics. Later the model is used to verify the transformation of the gait variables to a dynamically time-invariant form.

2.2 Lyapunov–Floquet Theory.

The stability and response of a linear time-varying system with periodic coefficients are evaluated using Floquet theory. Consider a linear time-varying system of the form indicated in (2). Let Φ(t) be the state transition matrix (STM), such that it satisfies (2) the initial condition Φ(0) = I. The solution (2) can be written as
x(t)=Φ(t)x(0),0tT
(3)
For tT, it can be calculated by
x(t+sT)=Φ(t)Φs(T)x(0)
(4)
where s = 1,2,3…. and Φ(T) is the Floquet transition matrix (FTM) or the Monodromy matrix. The stability criteria for periodic systems depend upon the eigenvalues of Φ(T), called the Floquet multipliers, and the system is stable only if all the Floquet multipliers lie on or inside the unit circle.
According to the Lyapunov–Floquet theorem, the STM can be expressed as
Φ(t)=Q(t)eRt,RRn×n,t0
(5)
where Q(t)Rn×n is known as the Lyapunov–Floquet (L-F) transformation matrix. Due to the periodic coefficients in A(t), a direct application of invariant manifold is not feasible. The requisite transformation of the linear matrix to a nonautonomous one was first accomplished via L-F transformation
x=Q(t)z
(6)
The elements of the L-F transformation matrix (Q(t)) are represented as a truncated Fourier series [33]. Applying the L-F transformation to Eq. (2) results in Eq. (7) with a time-invariant linear part, as detailed in Ref. [34]
z˙=Rz
(7)
where R is a constant n × n matrix that generally tends to be complex. The eigenvalues of the R matrix are known as the Floquet exponents, which can also be used as an indicator for the stability of the dynamical system. If the Floquet exponents’ real parts are positive, it is regarded as an unstable system, and if they are negative, the system is considered asymptotically stable. If the Floquet exponents are purely imaginary, the system is considered to be marginally stable.

2.3 Time-Delayed Embedding.

The concept of time-delayed embedding is explained briefly in this subsection. A comprehensive explanation of this technique is provided in Refs. [35,36]. The time-delayed embedding approach has been applied to time series processing algorithms, predictions, noise, system identification, chaos synchronization, and control [37]. In the case of an unknown dynamical system, the time-delayed embedding aids in the reconstruction of system states and dynamics from the observations or the system measurements over time. Here, the system’s true state wt is assumed to lie on the attractor A ⊂ ℝk in a new phase space of dimension k. The observed time series data (dt=f(wt)) is considered to be a smooth (continuously differentiable) map from the phase space to the scalar value (f:ℝk → ℝ). A time-delayed reconstruction in m dimensions with time delay τ would result in a vector dtRm i.e., dt=(dt,dt+τ,dt+2τ,,dt+(m1)τ). A one to one map from the attractor A to the reconstruction space ℝm is embedded and preserves the dynamics.

For the optimal value τ, the first minima of the average mutual information (AMI) function is considered. AMI evaluates the amount of information shared between two data sets over a broad range of time delays [38]. The Takens theorem suggests that the value of the embedded dimension maintains, m > 2 × dimension(A). However, it is not feasible to wait for as many data points for fast or real-time applications, delaying the control output. Therefore, a global false nearest neighbor (GFNN) analysis is performed to compute the optimal value of m [39]. In this analysis, the distances between the neighboring trajectories at successively higher dimensions are compared, and when the trajectories overlap, false neighbors are identified. As the dimension increases, the total percentage of false neighbors declines, and m is selected when this percentage approaches zero [39].

The application of the techniques, as mentioned earlier, toward the analysis of human walking gait pattern is detailed in Sec. 3.

3 Methodology

In this section, the sequence of operations performed to study the human walking gait is detailed. The time series of kinematic data collected from healthy subjects are mapped to a single degree of freedom (SDOF) oscillator. Then it is mapped back to verify the accurate representation of the original system. A detailed description of the methodology adopted is described in the Secs. 3.13.3.

3.1 Data Collection.

The normal walking tests were performed on seven human subjects. These subjects were selected to have had no previous neurological or physical injury to ensure normal gait. The physical parameters measured for each subject are listed in Table 1. The test protocol was created as per the Institutional Review Board at Arizona State University (Study 00009416). The kinematic and dynamic data were collected for all subjects using a ten camera Vicon Motion capture system (100 Hz sampling rate) and a Bertec instrumented treadmill with independent force measuring belts for each foot (1.75 × 0.5 m each). The marker setup used was based on the Newington-Helen Hayes model. The subjects walked on the treadmill at 1.2 m/s for 3 min. The GRF data were collected from the instrumented treadmill at 1000 Hz but was downsampled later to 100 Hz to improve computation time.

Table 1

Physical parameters for the subjects involved

SubjectAge (years)Weight (kg)Height (m)
143108.11.82
25065.41.70
34786.81.77
42175.61.77
52882.191.77
62473.91.80
71976.921.77
SubjectAge (years)Weight (kg)Height (m)
143108.11.82
25065.41.70
34786.81.77
42175.61.77
52882.191.77
62473.91.80
71976.921.77

The lower body joint angles were measured relative to the standing position. The output of the Vicon Nexus system and Bertec system were postprocessed using matlab 2018a. The noise was filtered out from the collected data using a low-pass Butterworth filter with a cutoff frequency of 10 Hz. The in-built plug-ins provided the instances of Heel-strike and Toe-off. The first and last 30 s of the joint kinematic data were discarded to account for the acceleration and deceleration artifacts. The sagittal, coronal, and transverse plane joint angles data were initially analyzed, and later the same approach was applied to the dynamic data of GRF.

3.2 Mapping and Re-Mapping.

As mentioned earlier, the kinematic data of lower extremity joint angles were initially considered for analysis in all three planes. The time series of joint angle data were differentiated with respect to time, using a linear approximation to compute the corresponding angular rate. Both the joint angle and angular rate time series were normalized by scaling with respect to their maximum value and then shifted such that the start point of the time series was either 0 or 1. The normalization of these data resulted in an accurate comparison between the data. A matrix x(t)2×n was constructed using the shifted joint angle time series and its derivative. In this matrix, each column corresponded to the state vector at each instant.

As discussed in Sec. 2.3, the optimal time delay (τ) for the shifted joint angle time series was calculated using the first minimum of the AMI function in matlab [40]. The optimal embedding dimension (m) was calculated using the GFNN analysis specified in Sec. 2.3. A new matrix (X2×m) was created using the values of m and τ, as follows:
X(t)=[x(t),x(t+τ),x(t+2τ),,x(t+(m1)τ)]
(8)
This matrix was generated from the initial point (t = 0) to the average principal period (t = T) of the time series (X(t)). These reconstructed matrices were utilized to evaluate the FTM (Φ(T)) of the system, using the Eq. (3) as indicated below
Φ(T)=X(T)X(0)
(9)
Since X(0) was not a square matrix, a pseudo-inverse (Moore–Penrose inverse), X(0) was utilized. Using this relation in Eq. (5), the constant matrix R was computed at the principal time period (T) using the relation
R=logm(Φ(T))T
(10)
The resulting 2 × 2 matrix (R) was generally complex in nature, and its eigenvalues, known as Floquet exponents (λ) was utilized to generate a time-invariant system, expressed as an SDOF harmonic oscillator, in Eq. (11)
z¨+||λ||2z=0
(11)

This SDOF harmonic oscillator formed an invariant manifold dynamical equation for the original system. The STM (Φ(t)) was computed from the relationship in the Eq. (3) and was further utilized to calculate an analytical expression for Q(t) using Eq. (5). The SDOF system was numerically integrated to find a solution in the z-domain. Using the relationship in Eq. (7), the solution in the z-domain is back-transformed/re-mapped to the original coordinates in the x-domain. The re-mapped solutions were compared with the time series in the original coordinates for validation.

3.3 Reinitialization Procedure.

While comparing the re-mapped solution to the original time series, it was observed that the correlation between both the data was good for the first few gait cycles, but they slowly started diverging later. This behavior indicated that a reinitialization operation was essential to update the system's initial conditions after every gait cycle. This reinitialization was achieved by re-integrating the SDOF system after every gait cycle. The initial conditions for the integration were obtained from the following relation
z(t)=Q1(t)x(t)
(12)
where t′ was the time at the end of any given gait cycle and resulted in improved correlation to the original time series data. The effectiveness of the reinitialization operation is demonstrated in Sec. 4.

4 Results and Discussion

The above-discussed approach was applied to the kinematic and kinetic data for all the subjects. The kinematic data considered were the lower body joint angles (hip, knee, and ankle) in all three planes (sagittal, coronal, and transverse), whereas the kinetic data considered were vertical and horizontal GRFs. Figure 1 shows the reconstruction operation, which was similar for all the cases. The time-invariant manifold was created using the system in Eq. (7). Table 2 lists the exponent values, i.e., the eigenvalues of the R matrix, for the sagittal plane data of the three joints. Appendix  A contains the other two planes (coronal and transverse) data and the two GRF cases.

Fig. 1
Re-mapping operation for the normalized knee joint angle for all three planes for subject 4: (a) sagittal plane original system, (b) sagittal plane manifold system, (c) sagittal plane reconstructed system, (d) coronal plane original system, (e) coronal plane manifold system, (f) coronal plane reconstructed system, (g) transverse plane original system, (h) transverse plane manifold system, and (i) transverse plane reconstructed system
Fig. 1
Re-mapping operation for the normalized knee joint angle for all three planes for subject 4: (a) sagittal plane original system, (b) sagittal plane manifold system, (c) sagittal plane reconstructed system, (d) coronal plane original system, (e) coronal plane manifold system, (f) coronal plane reconstructed system, (g) transverse plane original system, (h) transverse plane manifold system, and (i) transverse plane reconstructed system
Table 2

Exponent values for the joint angle data in the sagittal plane

SubjectHipKneeAnkle
10.0004 ± 0.0437i−0.0004 ± 0.0073i0.0002 ± 0.0155i
2−0.0041 ± 0.0978i−0.0072 ± 0.0159i−0.0026 ± 0.0089i
30.0001 ± 0.0517i0.0032 ± 0.0467i−0.0323 ± 0.0711i
40.0001 ± 0.0007i−0.0005 ± 0.0705i0.0003 ± 0.0002i
50 ± 0.0144i−0.0001 ± 0.0725i0.0032 ± 0.0045i
60 ± 0.036i0.0003 ± 0.0188i−0.097 ± 0i
70.0001 ± 0.0587i0.0008 ± 0.0062i0.0035 ± 0i
SubjectHipKneeAnkle
10.0004 ± 0.0437i−0.0004 ± 0.0073i0.0002 ± 0.0155i
2−0.0041 ± 0.0978i−0.0072 ± 0.0159i−0.0026 ± 0.0089i
30.0001 ± 0.0517i0.0032 ± 0.0467i−0.0323 ± 0.0711i
40.0001 ± 0.0007i−0.0005 ± 0.0705i0.0003 ± 0.0002i
50 ± 0.0144i−0.0001 ± 0.0725i0.0032 ± 0.0045i
60 ± 0.036i0.0003 ± 0.0188i−0.097 ± 0i
70.0001 ± 0.0587i0.0008 ± 0.0062i0.0035 ± 0i

The reinitialization process was necessary for a well-correlated reconstruction. As mentioned in Sec. 3.3, the data diverged at the end. This phenomenon can be seen clearly in the example presented in Fig. 2. The error increased with time in Fig. 2(a), while the amount of error was manageable in Fig. 2(b). Hence, utilizing the reinitialization technique, the amount of error over multiple gait cycles was mitigated.

Fig. 2
Normalized original GRF data (dashed line) compared with the reconstructed data (solid line) with and without reinitialization for subject 6 (samples versus normalized force/weight): (a) without re-initialization and (b) with re-initialization
Fig. 2
Normalized original GRF data (dashed line) compared with the reconstructed data (solid line) with and without reinitialization for subject 6 (samples versus normalized force/weight): (a) without re-initialization and (b) with re-initialization

We opted to showcase the comparison between data for only two subjects to keep the article shorter. Figure 3 shows a comparison between the normalized joint angle data and the reconstructed data using time-invariant manifolds. Figure 4 shows the comparison between the normalized GRF data and the reconstructed data for subject 4. Similar information is presented in Appendix  C for subject 5. The real parts of the exponents (from Table 2 and Appendix  A) are either negative or close to zero, denoting that the data is bounded. As the normal gait is, by definition, stable, and the exponent data is in support of this observation.

Fig. 3
Normalized joint angle data (dashed line) compared with the reconstructed angular data (solid line) for subject 4 (samples versus normalized angle in radians): (a) sagittal hip angle (R2 = 0.999), (b) coronal hip angle (R2 = 0.998), (c) transverse hip angle (R2 = 0.999), (d) sagittal knee angle (R2 = 0.999), (e) coronal knee angle (R2 = 0.998), (f) transverse knee angle (R2 = 0.997), (g) sagittal ankle angle (R2 = 0.995), (h) coronal ankle angle (R2 = 0.999), and (i) transverse ankle angle (R2 = 0.975)
Fig. 3
Normalized joint angle data (dashed line) compared with the reconstructed angular data (solid line) for subject 4 (samples versus normalized angle in radians): (a) sagittal hip angle (R2 = 0.999), (b) coronal hip angle (R2 = 0.998), (c) transverse hip angle (R2 = 0.999), (d) sagittal knee angle (R2 = 0.999), (e) coronal knee angle (R2 = 0.998), (f) transverse knee angle (R2 = 0.997), (g) sagittal ankle angle (R2 = 0.995), (h) coronal ankle angle (R2 = 0.999), and (i) transverse ankle angle (R2 = 0.975)
Fig. 4
Normalized original GRF data (dashed line) compared with the reconstructed data (solid line) for subject 4 (samples versus normalized force/weight): (a) horizontal GRF (R2 = 0.999) and (b) vertical GRF (R2 = 0.989)
Fig. 4
Normalized original GRF data (dashed line) compared with the reconstructed data (solid line) for subject 4 (samples versus normalized force/weight): (a) horizontal GRF (R2 = 0.999) and (b) vertical GRF (R2 = 0.989)

The r2 error value was calculated to quantify the error for each case. Table 3 lists the r2 values for the error between the hip joint angle data and the reconstructed data. Appendix  B contains the r2 values for the other cases. The lowest r2 value seen among all the data was 0.8332 (coronal plane ankle joint angle). Oh et al. [41] used a neural network framework to predict GRF for 48 subjects. Their experiments resulted in the lowest r value of 0.841 (r2 = 0.7073). Ren et al. [42] created a model to estimate GRF. Oh et al. applied the methods presented in Ref. [42] to obtain the lowest r value of 0.677 (r2 = 0.4583). The method presented in the article yielded a better correlation than both these methods, showing that the algorithm is useful in reconstructing the data with high accuracy. It should also be noted that the Floquet theory method requires less computational power than neural networks.

Table 3

R-squared error value of the original hip joint angle versus the reconstructed data

SubjectSagittalCoronalTransverse
11.00000.99901.0000
20.99801.00000.9940
30.99900.99801.0000
41.00001.00000.9430
51.00001.00000.9980
60.99900.99900.9990
70.99900.98700.9310
SubjectSagittalCoronalTransverse
11.00000.99901.0000
20.99801.00000.9940
30.99900.99801.0000
41.00001.00000.9430
51.00001.00000.9980
60.99900.99900.9990
70.99900.98700.9310

Since the complex dynamical systems describing the gait kinematics were transformed into a simple SDOF system, creating a time-dependent controller would be easier. The discussed algorithm can be utilized in mobile wearable devices to estimate kinematic data. Since the algorithm is based on previously collected data, the gait data can be estimated with a high degree of accuracy under similar conditions. The algorithm can be applied in active prostheses and orthoses. Utilizing the algorithm in prostheses design and control, such as those used in Ref. [9] or [16], could lead to better performance. Real-time applications for the algorithm are being researched and developed.

The results demonstrate the successful transformation of the time series data to an SDOF harmonic oscillator system that would essentially serve as an invariant manifold. However, this SDOF harmonic oscillator approach can be extended further toward system identification. For the cases where Q1(t) exists or Q(t) is invertible and differentiable, one could obtain the analytical expression for the original dynamical system using the expression below
x˙=[Q˙(t)+Q(t)R]Q1(t)x
(13)
Equation (13) can also be rewritten as
x˙=M(t)x
(14)
where M(t)=[Q˙(t)+Q(t)R]Q1(t) and by comparing Eq. (12) with Eq. (2), it is observed that M(t) is equivalent to A(t) in the original system coordinates. In the event of invertibility issues for the Q(t) matrix, the ZNN approach could compute the inverse [43]. From the elements of M(t), the dynamical system parameters (such as time-varying stiffness, mass, and damping) for a human walking gait system could be estimated, aiding in designing assistive devices for rehabilitation and humanoid robots efficiently. The extension of the discussed approach toward the system identification for human walking gait will be detailed in future work.

5 Conclusion

In this work, the authors have demonstrated the application of Floquet theory toward the analysis of human walking gait behavior. The kinematic and dynamic features of human walking gait were considered to be time-periodic systems and were validated using the data collected from seven healthy subjects. A time-delayed embedding of the observed time series, with the aid of AMI and GFNN, was proven to preserve the dynamics of the original system, facilitating the application of Floquet theory, and the human walking gait was observed to be a marginally stable time-periodic system from the Floquet multiplier and exponents.

This work also demonstrated a conformal mapping approach that successfully converted the time series of the human walking gait features to an SDOF harmonic oscillator system. A re-mapping of the oscillator system, with the aid of a reinitialization procedure, was successful in reproducing the system state variations close to the observed time series. In this work, the re-mapping approach was validated for both kinematic and dynamic features of walking gait.

The time-invariant form could be used in building a controller for a phase-based-oscillator for exoskeletons, rehabilitation devices, or humanoid robots. As mentioned earlier, the system identification for human walking gait and its application toward efficient control design will be demonstrated in the forthcoming paper. It is important to note that we can predict at least two gait cycle curves accurately into the future with this method.

Acknowledgment

The authors would like to acknowledge the Arizona State University for their support.

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. Data provided by a third party listed in Acknowledgment.

Nomenclature

     
  • t =

    time

  •  
  • m =

    embedding dimension

  •  
  • T =

    principal time period

  •  
  • x(t) =

    the original system

  •  
  • A(t) =

    linear system matrix

  •  
  • E(t) =

    nonlinear system matrix

  •  
  • G(t) =

    external forcing term

  •  
  • Q(t) =

    L-F transformation matrix

  •  
  • X(t) =

    time-delay embedded system

  •  
  • x˙(t) =

    the original system derivative

  •  
  • λ =

    Floquet exponents

  •  
  • τ =

    time delay

  •  
  • ϕ(t) =

    state transition matrix

  •  
  • Φ(T) =

    Floquet transition matrix

Appendix A: Exponent Values

Table 4

Exponent values for the joint angle data in the coronal plane

SubjectHipKneeAnkle
−0.0229 ± 0.0147i −0.0029 ± 0.0096i 0.0004 ± 0.0189i 
0.0001 ± 0.0033i −0.1331 ± 0i −0.1810 ± 0i 
0.0026 ± 0.0395i −0.0160 ± 0.0554i 0.005 ± 0i 
0.0002 ± 0.0001i −0.0310 ± 0.2015i −0.0002 ± 0.0017i 
0.0001 ± 0i 0.0000 ± 0.0029i −0.0421 ± 0.0965i 
0.0006 ± 0.0420i −0.0041 ± 0.1121i −0.5312 ± 0.1532i 
−0.0653 ± 0.1293i −0.0002 ± 0.0141i −0.0377 ± 0i 
SubjectHipKneeAnkle
−0.0229 ± 0.0147i −0.0029 ± 0.0096i 0.0004 ± 0.0189i 
0.0001 ± 0.0033i −0.1331 ± 0i −0.1810 ± 0i 
0.0026 ± 0.0395i −0.0160 ± 0.0554i 0.005 ± 0i 
0.0002 ± 0.0001i −0.0310 ± 0.2015i −0.0002 ± 0.0017i 
0.0001 ± 0i 0.0000 ± 0.0029i −0.0421 ± 0.0965i 
0.0006 ± 0.0420i −0.0041 ± 0.1121i −0.5312 ± 0.1532i 
−0.0653 ± 0.1293i −0.0002 ± 0.0141i −0.0377 ± 0i 
Table 5

Exponent values for the Joint angle data in the transverse plane

SubjectHipKneeAnkle
0.0005 ± 0.0078i −0.5846 ± 0.4916i −0.0001 ± 0.0249i 
−0.0598 ± 0.1300i −0.0290 ± 0.0567i −0.4669 ± 0i 
−0.0006 ± 0.0016i −0.0043 ± 0.0976i −0.0905 ± 0i 
−0.1817 ± 0.1344i −0.0750 ± 0.2112i −0.0001 ± 0.0018i 
0.0119 ± 0.0720i −0.0001 ± 0.0057i 0.0044 ± 0.0034i 
0.0262 ± 0.0545i 0.0051 ± 0.0960i −0.5866 ± 0.1358i 
−0.3306 ± 0.1093i −0.0033 ± 0.0071i −0.0368 ± 0.0210i 
SubjectHipKneeAnkle
0.0005 ± 0.0078i −0.5846 ± 0.4916i −0.0001 ± 0.0249i 
−0.0598 ± 0.1300i −0.0290 ± 0.0567i −0.4669 ± 0i 
−0.0006 ± 0.0016i −0.0043 ± 0.0976i −0.0905 ± 0i 
−0.1817 ± 0.1344i −0.0750 ± 0.2112i −0.0001 ± 0.0018i 
0.0119 ± 0.0720i −0.0001 ± 0.0057i 0.0044 ± 0.0034i 
0.0262 ± 0.0545i 0.0051 ± 0.0960i −0.5866 ± 0.1358i 
−0.3306 ± 0.1093i −0.0033 ± 0.0071i −0.0368 ± 0.0210i 
Table 6

Exponent values for the ground reaction force data

SubjectHorizontalVertical
−0.0113 ± 0.0089i −0.0006 ± 0.0153i 
−0.0676 ± 0.0727i −0.0105 ± 0.0336i 
−0.0174 ± 0.0172i 0.0010 ± 0.0086i 
−0.0170 ± 0.0917i −0.0050 ± 0.0747i 
−0.0042 ± 0.1816i −0.0093 ± 0.1409i 
0.0307 ± 0.1917i −0.0606 ± 0.0772i 
−0.0320 ± 0.0425i −0.0052 ± 0.0246i 
SubjectHorizontalVertical
−0.0113 ± 0.0089i −0.0006 ± 0.0153i 
−0.0676 ± 0.0727i −0.0105 ± 0.0336i 
−0.0174 ± 0.0172i 0.0010 ± 0.0086i 
−0.0170 ± 0.0917i −0.0050 ± 0.0747i 
−0.0042 ± 0.1816i −0.0093 ± 0.1409i 
0.0307 ± 0.1917i −0.0606 ± 0.0772i 
−0.0320 ± 0.0425i −0.0052 ± 0.0246i 

Appendix B: R-Squared Values

Table 7

R-squared error value of original knee joint angle versus the reconstructed data

SubjectSagittalCoronalTransverse
1.0000 0.9999 0.9999 
1.0000 0.9836 0.9995 
0.9999 0.9995 0.9979 
0.9998 0.9876 0.9992 
0.9994 1.0000 1.0000 
0.9998 0.9978 1.0000 
1.0000 0.9999 1.0000 
SubjectSagittalCoronalTransverse
1.0000 0.9999 0.9999 
1.0000 0.9836 0.9995 
0.9999 0.9995 0.9979 
0.9998 0.9876 0.9992 
0.9994 1.0000 1.0000 
0.9998 0.9978 1.0000 
1.0000 0.9999 1.0000 
Table 8

R-squared error value of original ankle joint angle versus the reconstructed data

SubjectSagittalCoronalTransverse
0.9999 0.9999 0.9993 
0.9995 0.8450 0.8717 
0.9958 0.9990 0.9917 
0.9999 0.9999 0.9999 
0.9999 0.9990 0.9998 
0.9897 0.8332 0.9997 
0.9997 0.9898 0.9975 
SubjectSagittalCoronalTransverse
0.9999 0.9999 0.9993 
0.9995 0.8450 0.8717 
0.9958 0.9990 0.9917 
0.9999 0.9999 0.9999 
0.9999 0.9990 0.9998 
0.9897 0.8332 0.9997 
0.9997 0.9898 0.9975 
Table 9

R-squared error value of original ground reaction force data versus the reconstructed data

SubjectHorizontalVertical
0.9999 0.9999 
0.9995 0.9999 
0.9996 0.9984 
0.9983 0.9843 
0.9974 0.9977 
0.9972 0.9999 
0.9999 0.9999 
SubjectHorizontalVertical
0.9999 0.9999 
0.9995 0.9999 
0.9996 0.9984 
0.9983 0.9843 
0.9974 0.9977 
0.9972 0.9999 
0.9999 0.9999 

Appendix C: Comparative Images for Subject 5

Fig. 5

Normalized joint angle data (dashed line) compared with the reconstructed angular data (solid line) for subject 5 (samples versus normalized angle in radians): (a) sagittal hip angle (R2 = 1), (b) coronal hip angle (R2 = 1), (c) transverse hip angle (R2 = 0.997), (d) sagittal knee angle (R2 = 0.999), (e) coronal knee angle (R2 = 0.986), (f) transverse knee angle (R2 = 0.999), (g) sagittal ankle angle (R2 = 1), (h) coronal ankle angle (R2 = 1), and (i) transverse ankle angle (R2 = 0.999)

Fig. 5

Normalized joint angle data (dashed line) compared with the reconstructed angular data (solid line) for subject 5 (samples versus normalized angle in radians): (a) sagittal hip angle (R2 = 1), (b) coronal hip angle (R2 = 1), (c) transverse hip angle (R2 = 0.997), (d) sagittal knee angle (R2 = 0.999), (e) coronal knee angle (R2 = 0.986), (f) transverse knee angle (R2 = 0.999), (g) sagittal ankle angle (R2 = 1), (h) coronal ankle angle (R2 = 1), and (i) transverse ankle angle (R2 = 0.999)

Fig. 6

Normalized original GRF data (dashed line) compared with the reconstructed data (solid line) for subject 5 (samples versus normalized force/weight): (a) horizontal GRF (R2 = 0.998) and (b) vertical GRF (R2 = 0.984)

Fig. 6

Normalized original GRF data (dashed line) compared with the reconstructed data (solid line) for subject 5 (samples versus normalized force/weight): (a) horizontal GRF (R2 = 0.998) and (b) vertical GRF (R2 = 0.984)

References

1.
Chinimilli
,
P. T.
,
Redkar
,
S.
, and
Sugar
,
T.
,
2019
, “
A Two-Dimensional Feature Space-Based Approach for Human Locomotion Recognition
,”
IEEE Sens. J.
,
19
(
11
), pp.
4271
4282
.
2.
Hurmuzlu
,
Y.
, and
Basdogan
,
C.
,
1994
, “
On the Measurement of Dynamic Stability of Human Locomotion
,”
J. Biomech. Eng.
,
116
(
1
), pp.
30
36
.
3.
Hurmuzlu
,
Y.
,
Basdogan
,
C.
, and
Stoianovici
,
D.
,
1996
, “
Kinematics and Dynamic Stability of the Locomotion of Post-Polio Patients
,”
Journal of Biomechanical Engineering
,
118
(
3
), pp.
405
411
.
4.
Rosenstein
,
M. T.
,
Collins
,
J. J.
, and
De Luca
,
C. J.
,
1993
, “
A Practical Method for Calculating Largest Lyapunov Exponents From Small Data Sets
,”
Phys. D Nonlinear Phenom.
,
65
(
1–2
), pp.
117
134
.
5.
Dingwell
,
J. B.
, and
Cusumano
,
J. P.
,
2000
, “
Nonlinear Time Series Analysis of Normal and Pathological Human Walking
,”
Chaos An Interdiscip. J. Nonlinear Sci.
,
10
(
4
), pp.
848
863
.
6.
Kang
,
H. G.
, and
Dingwell
,
J. B.
,
2008
, “
Effects of Walking Speed, Strength and Range of Motion on Gait Stability in Healthy Older Adults
,”
J. Biomech.
,
41
(
14
), pp.
2899
2905
.
7.
Bruijn
,
S. M.
,
Bregman
,
D. J. J.
,
Meijer
,
O. G.
,
Beek
,
P. J.
, and
van Dieën
,
J. H.
,
2012
, “
Maximum Lyapunov Exponents as Predictors of Global Gait Stability: A Modelling Approach
,”
Med. Eng. Phys.
,
34
(
4
), pp.
428
436
.
8.
Wagenaar
,
R. C.
, and
van Emmerik
,
R. E. A.
,
1994
, “
Dynamics of Pathological Gait
,”
Hum. Mov. Sci.
,
13
(
3–4
), pp.
441
471
.
9.
De La Fuente
,
J.
,
Subramanian
,
S. C.
,
Sugar
,
T. G.
, and
Redkar
,
S.
,
2020
, “
A Robust Phase Oscillator Design for Wearable Robotic Systems
,”
Rob. Auton. Syst.
,
128
, p.
103514
.
10.
Bhounsule
,
P. A.
, and
Zamani
,
A.
,
2017
, “
A Discrete Control Lyapunov Function for Exponential Orbital Stabilization of the Simplest Walker
,”
ASME J. Mech. Robot.
,
9
(
5
), p.
051011
.
11.
la Fuente
,
J.
,
Sugar
,
T. G.
, and
Redkar
,
S.
,
2017
, “
Nonlinear, Phase-Based Oscillator to Generate and Assist Periodic Motions
,”
ASME J. Mech. Robot.
,
9
(
2
), p.
024502
.
12.
Várkonyi
,
P. L.
,
2015
, “
On the Stability of Rigid Multibody Systems With Applications to Robotic Grasping and Locomotion
,”
ASME J. Mech. Robot.
,
7
(
4
).
13.
Hao
,
M.
,
Zhang
,
J.
,
Chen
,
K.
,
Asada
,
H.
, and
Fu
,
C.
,
2020
, “
Supernumerary Robotic Limbs to Assist Human Walking With Load Carriage
,”
ASME J. Mech. Robot.
,
12
(
6
), p.
061014
.
14.
Spyrakos-Papastavridis
,
E.
,
Dai
,
J. S.
,
Childs
,
P.
, and
Tsagarakis
,
N. G.
,
2018
, “
Selective-Compliance-Based Lagrange Model and Multilevel Noncollocated Feedback Control of a Humanoid Robot
,”
ASME J. Mech. Robot.
,
10
(
3
), p.
031009
.
15.
Luxman
,
R.
, and
Zielinska
,
T.
,
2017
, “
Robot Motion Synthesis Using Ground Reaction Forces Pattern: Analysis of Walking Posture
,”
Int. J. Adv. Robot. Syst.
,
14
(
4
), p.
1729881417720873
.
16.
Azimi
,
V.
,
Nguyen
,
T. T.
,
Sharifi
,
M.
,
Fakoorian
,
S. A.
, and
Simon
,
D.
,
2018
, “
Robust Ground Reaction Force Estimation and Control of Lower-Limb Prostheses: Theory and Simulation
,”
IEEE Trans. Syst. Man, Cybern. Syst.
,
50
(
8
), pp.
3024
3035
.
17.
Almeida
,
L.
,
Santos
,
V.
, and
Silva
,
F.
,
2018
, “
A Novel Wireless Instrumented Shoe for Ground Reaction Forces Analysis in Humanoids
,”
2018 IEEE International Conference on Autonomous Robot Systems and Competitions (ICARSC)
,
Torres Vedras, Portugal
,
Apr. 25–27
, pp.
36
41
.
18.
Wu
,
X. A.
,
Huh
,
T. M.
,
Mukherjee
,
R.
, and
Cutkosky
,
M.
,
2016
, “
Integrated Ground Reaction Force Sensing and Terrain Classification for Small Legged Robots
,”
IEEE Robot. Autom. Lett.
,
1
(
2
), pp.
1125
1132
.
19.
Baumgart
,
C.
,
Schubert
,
M.
,
Hoppe
,
M. W.
,
Gokeler
,
A.
, and
Freiwald
,
J.
,
2017
, “
Do Ground Reaction Forces During Unilateral and Bilateral Movements Exhibit Compensation Strategies Following ACL Reconstruction?
,”
Knee Surg. Sport. Traumatol. Arthrosc.
,
25
(
5
), pp.
1385
1394
.
20.
Baumgart
,
C.
,
Hoppe
,
M. W.
, and
Freiwald
,
J.
,
2017
, “
Phase-Specific Ground Reaction Force Analyses of Bilateral and Unilateral Jumps in Patients With ACL Reconstruction
,”
Orthop. J. Sport. Med.
,
5
(
6
), p.
2325967117710912
.
21.
Lanier
,
A. S.
,
Knarr
,
B. A.
,
Stergiou
,
N.
,
Snyder-Mackler
,
L.
, and
Buchanan
,
T. S.
,
2020
, “
ACL Injury and Reconstruction Affect Control of Ground Reaction Forces Produced During a Novel Task That Simulates Cutting Movements
,”
J. Orthop. Res.
,
38
(
8
), pp.
1746
1752
.
22.
Anand
,
M.
,
Seipel
,
J.
, and
Rietdyk
,
S.
,
2017
, “
A Modelling Approach to the Dynamics of Gait Initiation
,”
J. R. Soc. Interface
,
14
(
128
), p.
20170043
.
23.
Chicone
,
C.
,
2006
,
Ordinary Differential Equations with Applications
,
Springer Science & Business Media
,
New York
.
24.
Morris
,
B.
, and
Grizzle
,
J. W.
,
2009
, “
Hybrid Invariant Manifolds in Systems With Impulse Effects With Application to Periodic Locomotion in Bipedal Robots
,”
IEEE Trans. Automat. Contr.
,
54
(
8
), pp.
1751
1764
.
25.
Frank
,
J.
,
Mannor
,
S.
, and
Precup
,
D.
,
2010
, “
Activity and Gait Recognition With Time-Delay Embeddings
,”
Twenty-Fourth AAAI Conference on Artificial Intelligence
,
Atlanta, GA
,
July 11–15
.
26.
Hidaka
,
S.
, and
Fujinami
,
T.
,
2013
, “
Topological Similarity of Motor Coordination in Rhythmic Movements
,”
Proceedings of the Annual Meeting of the Cognitive Science Society
,
35
(
35
), pp.
2548
2553
.
27.
Perc
,
M.
,
2005
, “
The Dynamics of Human Gait
,”
Eur. J. Phys.
,
26
(
3
), p.
525
.
28.
Miller
,
D. J.
,
Stergiou
,
N.
, and
Kurz
,
M. J.
,
2006
, “
An Improved Surrogate Method for Detecting the Presence of Chaos in Gait
,”
J. Biomech.
,
39
(
15
), pp.
2873
2876
.
29.
Scafetta
,
N.
,
Marchi
,
D.
, and
West
,
B. J.
,
2009
, “
Understanding the Complexity of Human Gait Dynamics
,”
Chaos An Interdiscip. J. Nonlinear Sci.
,
19
(
2
), p.
26108
.
30.
Decker
,
L. M.
,
Cignetti
,
F.
, and
Stergiou
,
N.
,
2010
, “
Complexity and Human Gait
,”
Rev. Andaluza Med. del Deport.
,
3
(
1
), pp.
2
12
.
31.
Sinha
,
S. C.
,
Pandiyan
,
R.
, and
Bibb
,
J. S.
,
1996
, “
Liapunov-Floquet Transformation: Computation and Applications to Periodic Systems
,”
Journal of Vibration and Acoustics
,
118
(
2
), pp.
209
219
.
32.
Redkar
,
S.
,
2012
, “
Lyapunov Stability of Quasiperiodic Systems
,”
Math. Probl. Eng.
,
2012
.
33.
Sinha
,
S. C.
, and
Pandiyan
,
R.
,
1994
, “
Analysis of Quasilinear Dynamical Systems With Periodic Coefficients via Liapunov-Floquet Transformation
,”
Int. J. Non. Linear. Mech.
,
29
(
5
), pp.
687
702
.
34.
Hale
,
J. K.
,
2009
,
Ordinary Differential Equations
,
Dover Publications
,
Mineola, NY
.
35.
Kantz
,
H.
, and
Schreiber
,
T.
,
2004
,
Nonlinear Time Series Analysis
,
Cambridge University Press
,
Cambridge, UK
.
36.
Takens
,
F.
,
1981
, “Detecting Strange Attractors in Turbulence,”
Dynamical Systems and Turbulence, Warwick 1980
,
Springer
,
New York
, pp.
366
381
.
37.
Eftekhari
,
A.
,
Yap
,
H. L.
,
Wakin
,
M. B.
, and
Rozell
,
C. J.
,
2018
, “
Stabilizing Embedology: Geometry-Preserving Delay-Coordinate Maps
,”
Phys. Rev. E
,
97
(
2
), p.
22222
.
38.
Fraser
,
A. M.
, and
Swinney
,
H. L.
,
1986
, “
Independent Coordinates for Strange Attractors From Mutual Information
,”
Phys. Rev. A
,
33
(
2
), p.
1134
.
39.
Kennel
,
M. B.
,
Brown
,
R.
, and
Abarbanel
,
H. D. I.
,
1992
, “
Determining Embedding Dimension for Phase-Space Reconstruction Using a Geometrical Construction
,”
Phys. Rev. A
,
45
(
6
), p.
3403
.
40.
Leontitsis
,
A.
,
2019
,
Mutual Average Information, MATLAB Central File Exchange
, https://www.mathworks.com/matlabcentral/fileexchange/880-mutual-average-information,
Accessed March 2019
.
41.
Oh
,
S. E.
,
Choi
,
A.
, and
Mun
,
J. H.
,
2013
, “
Prediction of Ground Reaction Forces During Gait Based on Kinematics and a Neural Network Model
,”
J. Biomech.
,
46
(
14
), pp.
2372
2380
.
42.
Ren
,
L.
,
Jones
,
R. K.
, and
Howard
,
D.
,
2008
, “
Whole Body Inverse Dynamics Over a Complete Gait Cycle Based Only on Measured Kinematics
,”
J. Biomech.
,
41
(
12
), pp.
2750
2759
.
43.
Zhang
,
Y.
,
Yi
,
C.
, and
Ma
,
W.
,
2009
, “
Simulation and Verification of Zhang Neural Network for Online Time-Varying Matrix Inversion
,”
Simul. Model. Pract. Theory
,
17
(
10
), pp.
1603
1617
.