Probabilistic Sensitivity Amplification Control for Lower Extremity Exoskeleton

Abstract

[...]the Gaussian process (GP) state-space model presented in [19] has not addressed the crucial issue of computation expense, which we will carefully discuss in this section. [...]the graphic explanation of the crucial difference between the GP state-space inference [17,18,19] and our proposed distributed hidden-state identification are shown in Figure 2. [...]we apply a more efficient method in this paper, i.e., expectation maximization with Monte Carlo sampling. [...]the general class of latent state-space models with the Gaussian process is given below: xk=h(xk−1,ak−1)+ωk,ωk∼N(0,Q), zk=g(xk)+υk,υk∼N(0,R), With k=1,...,T . The x∈RD is a latent state that evolves over time, while z∈RD can be read from actual measurement data. In light of the work in [22], the updated posterior distribution of the latent state for the m-th (m = 1,⋯,M ) GP can be obtained: (μk|kx)m=(μk|k−1x)m+(Ck|k−1xz)m ((Ck|k−1z)m)−1(zk−(μk|k−1z)m), (Ck|kx)m=(Ck|k−1x)m−(Ck|k−1xz)m ((Ck|k−1z)m)−1 ((Ck|k−1zx)m)T, (μk−1|Tx)m=(μk−1|k−1x)m+(Jk−1)m((μk|T)m−(μk|kx)m), (Ck−1|Tx)m=(Ck−1|k−1x)m+(Jk−1)m((Ck|Tx)m−(Ck|kx)m)((Jk−1)m)T. Not akin to [22], to infer a hidden state with a corresponding observation input, we should combine all the M GPs to fusion of the M predictions. [...]in this paper, we compare four fusion algorithms, i.e., product of expert (PoE) [23], generalized product of expert (gPoE)[24], Bayesian committee machine (BCM) [25] and robust Bayesian committee machine (rBCM) [20]. The gPoE distributed scheme provides more reasonable identification of the hidden states than any other models in Figure 4. [...]Figure 5 shows two aspects of the distributed identification framework with the different number of GP experts: (1) a comparison of approximation quality; (2) computation expense.

1. Introduction

In recent decades, with the dramatic progress in computing, control and sensing, significant types of prostheses, orthoses and exoskeletons have been designed and developed to assist and support human limbs for various tasks. In order to increase the operator’s strength and endurance, the exoskeleton is a typical class of wearable mechatronic robotic systems, which are designed to meet this demand. A particular feature of this kind of exoskeleton is that they are fully driven by the pilot intention and follow the pilot movements in synchrony. The application feasibility includes functional rehabilitation for the injured, physical assistance for the disabled and elderly and fatigue relief for heavy-duty workers. To recognize human intention, two main groups of sensors, i.e., cognitive-based sensors, as well as physical-based sensors, are usually applied to capture the information of the human interaction or the muscle signals.

By measuring the electric signals from the musculoskeletal system or the nervous system, the operator’s intention read by the cognitive-based sensors can be interpreted as the input to the exoskeleton system controller. A successful practical application of these kinds of exoskeletons is the hybrid assistive limb (Tsukuba University and Cyberdyne, Japan) [1,2], which is structured for whole-body assistance. The muscular activity signal is collected using surface electromyography (sEMG) sensors. In addition, the exploration for brain-machine interface (BMI) control is implemented at the University of Tubingen in Germany [3]. The hybrid BMI system consists of an EMG continuous decoder, as well as an electroencephalogram (EEG) binary classifier. Consequently, the crucial advantage of cognitive-based sensors is that the pilot or the patient intention can still be read, even if the subject cannot provide necessary joint torques. Thus, these types of the sensor systems are widely utilized in rehabilitation and medical application.

An alternative sensor system recognizes the human intentions by applying physical-based sensors. The Body Extender [4] is designed to handle and transport heavy loads up to 100 kg. Consequently, the human intentions are read by the force sensors mounted on the interaction part between the pilot and exoskeleton. The velocity input to the controller is converted by the impedance control algorithm [5,6]. Moreover, for low impedance, back drivability and precise control, a compact rotary series elastic actuator is designed in [7] to fulfill these requirements. The input to the system controller is derived from the deformation of the series elastic module. Also note that an improvement of system stability is introduced in [8], which addresses the time-delay effect of the series elastic module.

Apart from the sensor-based human intention recognition, the sensitivity amplification control is proposed in [9,10] for the Berkeley Lower Extremity Exoskeleton. The design conception of control algorithm aims to enable a pilot’s locomotion with little effort on any type of terrain. In addition, the control scheme demands on interaction sensors between the pilot, as well as the exoskeleton. Nevertheless, the controller has little robustness to model uncertainties or parameter variations.

Although the sensors mounted between the pilot and the exoskeleton can simplify the control algorithm, a crucial issue that should be addressed is that the time-delay effect cannot be avoided. Consequently, the control performance is limited due to this issue, even if the robustness can be enhanced by the compensation of the disturbance observer [11] or the internal model control [12]. The problem of time-delay effect may be solved by the cognitive-based sensors [13,14], since the human intention can be acquired before the movement primitives. Nevertheless, these signals derived from the sensors suffer from unwanted noise, as well as the overlap of the spectrum with other signals.

The sensitivity amplification control (SAC) control scheme has provided a solution for the human performance augmentation exoskeletons. However, such an algorithm suffers from weak robustness, and it barely destroys the control system. Thus, a sophisticated identification process is required to improve system stability, as detailed in [15]. The second drawback is that the fixed sensitivity factor cannot manage the uncertainties of the interaction model.

In this paper, we retain the advantages of the SAC control scheme and propose a novel control scheme, so-called probabilistic sensitivity amplification control, to address the two crucial issues, as shown in Figure 1. Rather than identifying the system dynamics, we infer the state-space model based on a novel distributed identification process. To address the challenge of dealing with the uncertainties of the interaction model, an online reinforcement learning algorithm based on the derandomized evolution strategy with covariance matrix adaptation (CMA-ES) is presented for the adaptation of sensitivity factors. Thus, the sensitivity factors can be effectively updated in a constraint domain of interest according to the torque tracking error penalty.

The remainder of the paper is organized as follows:

After the Introduction, Section 2 is devoted to outlining the main routine of our proposed distributed hidden-state identification to address the issue of weak robustness. In Section 3, the novel evolving learning of sensitivity factors is presented to deal with the variational HRI. The implementation of the experiments and the comparison of the control algorithms are presented in Section 4. Finally, Section 5 concludes with a summary.

2. Distributed Hidden-State Identification

For realizing efficient control of the exoskeleton, it is essential to understand the human intention. According to our control scheme, the intention can be interpreted as the system state variables. Although tremendous attention has been paid to filtering and smoothing [16,17,18], these algorithms do not focus on the identification process, but on inference. In addition, the Gaussian process (GP) state-space model presented in [19] has not addressed the crucial issue of computation expense, which we will carefully discuss in this section. Therefore, the graphic explanation of the crucial difference between the GP state-space inference [17,18,19] and our proposed distributed hidden-state identification are shown in Figure 2.

In this section, firstly, we give a comprehensive explanation of the distributed model learning for the identification process based on the expectation maximization. Then, we introduce the online inference under a distributed framework and Markov chain Monte Carlo, also along with the fusion method, so-called generalized product-of-expert (PoE).

2.1. Model Learning

Although the Gaussian process has a substantial impact in many fields, the crucial limit is that the computation expense scales in O(^N3) for training, while regarding prediction, the computed expense is O(^N2) if caching the kernel matrix. Thus, we employ the distributed Gaussian process [20] to address this issue, since this method can handle arbitrarily large datasets with sufficient hierarchical product-of-expert (PoE) models.

The Gaussian process, fully specified by a mean function _μf and a covariance function _kf , is a probabilistic distribution over functions f. The latent function _yi=f(_xi)+ω,ω∼N(0,_σω2) can be learned from a training dataset {X,z} , where X=[_x1,...,_xn] and z=[_z1,...,_zn] . In addition, the posterior predictive distribution of the function _f*=f(_x*), by a given set of hyper-parameters θ is defined with mean and covariance, respectively

E[_f*]=m(_x*)=_k*T ^(K+σ2I)−1z=_k*β,

var[_f*]=^σ2(_x*)=_k**−_k*T ^(K+σ2I)−1 _k*,

With _x* defining input, _k*=k(X,_x*) and _k**=k(_x*,_x*) . Note that the hyper-parameters θ of the covariance function consist of the length-scale l , the signal variance _σf2 and the noise variance _σω2 . Besides, the hyper-parameters are of significance to the covariance function defining the non-parametric model behavior [21].

According to the independent distributed assumption, the marginal likelihood p(y|X,θ) can be decomposed into M individual products:

logp(y|X,θ)≈∑m=1Mlogp(^y(m)|^X(m),θ),

With M the total number of experts. In the model learning process, the hyper-parameters θ are shared for all M GPs. Moreover, a usual method for optimizing the hyper-parameters is to maximize the evidence maximization. Nevertheless, we apply a more efficient method in this paper, i.e., expectation maximization with Monte Carlo sampling.

Therefore, the general class of latent state-space models with the Gaussian process is given below:

_xk=h(_xk−1,_ak−1)+_ωk,_ωk∼N(0,Q),

_zk=g(_xk)+_υk,_υk∼N(0,R),

With k=1,...,T . The x∈^RD is a latent state that evolves over time, while z∈^RD can be read from actual measurement data. Both additive system noise, as well as measurement noise are assumed i.i.d. If not stated otherwise, both the measurement function and dynamic function are implicit, but instead can be described with the probability distributions over these Gaussian functions. Here, the functions are a Gaussian process non-parametric model, defined as h∼_GPh and g∼_GPg , separately.

Limited by the space of this paper, we briefly explain the learning process of the dynamic function h∼_GPh . In the expectation step, we compute the posterior distribution of the system states q(X)=p(^X(m)|^y(m),θ) , where ^X(m)=(_xk−1,_ak−1) and ^y(m)=_xk . Nevertheless, the posterior distribution is not a Gaussian process and hence is not easy to obtain analytically. Consider the situation in which we draw samples from the posterior distribution and the marginal likelihood is approximated by Monte Carlo integration, as shown in Figure 3. Consequently, in the maximization step, the hyper-parameters are given by maximizing the obtained likelihood _Eq[p(^X(m)|^y(m),θ] with stochastic conjugate gradient optimization. The whole learning routine is give in Algorithm 1. Note that though the learning procedure is time consuming, this is nota problem since the training process of hyper-parameters is off-line.

Algorithm 1: Calculate hyper-parameter values with Monte Carlo expectation maximization.

[Image omitted. See PDF.]

where p(_xk|z+1;k) is defined as a prior of the posterior distribution and p(_zk|_xk) is determined by the measurements. In addition, to improve the identification process, we reuse the observation data, i.e., current, previous and future:

γ(_xk)=p(_xk|_z1:T).

For the recursive process, the smooth hidden states γ(_xk−1) are given from k=T to 2:

γ(_xk−1)=∫p(_xk−1|_xk,_z1:k−1)γ(_xk)d_xk.

In light of the work in [22], the updated posterior distribution of the latent state for the m-th (m = 1,⋯,M ) GP can be obtained:

^(_μk|kx)m=^{(_μk|k−1x)m}+^{(_Ck|k−1xz)m ((_Ck|k−1z)m)−1}(_zk−^{(_μk|k−1z)m}),

^(_Ck|kx)m=^{(_Ck|k−1x)m}−^{(_Ck|k−1xz)m ((_Ck|k−1z)m)−1 ((_Ck|k−1zx)m)T},

^{(_μk−1|Tx)m}=^{(_{μk−1|k−1x})m}+^(_Jk−1)m(^(_μk|T)m−^(_μk|kx)m),

^{(_Ck−1|Tx)m}=^{(_{Ck−1|k−1x})m}+^(_Jk−1)m(^{(_Ck|Tx)m−(_Ck|kx)m})^{((_Jk−1)m)T}.

Not akin to [22], to infer a hidden state with a corresponding observation input, we should combine all the M GPs to fusion of the M predictions. Therefore, in this paper, we compare four fusion algorithms, i.e., product of expert (PoE) [23], generalized product of expert (gPoE)[24], Bayesian committee machine (BCM) [25] and robust Bayesian committee machine (rBCM) [20]. The only difference between the gPoE and PoE is the different precisions, since the important weights are set to _βk=1/M to avoid an unreasonable error bar. Moreover, the improvement made in the BCM is that it incorporates the GP prior explicitly when predicting the joint distribution. Hence, the prior precision should be added to the posterior predictive distribution. The rBCM seeks to add the flexibility of gPoE to BCM to enhance predictive performance.

Moreover, we conduct an experiment with a one dimension problem [26]. Figure 4 shows the performance comparison of four distributed identification process (PoE , gPoE, BCM, rBCM) and robust Rauch-Tung-Striebel smoother for GP dynamic systems (GP-RTSS) [18], which is implemented under the full GP scheme. The PoE and BCM distributed schemes overestimate the variances in the regime of the series hidden states, which leads to overconfident precision compared with GP-RTSS. Moreover, the prediction of the rBCM is too conservative, although the rBCM aims to provide a more robust solution. The gPoE distributed scheme provides more reasonable identification of the hidden states than any other models in Figure 4.

Moreover, Figure 5 shows two aspects of the distributed identification framework with the different number of GP experts: (1) a comparison of approximation quality; (2) computation expense. According to the information in Figure 5, for a certain set of training data, with the increase of the GP experts, the computation time significantly decreases. However, the computation expense reaches a minimum when there are five GP experts. Besides, the interesting properties are that there is no significant difference (root-mean-square-error (RSME): 0.01 and negative log-likelihood (NLL): 0.2) between the identification framework with one GP expert (full GP) and more than one expert. Furthermore, note that a better approximation quality can be expected with ten GP experts with the minima of RSME or with eight GP experts according to the NLL criterion. Thus, although we do not argue that our proposed hidden-state identification scheme outperforms the full GP identification algorithm concerning precision, the computation expense has been dramatically reduced. Moreover, the overall identification procedure is concluded in Algorithm 2.

Algorithm 2: Inference with Monte Carlo Markov chain.

[Image omitted. See PDF.]

The interaction between the pilot and the exoskeleton is given:

d=τ−P(θ)−(1−^α−1)[M(θ)θ¨−C(θ,θ˙)θ˙],

With perfect estimation M^=M,C^=C and P^=P and sufficient large α . The interaction d can approach zero. Nevertheless, the sensitivity factor should not be selected with a sufficiently large parameter due to the weak robustness. Moreover, in an actual situation, the interaction varies owing to stochastic environment factors. Therefore, choosing a fixed sensitivity factor is not a good option. An improvement has been made in [27], where the sensitivity factor is learned with Qlearning. Compared with the classic SAC, the interaction force has been reduced. However, the learning algorithm mainly copes with the discrete region of interest. In addition, the optimization of the constraint domain of the sensitivity factor is not addressed in the reference.

In this paper, the sensitivity factor is evolved in real time in a continuous constraint domain. By applying (1+1) covariance matrix adaptation evolution strategy ( (1+1) -CMA-ES) [28], the sensitivity factor seeks the optimal value owing to tracking error penalty.

3.1. (1 + 1)-CMA-ES

Derived from the CMA-ES, a crucial feature of the (1+1) -CMA-ES is that it directly operates on Cholesky decomposition and is about 1.5 -times faster than the classic CMA-ES. The (1+1) -CMA-ES adapts not only the entire covariance matrix of its offspring candidates, but also the entire step size. Since the covariance matrix C is definitely positive, owing to Cholesky decomposition, C=^AAT , with A an n×n matrix. Thus, the offspring candidate y is written:

y=x+σAz,

where σ is the step size, x is the corresponding parental candidate and z∈^Rn is derived from normal distribution. In addition, the adaptation strategy maintains a search path in order to save an exponentially-fading record s∈^Rn . Thus, a successful offspring candidate can be obtained from the covariance matrix:

C←(1−_ccov+)C+_ccov+ ^ssT,

With _ccov+∈(0,1) . Rather than dealing with the Cholesky decomposition operation in every loop, the matrix A can be updated as follows:

A←1−_cco^v+A+1−_ccov+^||w||2(1+_ccov+ ^||w||21−_ccov+−1)s^wT,

With w=^A−1s . To enhance the learning efficiency, if the offspring candidate solution is extremely unsuccessful, the algorithm should also learn from failure cases. Therefore, the covariance matrix update is written:

C←(1+_ccov−)C−_ccov−(Az)^(Az)T,

where _ccov−∈(0,1) ; similarly, the equivalent update of matrix A is given:

A←1−_cco^v−A+1−_ccov−^||w||2(1−_ccov+ ^||w||21−_ccov−−1)s^wT.

The graphic interpretation of CMA-ES is shown in Figure 6. More specifically, the contributions of each sample are ranked according to corresponding penalties, which will be explained in the next subsection. Not all samples, but several high-contribution points are used during the update procedure. Thus, the new distribution is evolved from the previous step distribution with an adaptive step size.

3.2. Constraint Optimization

Since the sensitivity factor ^α−1∈(0,1) is constrained regarding the optimization domain, to deal with this problem, the variance of the offspring distribution is reduced in the normal vector direction in the vicinity of the corresponding parental candidate. For the constraints _vj,j=1,⋯,m , an exponential fading vector is defined following [29] and updated when the optimization is lessened by the boundary:

_vj←(1−_cc)_vj+_ccAz,

where _cc∈(0,1) is a parameter that determines the convergence speed. Equation (20) can be seen as a low-pass filter, which guarantees that the tangential components of Az have been canceled out.

In these situations, the Cholesky factor A should also be rewritten to meet the demand of the restrained boundary:

A←A−β_{∑j=1m Γgj(y)>0}∑j=1m_{Γgj(y)>0vi wjTwjT wj},

With _wj=^A−1 _vj and _Γgj(y)=1 if _gj(y)>0 , otherwise zero. The update size is controlled by the parameter β . Moreover, compared with Equations (17) and (19), Cholesky decomposition A without a scalar greater than one aims to decrease the update step size to approach the boundary of constraints.

The solution to generating offspring candidates is explained in Figure 7. In the direction of the normal vector, the variance of the offspring is reduced as shown in Figure 7a. After the operation, the evolved variance is tangential to the constraint boundary in Figure 7b. Thus, the offspring candidates can be sampled and derived from the evolved variance A . Note that, in Figure 7c, after constraint optimization, the constraint angle has been increased due to the change of the offspring distribution.

In terms of the sequence hidden-states (_q¨ides,_q˙ides,_qides) , as well as actual data, the penalty function for updating of the sensitivity factor is given:

J(τ)=∑i=1i=k|_q¨ides|+^{(_q˙ides−_q˙i)2}+^{(_qides−_qi)2},

where _q˙i and _qi are the joint velocity and angle, separately. Note that not akin to the control algorithm in [27], where the penalty function consists of the previous loop trajectory data, the desired control input (_q¨ides,_q˙ides,_qides) is computed and predicted from the identification process. Therefore, our scheme is more suitable for improving real-time control performance. The first term in Equation (22) is to avoid high jerk. Furthermore, the last two terms are penalized for tracking error. From our intuition, with less interaction force between the human and the exoskeleton, the pilot will feel more comfortable. According to this point view, punishing the last two terms may enhance the tracking performance. The main routine of the evolving learning of sensitivity factors scheme is presented in Algorithm 3.

Algorithm 3: Evolving learning of the sensitivity factor.

[Image omitted. See PDF.]

Figure 1.Proposed control framework. In the feedforward control loop, the sensitivity amplification control (SAC) controller with the evolving sensitivity factor (in black dashed arrow) provides the necessary driven torque according to the variational HRI (human-robot interaction, given in double blue arrows). Moreover, to improve the control performance, the hidden-states are identified in the feedback control loop according to the observation information. Note that, the red color arrows stand for input and output, while the black arrows present the signal direction.

[Image omitted. See PDF.]

Figure 2.Graphic explanation of the Gaussian process (GP) state-space inference and distributed hidden-state identification. Derived from the distributed GP framework, the hidden state variables (embedded in the gray circles) are identified by the dynamic models_hi()and the measurement models_gi(), separately. Also, the connection is presented in the red arrows and the influence of the input to the control plant is shown in the blue arrows. Note that the hyper-parameters_ϕh,_σhrelated to the dynamic models are shared in the training process (in green arrows), along with the hyper-parameters_ϕg,_σg for measurement models. (a) GP state-space inference; (b) proposed distributed hidden-state identification.

[Image omitted. See PDF.]

Figure 3.The graphic interpretation of the Monte Carlo sampling. The real distribution in the dashed blue line is not the GP type, but approximated as GP in the solid blue line with the mean and covariance of the corresponding data.

[Image omitted. See PDF.]

Figure 4.Comparison of four data fusion methods for distributed hidden-state identification. In each subfigure, the confidence interval of the comparison algorithm is filled between two green lines. Furthermore, the GP-RTSS (full GP) to be approximated is shown in the blue shaded area, representing95% of the confidence interval, and the ground truth is given in the red dashed line. (a) PoE; (b) gPoE ; (c) BCM; (d) rBCM.

[Image omitted. See PDF.]

Figure 5.Computation time, RSME and NLL for the different number of GP experts. The computation time is measured in MATLAB 2015b (MathWorks Company, USA) on a laptop ( 2.6 GHz cores and 8 GB RAM) with a Ubuntu 16.04 system (Canonical Ltd., UK). Each case of different numbers of GP experts is tested 20 times for the purpose of reducing the stochastic data influence.

[Image omitted. See PDF.]

Figure 6.Visualization of CMA-ES exploration and updating rules. In (a), the new distribution is shown in the red ellipse, while the last-step distribution is in blue. The quantity of the contributions of each sample is defined with the area of the green circles. Consequently, these contributions are also ranked in (b). (a) Sampling, estimation and new distribution; (b) contributions of the samples.

[Image omitted. See PDF.]

Figure 7.Evolved variance of the offspring distribution along the normal direction of the boundary. (a) Original distribution; (b) reduced distribution; (c) inverse Cholesky decomposition. Additionally, the parent is donated as red point, while the offspring distribution is given in green circle (ecllipse). The gradient direction is define with blue dashed line.

[Image omitted. See PDF.]

Figure 8.Exoskeleton robotic platform. The safety concern for our pilot is from two aspects, i.e., hardware precaution, as well as programming. Accordingly, exceeding of the motion is prohibited by the mechanical limit, and the system will be shut down if an error occurs.

[Image omitted. See PDF.]

Figure 9.Evaluation of m GPs (m=10 ). According to the distributed framework, the ten GPs share the same hyper-parameters. (a) GP1; (b) GP2; (c) GP3; (d) GP4; (e) GP5; (f) GP6; (g) GP7; (h) GP8; (i) GP9; (j) GP10.

[Image omitted. See PDF.]

Figure 10.Cross-validation of the learned identification model. For the purpose of verifying the learning result of the identification model, we conduct a cross-validation with 600 successive data pairs. The state-space inferences of the hip joint, ground truth and the prediction error are presented in the above figure.

[Image omitted. See PDF.]

Figure 11.Flat walking experiment. (a) Initial ground contact; (b) heel rise; (c) initial swing; (d) middle swing; (e) terminal swing; (f) heel strike; (g) full ground contact. The interaction between the pilot and the exoskeleton is made possible with the compliant cuff. Moreover, the pilot is asked to naturally sway his left leg.

[Image omitted. See PDF.]

Figure 12.Tracking performances and evolution of sensitivity factors. (a) Tracking performance of the hip joint; (b) tracking performance of the knee joint; (c) evolution sensitivity factor of the hip joint; (d) evolution sensitivity factor of the knee joint.

[Image omitted. See PDF.]

Figure 13.The HRI comparison of four control algorithms. (a) SAC; (b) SAC + Qlearning; (c) SAC + CMA-ES; (d) proposed algorithm. Note that the HRI force is estimated by the disturbance observer since there is no sensor mounted between the human and exoskeleton.

[Image omitted. See PDF.]

Table 1.Average performance of HRI (RMSE).

RMSE SAC [10] SAC + Qlearning [27] SAC + CMA-ES Proposed Algorithm

Hip 9.493±1.444 8.245±1.301 7.822±1.522 3.600±0.639

Knee 1.598±0.335 1.253±0.281 1.039±0.265 0.582±0.152

Acknowledgments

Part of this work has received funding from National Science Foundation of China under Grant No. 51521003. We gratefully acknowledge the constructive comments and suggestions of the reviewers.

Author Contributions

Likun Wang and Yi Shen conceived of and designed the experiments. Likun Wang, Wei Dong and Zhijiang Du performed the experiments. Likun Wang analyzed the data; Guangyu Zhao contributed experiment platform and devices. Likun Wang and Wei Dong wrote the paper.

Conflicts of Interest

The authors declare no conflict of interest.

© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

References

1. Kawamoto, H.; Sankai, Y. Power assist method based on phase sequence and muscle force condition for HAL. Adv. Robot. 2005, 19, 717-734.

2. Cruciger, O.; Schildhauer, T.A.; Meindl, R.C.; Tegenthoff, M.; Schwenkreis, P.; Citak, M.; Aach, M. Impact of locomotion training with a neurologic controlled hybrid assistive limb (HAL) exoskeleton on neuropathic pain and health related quality of life (HRQoL) in chronic SCI: a case study. Disabil. Rehabilit. Assist. Technol. 2016, 11, 529-534.

3. Sarasola-Sanz, A.; Irastorza-Landa, N.; López-Larraz, E.; Bibián, C.; Helmhold, F.; Broetz, D.; Birbaumer, N.; Ramos-Murguialday, A. A hybrid brain-machine interface based on EEG and EMG activity for the motor rehabilitation of stroke patients. In Proceedings of the 2017 IEEE International Conference on Rehabilitation Robotics (ICORR), London, UK, 17-20 July 2017; pp. 895-900.

4. Fontana, M.; Vertechy, R.; Marcheschi, S.; Salsedo, F.; Bergamasco, M. The body extender: A full-body exoskeleton for the transport and handling of heavy loads. IEEE Robot. Autom. Mag. 2014, 21, 34-44.

5. Papini, G.R.; Avizzano, C.A. Transparent force control for Body Extender. In Proceedings of the 2012 IEEE RO-MAN, Paris, France, 9-13 September 2012; pp. 138-143.

6. Anam, K.; Al-Jumaily, A.A. Active exoskeleton control systems: State of the art. Procedia Eng. 2012, 41, 988-994.

7. Kong, K.; Bae, J.; Tomizuka, M. Control of rotary series elastic actuator for ideal force-mode actuation in human-robot interaction applications. IEEE/ASME Trans. Mech. 2009, 14, 105-118.

8. Kim, S.; Bae, J. Force-mode control of rotary series elastic actuators in a lower extremity exoskeleton using model-inverse time delay control (MiTDC). IEEE/ASME Trans. Mech. 2017, 22, 1392-1400.

9. Kazerooni, H.; Racine, J.L.; Huang, L.; Steger, R. On the control of the berkeley lower extremity exoskeleton (BLEEX). In Proceedings of the 2005 IEEE International Conference on Robotics and Automation (ICRA), Barcelona, Spain, 18-22 April 2005; pp. 4353-4360.

10. Kazerooni, H.; Steger, R.; Huang, L. Hybrid control of the Berkeley lower extremity exoskeleton (BLEEX). Int. J. Robot. Res. 2006, 25, 561-573.

11. Kong, K.; Bae, J.; Tomizuka, M. A compact rotary series elastic actuator for human assistive systems. IEEE/ASME Trans. Mech. 2012, 17, 288-297.

12. Wang, L.; Du, Z.; Dong, W.; Shen, Y.; Zhao, G. Intrinsic Sensing and Evolving Internal Model Control of Compact Elastic Module for a Lower Extremity Exoskeleton. Sensors 2018, 18, 909.

13. Sankai, Y. HAL: Hybrid assistive limb based on cybernics. In Robotics Research; Springer: London, UK, 2010; pp. 25-34.

14. Suzuki, K.; Mito, G.; Kawamoto, H.; Hasegawa, Y.; Sankai, Y. Intention-based walking support for paraplegia patients with Robot Suit HAL. Adv. Robot. 2007, 21, 1441-1469.

15. Ghan, J.; Kazerooni, H. System identification for the Berkeley lower extremity exoskeleton (BLEEX). In Proceedings of the 2006 IEEE International Conference on Robotics and Automation (ICRA), Orlando, FL, USA, 15-19 May 2006; pp. 3477-3484.

16. Ko, J.; Fox, D. GP-BayesFilters: Bayesian filtering using Gaussian process prediction and observation models. Auton. Robots 2009, 27, 75-90.

17. Deisenroth, M.P.; Huber, M.F.; Hanebeck, U.D. Analytic moment-based Gaussian process filtering. In Proceedings of the 26th Annual International Conference on Machine Learning, Montreal, QC, Canada, 14-18 June 2009; pp. 225-232.

18. Deisenroth, M.P.; Turner, R.D.; Huber, M.F.; Hanebeck, U.D.; Rasmussen, C.E. Robust filtering and smoothing with Gaussian processes. IEEE Trans. Autom. Control 2012, 57, 1865-1871. [Green Version]

19. Eleftheriadis, S.; Nicholson, T.F.; Deisenroth, M.P.; Hensman, J. Identification of Gaussian Process State Space Models. arXiv, 2017; arXiv:1705.10888.

20. Deisenroth, M.P.; Ng, J.W. Distributed gaussian processes. arXiv, 2015; arXiv:1502.02843.

21. Rasmussen, C.E. Gaussian processes in machine learning. In Advanced Lectures on Machine Learning; Springer: Berlin, Germany, 2004; pp. 63-71.

22. Turner, R.; Deisenroth, M.; Rasmussen, C. State-space inference and learning with Gaussian processes. In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, Sardinia, Italy, 13-15 May 2010; pp. 868-875.

23. Ng, J.W.; Deisenroth, M.P. Hierarchical mixture-of-experts model for large-scale Gaussian process regression. arXiv, 2014; arXiv:1412.3078.

24. Cao, Y.; Fleet, D.J. Generalized product of experts for automatic and principled fusion of Gaussian process predictions. arXiv, 2014; arXiv:1410.7827.

25. Tresp, V. A Bayesian committee machine. Neural Comput. 2000, 12, 2719-2741.

26. Kitagawa, G. Monte Carlo filter and smoother for non-Gaussian nonlinear state space models. J. Comput. Graph. Stat. 1996, 5, 1-25.

27. Huang, R.; Cheng, H.; Guo, H.; Chen, Q.; Lin, X. Hierarchical interactive learning for a human-powered augmentation lower exoskeleton. In Proceedings of the 2016 IEEE International Conference on Robotics and Automation (ICRA), Stockholm, Sweden, 16-21 May 2016; pp. 257-263.

28. Arnold, D.V.; Hansen, N. Active covariance matrix adaptation for the (1 + 1)-CMA-ES. In Proceedings of the 12th Annual Conference on Genetic and Evolutionary Computation, Portland, OR, USA, 7-11 July 2010; pp. 385-392.

29. Arnold, D.V.; Hansen, N. A (1 + 1)-CMA-ES for constrained optimization. In Proceedings of the 14th Annual Conference on Genetic and Evolutionary Computation, Philadelphia, PA, USA, 7-11 July 2012; pp. 297-304.

30. Long, Y.; Du, Z.; Cong, L.; Wang, W.; Zhang, Z.; Dong, W. Active disturbance rejection control based human gait tracking for lower extremity rehabilitation exoskeleton. ISA Trans. 2017, 67, 389-397.