Bio-Inspired Real-Time Prediction of Human Locomotion for Exoskeletal Robot Control

Human motion detection is of fundamental importance for control of human–robot coupled systems such as exoskeletal robots. Inertial measurement units have been widely used for this purpose, although delay is a major challenge for inertial measurement unit-based motion capture systems. In this paper, we use previous and current inertial measurement unit readings to predict human locomotion based on their kinematic properties. Human locomotion is a synergetic process of the musculoskeletal system characterized by smoothness, high nonlinearity, and quasi-periodicity. Takens’ reconstruction method can well characterize quasi-periodicity and nonlinear systems. With Takens’ reconstruction framework, we developed improving methods, including Gaussian coefficient weighting and offset correction (which is based on the smoothness of human locomotion), Kalman fusion with complementary joint data prediction and united source of historical embedding generation (which is synergy-inspired), and Kalman fusion with the Newton-based method with a velocity and acceleration high-gain observer (also based on smoothness). After thorough analysis of the parameters and the effect of these improving techniques, we propose a novel prediction method that possesses the combined advantages of parameter robustness, high accuracy, trajectory smoothness, zero dead time, and adaptability to irregularities. The proposed methods were tested and validated by experiments, and the real-time applicability in a human locomotion capture system was also validated.


Introduction
Human motion detection is of fundamental importance for control of human-robot coupled systems such as exoskeletal robots, and for many other applications.In recent years, inertial measurement unit (IMU) technologies have quickly evolved, and complex mechanisms for tracking the operator's movement [1,2] have gradually been replaced by motion capture systems comprising IMUs [3][4][5][6][7].IMU-based motion capture systems have unique advantages, such as multidimensional information richness, compactness, and light weight.However, compared with biosensors for which signals appears earlier than in human body motion [8][9][10], human motion capture systems bear the apparent disadvantage of asynchronization with execution delay, which is a critical challenge for applications in exoskeletal robot control.
Human locomotion is a nonlinear quasi-periodic process with chaos [11] resulting from latent human musculoskeletal dynamics.During human walking, the two lower limbs operate synergistically [12][13][14][15].Both the quasi-periodicity and the synergy lead to the development of several efficient prediction methods.Firstly, by drawing evidence from the lower limb synergy, echo control has long been used in powered prosthesis for generating joint trajectories of one limb by a phase shift of the opposite limb [16].To apply echo control in prediction, the time shift needs to be correctly identified.Secondly, the complementary limb motion estimation (CLME) method [8] improves the echo control method by assuming a coupling relationship between the two lower limbs.However, its performance is not always consistent due to a significant influence from the underlying regression method [17].Thirdly, the aforementioned quasi-periodic nonlinear behaviour of walking has led to Takens' reconstruction method, which is closely related to Takens' reconstruction theorem.It is essentially a nonlinear time series dynamic factor analysis method for predicting joint trajectories [18,19].However, it requires historical system data to achieve a better performance which brings in inevitable dead time.There are also some other prediction methods, such as the autoregressive integrated model method, autoregressive integrated moving average method, autoregressive moving average model method [20], and modified recursive least square methods [21], which are also time series analysis methods but with undesired computational expense, dead time, or complex parameter tuning processes.Besides, the consistency and smoothness of human joint motion should not be ignored as human joints possess apparent biological damping [22][23][24] which damps out the high frequency vibration.Such observations indicate that the traditional Newton's method, in the short term for smoothing or correction techniques, might be beneficial.For Newton's method, estimation of the velocity and acceleration with noise rejecting ability is needed.High-gain observers could be exploited, which are usually used in robot or automation applications [25][26][27].
In this work, we consider, without loss of generality, a series of prediction methods and its improving techniques utilizing the biological features of human locomotion.Based on progressive analysis, we introduce a novel method which fuses the advantages of different approaches to achieve an even better performance verified by a daily motion dataset of a normal subject, including irregular motion.The method is also applicable to real-time embedded systems.
This paper is organized as follows.Section 2 introduces the IMU-based motion capture system.Section 3 describes different prediction methods, assessment measures, and the used dataset.Section 4 evaluates the influence of several common parameters and the prediction performance with different prediction methods by simulation, elaborates their characteristics, and then proposes our best method.Also, the real-time applicability of the proposed algorithm is demonstrated.Section 5 concludes the paper.

The IMU-Based Motion Capture System
A motion capture system using IMUs has been developed by Xeno Dynamics Co., Ltd.(Shenzhen, Guangdong, China), as shown in Figure 1.In each IMU, the gyroscope, accelerometer and magnetometer data are fused with a gradient descent algorithm with magnetic distortion compensation [28].The compact rigid container with elastic straps could be easily attached to human lower limbs.The data collected by serial communication could either be stored and processed in a standalone central control device as shown in Figure 1b or sent into the developing exoskeletal robots with a sampling frequency of 50 Hz.The test subject can move freely when wearing the motion capture system.In the experiment designed in this research, only the left hip, right hip, left knee, and right knee joint data in the human sagittal plane were collected, as illustrated in Figure 1b.

Methods
The prediction methods starts with Takens' reconstruction prediction and its enhancement with Gaussian coefficient weighting, offset correction, and a feedback scheme.Synergy of a complementary joint is exploited with two methods producing the synergy-enhanced Takens' reconstruction predictions.To further improve the prediction performance during the short term, and adaptability to irregularities, the classic Newton's method with a high-gain observer is brought in and some new methods appear.

Takens' Reconstruction Prediction and Modifications Based on Motion Smoothness
Takens' reconstruction method is correct for predicting the joint trajectories of human lower limbs, which present quasi-periodicity and nonlinearity.According to Takens' reconstruction theorem, it is possible to reconstruct the dynamics of a nonlinear system based on its past measurements.For this purpose, a current embedding vector is defined: where y(t) stands for the current state which might be position, velocity, acceleration, etc., t is the current sampling instant, p is the embedding vector length, and the sampling constant is ∆t.Takens' original reconstruction algorithm proposed by Herrmann [19] is illustrated in Figure 2:  At each sampling instant t with t ≥ l∆t, an embedding vector D(t) is acquired.The current embedding vector D(t) is then compared to all collected historical embeddings D(t i ) with l∆t ≤ t i < t by determining the Euclidean distances In Equation (2) i is the index of the ith historical embedding D(t i ) and 1 ≤ i ≤ l. l is also called the history data length.As D(t i ) is a vector composed of a series of states, the first component of D(t i ) is labeled as y(t i ) and the components are defined as y(t i − ∆t x−1 ) where ∆t x−1 = (x − 1)∆t and 1 ≤ x ≤ p.
The M best matching embeddings D( t1 ), D( t2 ), . . ., D( tM ) from the sampling instants tj with 1 ≤ j ≤ M are used to calculate a k-step prediction ŷ(t + ∆t k |t): with the weightings and In Equation ( 3) ∆t k = k∆t.When k = 1, the subscript is dropped as ∆t.k is also called prediction horizon.Although Takens' reconstruction method bears many practical advantages such as easy parameter calibration, result smoothness, accuracy and so on, its improvement is still quite possible.The natural smoothness of motion brought by the human body indicates whether high frequency jumps and sharp spikes in trajectories were initially predicted, as correction and smoothing techniques would be used to produce a more realistic result.
First, we assumed that the y(t i − ∆t x−1 ) in D(t i ) closer to y(t i ) had more of an influence on Euclidean distance determination.Thus, the Euclidean distances was calculated with the addition of Gaussian coefficients G(x): and where Q is constant and p is the embedding vector length.As the prediction result was not very sensitive to the value of Q, it was chosen to be 4 in this research, which was justified by several items of gait data.Second, the offset between y( tj ) in a best matching embedding of D( tj ) and y(t) in the current embedding vector was decided during each vector match by where y( tj ) means the past observation of the current state y(t) in the jth best matching vector D( tj ).
The obtained ∆y j is used for correction during the weighting procedure with the equation Third, the feedback scheme is integrated and tested as illustrated in Figure 3.At time t, the 1 step prediction ŷ(t|t − ∆t) of the current time is achieved at t − ∆t.The difference between the current position y(t) and ŷ(t|t − ∆t), ∆ f bk, is added to the current k step prediction ŷ(t + ∆t k |t) with a factor of k − 1.At the same time a new 1 step prediction is completed as ŷ(t + ∆t|t), which will be used in the feedback determination in the next time step.

Takens' Reconstruction Method's Synergy Improvement
The synergy feature of the lower limb joints could be used in a simple way for Takens' reconstruction-based methods, making use of the complementary joint motion data of the other limb.Basically, considering the complementary joint data as another source of history or taking it as an indirect measurement of the original joint motion should be able to provide help.This results in two distinct methods: 1.While searching for the best matching vectors, the embedding vectors are generated with the united source of the two joints' history data.The detailed process is shown in Figure 4.The two means represented different directions.They have limited theoretical foundation as the synergy mechanism is still under research and could not be commented by comparative analysis without being applied to assessment on real data.

Takens' Reconstruction Adaptability Enhancement with Newton's Method
Motion smoothness also means the assumption of local linearity is appropriate, thus the traditional Newton's method in the short term might be beneficial.A high-gain observer [25][26][27] is formulated to obtain the required velocity and acceleration for a second-order Newton's method.
where ŷ and ŷ are computed with ŷ = s and ŷ = s where α 1 = 100, α 2 = 80 and = sampling constant which was 0.02 s in this context.As a high-gain observer could be taken as a non-linear filter, α 1 and α 2 were chosen so that the amplitude of the data after the filter did not change and the phase delay in time domain was constant.It usually brought in a delay of about one sampling cycle.Equation (10) roughly compensates the delay by taking a one step longer prediction at k + 1 as the k step result.Finally, the Newton-based and Takens' reconstruction-based methods could be fused to achieve some methods with unique adaptability enhancement.

Methods Summary
Summarizing the above methods, besides the fundamental Takens' reconstruction and Newton's method, there are 12 new prediction methods that are introduced and analysed which are listed in Table 1.The relationship between different prediction methods with different improving approaches is illustrated in Figure 6.

Assessment Data
The dataset used to assess the performance of different prediction methods is 208.5 s long with a sampling rate of 50 Hz, including the joint angles of the knees and hips of a test subject in a sagittal plane.It is comprised of the following continuous motion process segmented and highlighted in Figure 7: 5 stand-sits, 5 steps walking, 2 left turns, 5 steps walking, 2 left turns, 5 steps walking, stepping onto a treadmill, 1 km/h, 2 km/h and 3 km/h walking on a treadmill, stopping and leaving the treadmill, 5 steps walking and sitting down.In order to assess the adaptability of different methods the following segments were defined as irregularities: stand-walk transition, turns, stepping on a treadmill, stopping and leaving the treadmill, and stand-sit.

Assessment Measures
To quantitatively evaluate the algorithms' prediction performance, the prediction ratio (PR) representing prediction errors, and the smoothness factor (SF) measuring predicted trajectory smoothness are calculated by the following equations [19]: where ŷ(t|t − ∆t k ) represents the predicted current position ŷ(t) at t − ∆t k , which is obtained by shifting the predicted results with ∆t k and where t End is the time when prediction ends and f (t) is the filtered result of y(t) by a zero-phase filter with a cut-off frequency of 5 Hz.A higher PR value means higher accuracy, while a lower SF value means smoother motion.In order to provide an overall relative comparison among different methods, PR and SF are combined into an overall performance index (OPI), which is computed with where i represents the index of a specific method involved in the comparison.It should be emphasized that such an index is not an absolute criterion but a relative one.The preference over prediction accuracy and smoothness could be modified by the coefficient α.
In addition, the adaptability of prediction methods to irregularities happening from time to time is very essential.Such capability could be assessed by isolating the irregularities, including the beginning and termination of walking and the transition between different motion speed and modes.

Implementation and Assessment
In practical applications, prediction horizon may change in real-time, so a good algorithm does not only perform well at a fixed prediction horizon, but also over a certain range.Before the performance evaluation over different prediction horizons, the common parameters used in Takens' reconstruction-based algorithms must be reasonably selected.

Parameter Selection for Takens' Reconstruction-Based Methods
For Takens' reconstruction-based methods common parameters were evaluated, including the history data length, embedding vector length, and the number of weighted best matching vectors.A group of basic parameters are listed in Table 2.During the evaluation the evaluated parameters were taken as variables while others were fixed.Only the left hip data were shown in the following paragraphs because the four joints results exhibited the same varying patterns which meant the ruling laws were the same for different methods.
Theoretically, a longer history is beneficial and Figure 8 gave the proof with continuous escalation of PRs.A steep rise ended at the value of 360, which was the chosen value in Table 2.The history data length had limited influence on SFs after 360, because a longer history could not generate better historical embedding vectors.
The length of embedding vectors represents the orders of the reconstructed system.Excessive higher order might cause overfitting, which reduces the prediction accuracy and increases unnecessary computational expense.Gaussian weightings mitigate the problem by reducing the influence of the history data far away and the offset correction and feedback scheme also provide powerful fixation.Figure 9 proved the analysis.For the SF values, the expending of embedding vectors could improve them more or less for all methods.Thus, considering that a length of 20 is not an burden to the embedded system and it has limited influence on most of the improved methods, the embedded vector length was chosen to be 20.As shown in Figure 10, when the weighted best matching vector number increased, the PRs of different methods rose to the highest at 5 and the trend then descended afterwards.Meanwhile, for the FSs after 5, they also entered an slow falling state.This was reasonable because the matched vector quality degraded when the number increased to a certain level, and their weighting values fell quickly.At the same time, the increasing number of weighted history embedding vectors worked as a smoother by averaging down the oscillations.Surprisingly, most newly-introduced Takens' reconstruction-based methods exhibited robust performance with parameter variation.Thus, the parameter robustness of the newly introduced methods was discovered.
Summarizing the previous analysis, the basic parameters are selected with enough confidence, and could be applied to algorithm performance assessment.

Results and Assessment
The prediction methods listed in Table 1 and evaluated (Figure 7) demonstrated the dataset with respect to an extending prediction horizon from 1 to 20, which was the requirement of exoskeletal robot sensing for control.For OPI calculation, α in Equation ( 16) was chosen as 0.8, which emphasized accuracy over smoothness.The evaluations were performed with the simulation of real-time operation by pushing the data one by one in a time sequence into the last position of a data pool, which was fair for all the methods.
Studying the results by comparison the pros and cons of different methods and improving techniques were clear: 1.It can seen more clearly in Figure 11 that when compared with the original Takens' reconstruction method, the PR and SF values of Gw and GwOC had significant improvements, especially during the lower prediction horizon range, which meant the Gaussian weighting and offset correction were effective and the assumption of smoothness and consistency of human motion was reasonable.2. Also, it can be inferred by comparing GwOc, GwOc-KC, and GwOc-US that for prediction accuracy the synergy improvement of united source of history embedding vector generation had a positive effect on smoothness improvement and a negative effect on accuracy.In contrast, the improvement of Kalman fusion with complementary joint prediction played an opposite role.Hence, the two synergy improvements would be chosen based on application preference.3. Figures 11-13 show that feedback scheme would push up the PR value and hold back the FS value.The FS curves of the methods with feedback scheme deviated a lot from the other methods with the extension of the prediction horizon.This was always true by comparative inspection.
Hence, if prediction accuracy is expected with highest priority, the feedback scheme should be considered.

Figure 11.
Comparison of prediction performance between the Takens', Gw, GwOc, GwOcFb, GwOc-KC, and GwOc-US methods (left hip joint data in the sagittal plane).OPI: overall performance index.4. In Figure 14, initially in the prediction horizon range of 1-3, the Newton's method with a high-gain observer had a medium accuracy and also adequate smoothness, but after 3, the Newton's method-based predictions degraded the fastest in both accuracy and smoothness.Meanwhile, for Takens' reconstruction-based methods, the curves changed moderately.The NHgo-fused methods, GwOC-US-NHgo and GwOC-US-NHgo, exhibited a slight drop in PR values compared to the methods without fusion, but the smoothness was more elevated.As a result, their overall performance was increased.5.By observing Figure 15 for short-term prediction such as in steps 1-3, the basic Takens' reconstruction method's performance was unsatisfactory, but its OPI value rose gradually, which meant its relative performance over other methods increased.The OPI curve with the KC-NHgo method changed gradually from an descending trend to ascending, which meant it was superior to other methods with respect to an extended prediction horizon.Considering the OPI values with designed weighting, the fully improved method, GwOc-KC-NHgo, was apparently better across the prediction horizon range.The OPI values are relative criterion and are significantly affected by the weighting value.If smoothness is preferred, the GwOc-KC-NHgo method is superior.The adaptability to irregularities could be assessed visually and quantitatively by local irregularity tracking performance observations.The segments considered as irregularities in the assessment dataset were isolated and the predictions with different methods based on the parameters in Table 2 were executed.Some typical local trajectories are shown in Figure 16.
It can be seen in Figure 16a that a flat curve with zero values from 0 to 0.4 s exists at the very beginning of the predicted curves of Takens' reconstruction methods.Noting that there are no human data (zero), this meant dead time existed.Zero dead time was gained by other methods via offset correction, feedback, and fusion with the Newton-based methods.
Figure 16b presents the performance of different methods during transitions.Except for the basic Takens' reconstruction method, other methods were able to provide more or less prediction.The NHgo method-predicted curve fluctuated severely.Although the methods fusing the NHgo method, GwOc-KC-NHgo and GwOcFb-KC-NHgo, were influenced by these fluctuations locally, the fusion provided a better inhibiting effect.As the NHgo method exhibited outstanding prediction capability during periodicity loss, the methods fused with it acquired better prediction ability than without it.Figure 16c,d were segments where periodicity was invisible, thus it was quite obvious that Takens' reconstruction method could hardly provide prediction, while other methods still worked more or less.
Quantitatively, the numerical criteria could gave a more objective judgement.Figure 17 compared the PR, OPI, and SF values between the whole dataset and the dataset during irregularities only.There are three phenomena worth noting in the figure : 1.The closing up of the assessment measure value gaps between the whole dataset and dataset during irregularities 2. The improvement in the values' changing trends.3. The signs of differences in assessment measure values between the whole dataset and dataset during irregularities.
Compared with the basic Takens' method, the GwOc method closed up the gaps of PR, SF, and OPI values between the whole dataset and the dataset during irregularities.Also, all the values were improved with a positive trend.Compared with the GwOc method, the GwOcFb method also made the gap even smaller, but the SF values rose, which brought down the OPI values.With a similar analysis, the GwOc-KC and GwOc-US methods elevated the smoothness with prediction accuracy degradation, and also narrowed the gaps, which means the synergy improvement was effective in improving the prediction performance.The PR value difference using the NHgo method with respect to the whole dataset and the dataset of irregularities is negative, which means the prediction accuracy adapted better to irregularities, while for the SF values the difference is positive, which means the smoothness of predicted trajectories also adapted to irregularities better than periodicity.In contrast, the Takens', GwOc, GwOc-KC, and GwOc-US methods adapted to the whole dataset better.
Thus, by fusing the two types of prediction methods, the GwOcFb-KC-NHgo and GwOc-KC-NHgo methods achieved acceptable accuracy and smoothness, and their adaptability to irregularities was improved.The OPI values for irregularities also demonstrated better adaptability of GwOcFb-KC-NHgo and GwOc-KC-NHgo over other methods, performing extremely well with the whole dataset and the irregular segments.
Based on the previous analysis, the GwOc-KC-NHgo method was the first choice for prediction because of its high prediction accuracy, smoothness, and adaptability to irregularities with no dead time.

Real-Time Demonstration
The GwOc-KC-NHgo methods were implemented on the human lower limb motion capture system in Figure 1 with the basic parameters designed as in Table 2.The algorithm's computational time was between 3.2 and 3.3 ms (Figure 18) with tiny fluctuations which could satisfy the real-time execution requirements.The computational expenses of the other methods were similar, but could not be shown here.

Conclusions
In exoskeletal robot control, it would be much easier to handle the human robot coordination problem if the human motion was known in advance.Prediction is a helpful approach.Based on the biological features of human walking, a series of bio-inspired prediction methods were proposed and implemented.The parameters of Takens' reconstruction-based methods were analysed and selected with the parameter robustness of newly introduced methods demonstrated as well.The comparative assessment demonstrated the offset correction and that Gaussian weighting techniques could effectively improve the prediction performance of Takens' reconstruction method in both prediction accuracy and trajectory smoothness.The feedback scheme was useful in prediction accuracy improvement but degraded the smoothness.The synergy improvement of the coupled joints' united source of historical embedding generation exhibited better performance in prediction accuracy than Kalman fusion with the complementary joint prediction and poorer performance in trajectory smoothness.The second-order Newton's method with a high-gain observer gave acceptable performance with short prediction horizons.When considering improving the adaptability to irregularities, the Takens' reconstruction-based methods were fused with Newton-based methods.On comparison with their contending methods, the new hybrid methods showed mixed features with nearly zero dead time, a slight deterioration in prediction accuracy and smoothness in short prediction horizons, and improvement in longer prediction horizons.The methods should be tailored for the detailed requirements of exoskeletal robot applications.When accuracy is emphasized, the GwOc-KC-NHgo method is recommended.The real-time computational expense would not be a burden for the methods introduced.

Figure 1 .
Figure 1.The inertial measurement unit (IMU)-based lower limb motion capture system.(a) One IMU module.(b) The system worn on a testing subject and the illustration of joint angles in the human sagittal plane collected in this research.

Figure 6 .
Figure 6.The relationship map of different methods and improving approaches.

Figure 7 .
Figure 7. Dataset for prediction method assessment (the joint angles of the knees and hips of a test subject in a sagittal plane collected by the inertial measurement unit-based motion capture system).

3 Figure 8 .
Figure 8. Influence of history data length used in embedding vector generation (left hip joint data in the sagittal plane).SF: smoothness factor; PR: prediction ratio.

Figure 9 .
Figure 9. Influence of embedding vector length (left hip joint data in the sagittal plane).

3 Figure 10 .
Figure 10.Influence of weighted historical embedding vector number (left hip joint data in the sagittal plane).

Figure 12 .
Figure 12.Comparison of the prediction performance of the GwOc-KC, GwOc-US, GwOcFb-KC, and GwOcFb-US methods (left hip joint data in sagittal plane).

Figure 16 .
Typical adaptability of different methods to irregularities (left hip joint angle in sagittal plane).(a) Dead time during initiation; (b) Speed transition from 1 km/h to 2 km/h; (c) Stand-walk transition; (d) Non-periodic motion during "stopping and leaving the treadmill".

Figure 17 .
Figure 17.Performance of different methods during irregularities (left hip data).
∆t k |t) with the Kalman filter and obtain ŷf use (t + ∆t k |t).The detailed process is shown in Figure5.
2. Predict k step position ahead as usual and get ŷ(t + ∆t k |t).Then take the complementary joint as the source of embedding vector generation and predict ŷc (t + ∆t k |t).Finally, fuse ŷ(t + ∆t k |t) and ŷc (t +

Table 1 .
Prediction methods in this article and their abbreviations.

Table 2 .
Basic parameters for parameter evaluation.