Next Article in Journal
A Biomimetic Design Method for 3D-Printed Lightweight Structures Using L-Systems and Parametric Optimization
Previous Article in Journal
Synchronization Sliding Mode Control of Closed-Kinematic Chain Robot Manipulators with Time-Delay Estimation
Previous Article in Special Issue
Optimal Assistance Timing to Induce Voluntary Dorsiflexion Movements: A Preliminary Study in Healthy Participants
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Human-Robot Interaction Torque Estimation Methods for a Lower Limb Rehabilitation Robotic System with Uncertainties

1
Electronics Engineering, Universidad Pontificia Bolivariana, Medellín 050031, Colombia
2
Electronics and Telecommunications Engineering Department, Universidad de Medellín, Medellín 050026, Colombia
3
Instituto de Ingeniería, Universidad Nacional Autónoma de México, Ciudad de México 04510, Mexico
4
ENSTA Bretagne, French State Graduate, 29200 Brest, France
5
School of Computing, Mathematics and Engineering, Charles Sturt University, Bathurst 2678, Australia
6
School of Science, Edith Cowan University, Joondalup 6027, Australia
*
Author to whom correspondence should be addressed.
Appl. Sci. 2022, 12(11), 5529; https://doi.org/10.3390/app12115529
Submission received: 20 March 2022 / Revised: 10 April 2022 / Accepted: 12 April 2022 / Published: 29 May 2022
(This article belongs to the Special Issue Assistive Technology: Biomechanics in Rehabilitation Engineering)

Abstract

:
Lower limb rehabilitation robot (LLRR) users, to successfully conduct isotonic exercises, require real-time feedback on the torque they exert on the robot to meet the goal of the treatment. Still, direct torque measuring is expensive, and indirect encoder-based estimation strategies, such as inverse dynamics (ID) and Nonlinear Disturbance Observers (NDO), are sensitive to Body Segment Inertial Parameters (BSIPs) uncertainties. We envision a way to minimize such parametric uncertainties. This paper proposes two human–robot interaction torque estimation methods: the Identified ID-based method (IID) and the Identified NDO-based method (INDO). Evaluating in simulation the proposal to apply, in each rehabilitation session, a sequential two-phase method: (1) An initial calibration phase will use an online parameter estimation to reduce sensitivity to BSIPs uncertainties. (2) The torque estimation phase uses the estimated parameters to obtain a better result. We conducted simulations under signal-to-noise ratio (SNR) = 40 d B and 20% BSIPs uncertainties. In addition, we compared the effectiveness with two of the best methods reported in the literature via simulation. Both proposed methods obtained the best Coefficient of Correlation, Mean Absolute Error, and Root Mean Squared Error compared to the benchmarks. Moreover, the IID and INDO fulfilled more than 72.2% and 88.9% of the requirements, respectively. In contrast, both methods reported in the literature only accomplish 27.8% and 33.3% of the requirements when using simulations under noise and BSIPs uncertainties. Therefore, this paper extends two methods reported in the literature and copes with BSIPs uncertainties without using additional sensors.

1. Introduction

According to the World Health Organization (WHO), in the last World Report on Disability, there are one thousand million people worldwide with some disability, and about 200 million have function disabilities. These people tend to have lower health and academic outcomes, lower economic participation, and higher poverty rates than people without disabilities [1]. Moreover, the Sustainable Development Goals (SDGs) reported that people living with disabilities are socially excluded. Therefore, they have fewer opportunities to improve their lives and reduce their poverty [2,3].
In the European Union, almost 45 million people aged between 15 and 64 had a disability during 2014, which corresponds to 14.1% of that age group [4]. Between 2001 and 2013, about 70 million people lived with disabilities in Latin America and the Caribbean, which is equivalent to a 12.5% of the regional population. The Pan American Health Organization (PAHO) reported that 6.3% of the people in Colombia had a state of disability, taking second place in Latin American countries after Brazil [5]. Additionally, 29% of those persons with disabilities have limitations to move or walk [6].
Multiple causes can trigger these types of limitations. In Colombia, one of the most disturbing causes is the high incidence of victims of antipersonnel mines (APM), improvised explosive devices (IED), and unexploded ordnance (UXO). There have been more than 11,700 APM and UXO victims since 1990. According to the Colombian National Development Plan 2018–2022, this situation is one of the most unfortunate consequences of the armed conflict in Colombia. These artifacts may have been activated for 20 years, and today, they are still hidden in several areas of the country, which makes it one of the main challenges to be faced in the post-conflict. Within the total of the victims, 20% died and the remaining 80% were injured, requiring physical rehabilitation during pre-prosthetic and prosthetic phases [7,8].
There are many traditional methods for rehabilitation purposes. However, the study of new technologies applied in automation and bioengineering has brought research and experimentation in robotic platforms for lower limb rehabilitation [9]. A lower limb rehabilitation robot (LLRR) may be used to conduct isotonic exercises for victims of APM, IED, and UXO. According to Hyunchoi et al. [10], an isotonic exercise is a type of resisted mode. A resistive method is an exercise where the patient is asked to overcome a given force during the movements imposed by the LLRR. Ideally, during isotonic exercises, the force of the subject is constant. In contrast, the muscle’s length changes, i.e., the joint angle changes, and the velocity changes. Such ideal isotonic exercise may be satisfied by quantifying the interaction torque that the subject performs. This variable is also defined as human–robot interaction torque. Force/torque sensors [11] could measure this torque, but it significantly adds to device cost. In the case of an LLRR for a developing country’s market, the costs of the sensors would be prohibitive. Moreover, these sensors add complexity and weight to the device unless miniaturized custom-made force/torque sensors are used [12]. Other approaches require using additional pressure sensors [13], Foot Contact Forces (FCFs) [14], or ground reaction forces (GRFs) [15]. Moreover, several approaches use Surface Electromyography (sEMG) sensors and signal processing techniques to estimate the human–robot interaction torque. However, using sEMG signals adds additional procedures and costs. In addition, they are also subject to variabilities inherent to the sEMG signal acquisition [16,17,18]. In addition, other studies report only the tests for a one-Degree-of-Freedom (DOF) LLRR, such as the estimation of ankle dynamic joint torque by a neuromusculoskeletal solver-informed Neural Network (NN) model [19]. These methods do not take into account the coupled dynamics of multiple joints. Saadatzi et al. [12] reported an Inverse Dynamics (ID) and Nonlinear Disturbance Observer (NDO) method to estimate the human–robot interaction torque without using additional sensors. As a result, their methods can produce accurate estimations. However, they rely on a precise model. Therefore, they reported the need to extend their methods for the cases when there are uncertainties. Aljuboury et al. [20] evaluated the behavior of three methods to control the angular position of an exoskeleton for knee assistance. They reported the controllers’ effectiveness and showed that a model reference adaptive control (MRAC) strategy based on an NDO obtained better performance than the others. Therefore, NDO-based methods are promising for human–robot interaction torque estimation when there are uncertainties.
A multi-DOF LLRR is a nonlinear system with coupled dynamics [21], which is subject to parametric uncertainties and external perturbations, e.g., there are uncertainties in the Body Segment Inertial Parameters (BSIPs) such as segment masses, Center of Mass (CoM), and inertias. Therefore, it is difficult to obtain an exact model of the system [22,23]. Moreover, according to Fahmy et al. [24], the dynamics of the human–robot system change from subject to subject. Additionally, it changes for one subject during the progress of their rehabilitation, i.e., between sessions. In cadaver studies, it is simple to locate the CoM by determining the center of balance of the segment [25]. However, the cross-sectional area and length profile are required to compute the CoM in vivo accurately, obtained using complex imaging systems. Therefore, several models are used to estimate the BSIPs indirectly. The models reported by Plagenhoef et al. [26] and De Leva et al. [27] contain each body segment weight and length as the percentage of total body weight and height, respectively. Moreover, the CoM location is given as a percentage of the segment length. Both methods may be helpful to obtain a BSIPs rough estimate. However, each body segment has a unique bone, muscle, fat, and other tissues combination [25]. Therefore, the density of a segment is not uniform, and models such as the reported by Plagenhoef et al. [26], and De Leva et al. [27] do not accurately model the LLRR with the user. Thus, developing a human–robot interaction torque estimation method is challenging to perform the isotonic rehabilitation exercises with the LLRR. Furthermore, subjects with partial lower limb amputations are not included in the aforementioned models. Therefore, it is even more challenging to compute their BSIPs. Moreover, adding a prosthesis to their limbs would also modify the BSIPs values in an unknown manner, adding uncertainties to the model.
The estimation problem concerns the use of the dynamic model of an LLRR to estimate the human–robot interaction torque τ ^ i n t despite BSIPs uncertainties. To the best of the authors’ knowledge, none of the previous research has examined how a human–robot interaction torque estimation method for an LLRR, usable for lower limb isotonic exercises, may be used as biofeedback to favor isotonicity, despite the BSIPs uncertainties. Therefore, this paper extends the ID-based and NDO-based methods reported by Saadatzi et al. [12]. We propose the Identified ID-based method (IID) and the Identified NDO-based method (INDO).
The following statements emphasize this study’s significant contributions:
  • We defined a list of 18 requirements of a human–robot torque estimation method for an LLRR, usable for isotonic. We defined the requirements with a literature review and a survey-based methodology.
  • This paper proposes two human–robot interaction torque estimation methods, the IID and INDO. Both methods must cope with two challenging issues: (a) BSIPs uncertainties exist in the subject model, and (b) no force or additional sensors are to be used. Both methods do not require a physiotherapist to make an exact measurement of the BSIPs of the patient’s limbs, but uses approximate values computed in terms of total height and body weight. Finally, the methods would take the data from the first iterations to reduce the sensitivity to BSIPs uncertainties.
  • Our proposal avoids the bidirectional coupling between identifier, estimator, and controller to guarantee convergence in each stage separately. For this purpose, the robot performs a persistent excitation trajectory to identify the parameters. The optimized trajectory may be used for any given patient regardless of their BSIPs. Then, we turn off the parameter identifier in a second phase, and the robot executes the rehabilitation therapy movements, estimating the torque exerted by the subject.
The organization of the paper is as follows. Section 2 describes the requirements for the torque estimation method. Section 3 and Section 4 describe the model of the system and the torque interaction method reported by Saadatzi et al. [12], respectively. Section 5 presents two methods to estimate the human–robot interaction torque despite uncertainties in the BSIPs. Section 6 shows the effectiveness of the IID and INDO by performing simulations on a 3-DOF LLRR called Nukawa, see Appendix A for a detailed explanation of the LLRR. Finally, some conclusions and future work are presented in Section 7.

2. Requirements Definition

The requirements for the human–robot interaction torque estimation method were defined by applying a survey to a group of physiotherapists and conducting a literature review.
The idea was to ask people in the health area about the requirements and to define the movements for the LLRR, i.e., the isotonic exercises to be used within the simulation study. Appendix B describes the study in more detail. In addition, it specifies the questions and answers used in the questionnaire. A total of 26 physiotherapists answered the anonymous questionnaire, and we asked about the Range of Motion (ROM), speed, and force to be used within the simulations and maximum tolerable error for the human–robot torque estimation methods. The survey results show that the physiotherapists defined the squat and press exercises to evaluate the performance of the human–robot interaction torque estimation methods. Moreover, they suggested four additional requirements for these methods.
The methods used during the information search strategy are described below. Additionally, the methodology used to select studies based on inclusion and exclusion criteria is presented. Finally, the steps that were used to extract information are reported.
A literature review was carried out to define the requirements of a human–robot torque estimation method for a robotic system usable for isotonic exercises. The review was conducted using the primary databases available, such as Scopus, ScienceDirect, IEEE Xplore, and Google Scholar, among others. This review was conducted, covering publications mainly between 2006 and 2022. The search was initially oriented to find the main metrics and characteristics, which are currently employed to evaluate the behavior of these methods. A list of keywords or phrases was defined to be used within the search. The search strategy included combining several Boolean operators to filter the results. Only literature in English and Spanish was reviewed. The keywords were: (a) Human-robot, (b) Interaction, (c) Torque, (d) Estimation, (e) Requirements, (f) Force, (g) Measurements (h) Sensors, (i) BSIPs, (j) Uncertainty, (k) Rehabilitation.
Inclusion and exclusion criteria were also defined for the study selection. Inclusion criteria were (a) studies published between January 2006 and March 2022, and (b) full access to the research text. Exclusion criteria included: (a) publications that do not fit the subject of study of this review, (b) publications in languages other than Spanish and English, (c) letters to the editor, and (d) opinion articles. We extracted the essential data of each of the publications in a spreadsheet.
Finally, let us enumerate the 18 chosen requirements of a human–robot interaction torque estimation method for an LLRR, which is usable for isotonic exercises. The literature review inspired the first 13, and the physiotherapists’ answers inspired the last five:
  • Non-dependence of additional sensors [11,13].
  • Low phase lag in the estimation [15].
  • Low sensor noise sensitivity [15].
  • Low sensitivity to BSIPs uncertainties [12].
  • The average percentage error must be lower than 20 to 22% when using accurate model parameters [12].
  • Small error band [10].
  • Approximately 0.5 s of settling time or lower [10].
  • Overshoot of approximately 25% or lower within the estimation [10].
  • Coefficient of Correlation ( R 2 ) greater or approximately 0.935 for the hip joint and 0.924 for the knee joint [11].
  • Root Mean Squared Percent Error (RMSPE) lower than 8.74% for the hip joint and 10.26% for the knee joint [11].
  • A maximum of 5% error when moving just one joint, i.e., the distal one [28].
  • Finite-time convergence [29].
  • It should not require calibration each time that the user wears the LLRR, i.e., it requires a maximum of one calibration per user. [13].
  • It works in all ROM, and the limiting angles of the joints must be configurable.
  • It works with slow exercises, executing between 1–25 repetitions per set.
  • It works within the following ranges of forces: 0 kg to 15 kg for the hip, 0 kg to 15 kg for the knee, and 0 kg to 10 kg for the ankle.
  • It works having a maximum percentage of error in the range from 1% to 3%.
  • It works within a squat exercise.

3. Dynamic Model

The process for the direct kinematics was based on the methodology reported by Siciliano et al. [30]. The LLRR’s computer model was derived through an analytical method, taking the kinematic and dynamic equations into account. We used the Euler–Lagrange method as a systematic energy-based approximation that uses Kinetic Energy (KE) and Potential Energy (PE) [31,32].
Figure 1 presents the LLRR configuration, where i are the segment’s lengths with i = { 1 , 2 , 3 } , b i are the position of the CoM for each limb, m i are the masses of each link, q i the angles of each segment with respect to the previous segment, and g 0 the gravity acting on the vertical plane. Finally, B r is the backrest angle, for simulation purposes. Table 1 reports the Denavit Harbenbert (D-H) parameters for the LLRR using the standard (classic or distal) notation [33].
We can write the dynamics for the LLRR as
τ a = M ( q ) q ¨ + C ( q , q ˙ ) + G ( q ) τ i n t τ f r i
where q R 3 × 1 , q ˙ R 3 × 1 , and q ¨ R 3 × 1 represent the joint positions, velocities, and accelerations, respectively. M ( q ) R 3 × 3 is the inertia matrix, C ( q , q ˙ ) R 3 × 1 is the coriolis and centripetal vector, G ( q ) R 3 × 1 gravity vector, F ( q ˙ ) R 3 × 1 the friction term [32], τ a R 3 × 1 is the torque of the actuators, and τ i n t R 3 × 1 are the disturbance torques, e.g., the torque exerted by the subject. Appendix C presents detailed information of the dynamics.

4. Original ID-Based and NDO-Based Methods

The objective of this section is to present the original ID-based and NDO-based methods and their limitations. Let us introduce the formulation of the ID-based method reported by Saadatzi et al. [12]:
τ ^ d = M ^ ( q ) q ¨ + C ^ ( q , q ˙ ) + G ^ ( q ) τ a τ ^ f r i = F v q ˙ τ ^ i n t = τ ^ d τ a τ ^ f r i
It is important to note that this ID-based method uses double derivatives to compute q ¨ , which significantly increases the sensitivity of the model to sensor noise. Moreover, this method uses the estimates M ^ ( q ) , C ^ ( q , q ˙ ) , and G ^ ( q ) of the actual M ( q ) , C ( q , q ˙ ) , and G ( q ) dynamic matrices, respectively. Taking into account that
M ( q ) = M ^ ( q ) + Δ M ( q ) ,
C ( q , q ˙ ) = C ^ ( q , q ˙ ) + Δ C ( q , q ˙ ) ,
and
G ( q ) = G ^ ( q ) + Δ G ( q )
where Δ M ( q ) , Δ C ( q , q ˙ ) , and Δ G ( q ) are deviations of the dynamic matrices due to BSIPs uncertainties, these estimates may turn into sensibility of the method to BSIPs uncertainties. These dynamic matrices M ^ ( q ) , C ^ ( q , q ˙ ) , and G ^ ( q ) may be computed using an initial estimate of the parameters of the LLRR. However, they may have uncertainties when an amputee subject is using the LLRR, e.g., in the masses, CoM, or Inertias.
Now, let us introduce the formulation of the NDO-based method reported by Saadatzi et al. [12,34,35]. This NDO has the following formulation:
z ˙ = L ( q , q ˙ ) z + L ( q , q ˙ ) { C ^ ( q , q ˙ ) + G ^ ( q ) τ a p ( q , q ˙ ) } τ ^ d = z + p ( q , q ˙ ) d d t p ( q , q ˙ ) = L ( q , q ˙ ) M ^ ( q ) q ¨
where L ( q , q ˙ ) is the observer gain matrix, τ ^ d is an estimation of the disturbances in a lumped term, M ^ ( q ) , C ^ ( q , q ˙ ) , and G ^ ( q ) are initial estimates of the dynamic matrices, z and p ( q , q ˙ ) are an auxiliary variable and vector to avoid the differentiation of position measurements more than once, since it may add problems as mentioned before.
In order to design the observer, the challenge is to determine L ( q , q ˙ ) and p ( q , q ˙ ) . According to Mohammadi et al. [34], the NDO gain matrix L ( q , q ˙ ) can be computed taking into account that it depends only on q as
L ( q ) = X 1 M ^ 1 ( q )
where X is an n × n invertible and constant matrix to be determined. Additionally, remember that the robot’s inertia matrix is symmetric and positive definite and, thus, invertible. Then, by substituting (7) into (6), we obtain p in terms of q ˙ , as shown below in the disturbance observer auxiliary vector:
p ( q ˙ ) = X 1 q ˙ .
Saadatzi et al. [12] evaluated five methods for human–robot interaction torque estimation. These methods ranged from relying only on commanded actuation torques, i.e., motor currents, to techniques that use the whole dynamics of the system, such as the ID-based method and the NDO-based method. The ID-based method obtained 20 to 22% average error, while the NDO-based method obtained 12 to 18% average error when using accurate model parameters. Both methods obtained better results than the other three methods evaluated, which got more than 25% average error. Moreover, they analyzed the sensibility of the ID and NDO methods to inaccuracies in model parameter estimations. This analysis shows considerable sensitivity of both approaches to model parameter variations. Moreover, they evidenced the behavior of both estimation methods with sensor noise. The noise in the estimations with the NDO method was smaller than the one obtained with the ID method. Therefore, one of the main differences between both methods is noise reduction. In addition, the NDO method is a filtered version since the nature of the observer can work without double derivates to compute the acceleration. To sum up, Saadatzi et al. [12] demonstrated how these systems work and state that their method could adequately estimate human–robot interaction. However, they also report that both methods require an accurate system model. Therefore, we wanted to take advantage of the fact that there was a benchmark to compare their study with our extended methods. Consequently, we based our approaches on their methods and used them for lower limbs. Moreover, we used more degrees of freedom and reduced the sensitivity to BSIPs uncertainties. In addition, they reported only the effects of constant offsets added to the elements of model matrices dynamic matrices, not the BSIPs by themselves, which is a more realistic approach. Finally, they only tested the sensitivity of the methods to variations in the parameter mass of the second limb, not all BSIPs with uncertainties at once. For those above, we selected these two approaches, i.e., the ID and NDO method reported by et al. [12].
Saadatzi et al. [12] evaluated the behavior of their ID-based method (2) and NDO-based method (6) in the presence of sensor noise and uncertainties. They reported that both methods are suitable to estimate the human–robot interaction torque. However, both approaches failed to estimate the human–robot interaction torque τ ^ i n t when there are uncertainties. Therefore, this paper extends the ID-based and NDO-based methods to cope with BSIPs uncertainties.

5. Proposed Methods

In this section, we propose two methods, the IID and INDO, that extend the work reported by Saadatzi et al. [12]. The IID and INDO estimate the human–robot interaction torque despite practically expected uncertainties in the BSIPs of an amputee model victim of APM, IED, and UXO who would use an LLRR. These methods are comprised of two phases. Moreover, these phases are executed one after the other, not simultaneously. (1) Phase one is the calibration phase, which allows an online parameter estimation to reduce BSIPs uncertainties. (2) Phase two is the estimation phase, where the human–robot interaction torque uses the estimated parameters under the assumption that the BSIPs uncertainties have been reduced. The objective of the first phase is to extend their method and reduce BSIPs uncertainties to improve the performance of the two human–robot torque estimation methods. This calibration phase is based on the online parameter estimation algorithm reported by Riani et al. [36]. The combinations will be labeled as the IID and the INDO. In the remaining of this section, we report the design process for the IID and INDO.
This section is divided into three subsections. Section 5.1 presents the model (1) transformed as a linear in the parameters system, that will be used for the estimation of the parameters. Section 5.2 presents the calibration phase, i.e., the online parameter estimation method and a discussion of it. Section 5.3 presents the IID and INDO, with the addition of the calibration phase, i.e., the online parameter estimation method.
In the IID and INDO, the controller works independently. Therefore, our proposal does not have two-way dynamic couplings between the identifier, estimator, and controller. That is, there is no feedback on the estimate to the controller. Therefore, it is not necessary to consider the principle of separation of the observed-controller system.

5.1. Transforming the Model

The objective of this section is to present the model (1) transformed as a linear in the parameters system that will be used for the identification of the parameters. Consider the LLRR reported in Section 3. Each link of the LLRR is characterized by ten dynamic parameters expressed as
m i , r C i = r x i r y i r z i , I C i = I C i , x x I C i , x y I C i , x z I C i , x y I C i , y y I C i , y z I C i , x z I C i , y z I C i , z z .
However, the LLRR’s dynamics depend in a nonlinear way on some of these parameters, e.g., through the combination I c i , z z + m i r x i 2 . Many standard parameters do not play a role in the dynamic model of the given LLRR. Therefore, one can isolate p < 10 independent groups of parameters. These grouped parameters are called dynamic coefficients. These parameters are the only ones that affect robot dynamics. The base inertial parameters or identifiable parameters are the minimum inertial parameter set used to completely characterize and obtain a dynamic model of a robot. These parameters can be obtained from the standard dynamic parameters by eliminating all the parameters that have no effect in the dynamic model and grouping some of them in linear combinations as reported by Khalil et al. [37]. The introduction of base inertial parameters χ i R p × 1 is a convenient regrouping of the dynamic parameters, i.e., a linear parametrization of LLRR’s dynamics. Therefore, we should obtain a vector χ that contains the base inertial parameters. Thus, the dynamic model (1) can be written such that the parameters enter in a linear fashion as
τ = W ( q , q ˙ , q ¨ ) χ
where W ( q , q ˙ , q ¨ ) R N × p is the regressor matrix. It is important to note that it will be denoted as W in the rest of this report to simplify the notation.
Next, we will present the adjustment for the mathematical model of the LLRR (1) to the form (10) indicated by the formulation reported by Riani et al. [36]. Therefore, we assumed a set of parameters χ i as
χ 1 = I c 1 , z z + I c 2 , z z + I c 3 , z z + L 1 2 m 2 + L 1 2 m 3 + L 2 2 m 3 + b 1 2 m 1 + b 2 2 m 2 + b 3 2 m 3 χ 2 = L 1 g 0 m 2 + L 1 g 0 m 3 + b 1 g 0 m 1 χ 3 = m 3 L 2 2 + m 2 b 2 2 + m 3 b 3 2 + I c 2 , z z + I c 3 , z z χ 4 = L 1 L 2 m 3 + L 1 b 2 m 2 χ 5 = L 2 g 0 m 3 + b 2 g 0 m 2 χ 6 = m 3 b 3 2 + I c 3 , z z χ 7 = L 2 b 3 m 3 χ 8 = L 1 b 3 m 3 χ 9 = b 3 g 0 m 3
See Table 2 for the definition of all these parameters. Subsequently, we obtained the transformed matrices M ( χ , q ) , C ( χ , q , q ˙ ) , and G ( χ , q ) defined as
M ( χ , q ) = m 11 ( χ , q ) m 12 ( χ , q ) m 13 ( χ , q ) m 21 ( χ , q ) m 22 ( χ , q ) m 23 ( χ , q ) m 31 ( χ , q ) m 32 ( χ , q ) m 33 ( χ , q )
C ( χ , q , q ˙ ) = C 1 ( χ , q , q ˙ ) C 2 ( χ , q , q ˙ ) C 3 ( χ , q , q ˙ )
and
G ( χ , q ) = g 1 ( χ , q ) g 2 ( χ , q ) g 3 ( χ , q )
by using the same process implemented in Section 3. See Appendix D for a detailed description of the matrices M ( χ , q ) , C ( χ , q , q ˙ ) , and G ( χ , q ) .
Subsequently, we combined these parameters to get the vector χ of the base inertial parameters as given in Gautier et al. [38]. Therefore, we obtained a vector of the base inertial parameters as
χ = χ 1 χ 2 χ 3 χ 4 χ 5 χ 6 χ 7 χ 8 χ 9 T
with χ i with i = 1 , 2 , , 9 defined in (11). This linear parametrization of robot dynamics is not unique in the chosen set of dynamic coefficients χ and the associated W . The W is a 3 × 9 matrix, where n = 3 is the number of DOF of the LLRR and, in this case, p = 9 is the number of the base inertial parameters. Therefore, the dynamic model (1) can be written as (10) where W ( q , q ˙ , q ¨ ) is the regressor matrix, and its notation may be simplified as
W ( q , q ˙ , q ¨ ) = w 11 w 12 w 13 w 14 w 15 w 16 w 17 w 18 w 19 0 0 w 23 w 24 w 25 w 26 w 27 w 28 w 29 0 0 0 0 0 w 36 w 37 w 38 w 39
See Appendix D for a detailed description of the elements of (16). In this equation, we can note that the regressor matrix has a block upper triangular structure as reported by Siciliano et al. [30]. This triangular block structure may enable us to compute the parameter estimation using a sequential procedure. In this case, the parameter estimation can be executed by an iterative method, taking measurements from the distal limb to the proximal. However, such a technique may have the drawback of accumulating any error due to ill-conditioning of the matrices involved step-by-step. It may then be worthwhile to use a global approach that imposes motions on all LLRR joints at once [30]. The calibration phase presented in Section 5.2 will show a multi-DOF approach for an online estimation to reduce BSIPS uncertainties.

5.2. Calibration Phase Design

The objective of this section is to present the calibration phase. This calibration phase reduces BSIPs uncertainties based on an online parameter estimation method. Taking into account the model presented in (10), based in the regressor matrix W ( q , q ˙ , q ¨ ) and the vector χ , we used the observer reported by Riani et al. [36] to perform an online parameter estimation. This method is expressed as
Γ ^ ˙ = W χ ^ + α 2 ( Γ Γ ^ ) χ ^ ˙ = K W T ( Γ Γ ^ ) + K W T ( τ W χ ^ ) K ˙ = 2 K W T W K + α K
where
  • Γ = τ d t R 3 × 1 is the integral of τ .
  • χ ^ R p × 1 and Γ ^ R 3 × 1 are, respectively, the estimates of the vectors χ and Γ .
  • K R p × p is a bounded symmetric positive definite matrix.
  • α > 0 is a positive gain.
The following assumptions are required for the estimator:
  • The vector τ is measurable.
  • W satisfies the persistent excitation condition [39].
  • χ is a constant vector.
Let us take into account the new dynamic model for our LLRR expressed as (15) and (16). Moreover, the observer reported by Riani et al. [36] is defined as (17). According to Riani et al. [36], if we define the vector estimation error as Γ ˜ = Γ Γ ^ and the vector parameters error as χ ˜ = χ χ ^ , then we can compute their dynamics as
Γ ˜ ˙ = W χ ˜ α 2 Γ ˜ χ ˜ ˙ = K W T Γ ˜ K W T W χ ˜
Please refer to Riani et al. [36] for detailed stability analysis. In their document, the Lyapunov-like function decreases exponentially to zero, and the errors Γ ˜ and χ ˜ tend exponentially to zero since W satisfies the persistent excitation condition. To ensure the stability of the online parameter estimation, we can define the LMI [36] as:
H = γ I p γ 2 I p γ 2 I p α 2 K 1 > 0
where γ R is a positive gain, K R p × p is a symmetric positive definite matrix, and I p R p × p is the identity matrix. Therefore, the design for the calibration phase reduces to finding the condition on K , γ , and α such that H > 0 .

5.3. Estimation Phase Design

Phase two is the estimation phase, where the human–robot interaction torque uses the estimated parameters under the assumption that the BSIPs uncertainties have been reduced. The main objective of this section is to present the IID and INDO with the addition of the online parameter estimation method.
On the one hand, we extend the ID-based method reported by Saadatzi et al. by reducing the need for an accurate model since BSIPs uncertainties are reduced using online parameter estimation. Figure 2 presents a block diagram for the IID. This figure depicts that the IID is composed of two phases. Moreover, these phases are executed one after the other, not simultaneously. Figure 2a presents the calibration phase that uses an online parameter estimation to reduce BSIPs uncertainties. Thus, estimating χ ^ ˙ . During this phase, the τ i n t = 0 , i.e., the subject should be relaxed without exerting any torque on the robot. The first phase is executed during t < t i , where t means the simulation time and t i is the time where the identification ends, i.e., after convergence. Figure 2b presents the estimation phase that uses a transformed version of the ID along with the estimated parameters χ ^ to estimate the human–robot interaction torque. The estimation phase is executed while t t i . During this phase, the τ i n t 0 , i.e., the subject should is exerting torque on the robot and the goal of the phase is to estimate τ ^ i n t . Here, we can observe that the IID is defined as
  • Calibration phase t < t i
    Γ ^ ˙ = W χ ^ + α 2 ( Γ Γ ^ ) χ ^ ˙ = K W T ( Γ Γ ^ ) + K W T ( τ W χ ^ ) K ˙ = 2 K W T W K + α K
  • Torque estimation phase t t i
    τ ^ d = M ^ ( χ ^ , q ) q ¨ + C ^ ( χ ^ , q , q ˙ ) + G ^ ( χ ^ , q ) τ a τ ^ f r i = F v q ˙ τ ^ i n t = τ ^ d τ a τ ^ f r i
On the other hand, we extend the NDO-based method reported by Saadatzi et al. by reducing the need for an accurate model since BSIPs uncertainties are reduced using online parameter estimation. Figure 3 presents a block diagram for the INDO. This figure depicts that the INDO is composed of two phases. Moreover, these phases are executed one after the other, not simultaneously. Figure 3a presents the calibration phase that uses an online parameter estimation to reduce BSIPs uncertainties. Thus, estimating χ ^ ˙ . During this phase the τ i n t = 0 , i.e., the subject should be relaxed without exerting any torque on the robot. The first phase is executed during t < t i , where t means the simulation time and t i is the time where the identification ends, i.e., after convergence. Figure 3b presents the estimation phase that uses a transformed version of the NDO along with the estimated parameters χ ^ to estimate the human–robot interaction torque. The estimation phase is executed while t t i . During this phase, the τ i n t 0 , i.e., the subject should is exerting torque on the robot and the goal of the phase is to estimate τ ^ i n t .
Here, we can observe that the INDO is defined as
  • Calibration phase t < t i
    Γ ^ ˙ = W χ ^ + α 2 ( Γ Γ ^ ) χ ^ ˙ = K W T ( Γ Γ ^ ) + K W T ( τ W χ ^ ) K ˙ = 2 K W T W K + α K
  • Torque estimation phase t t i
    z ˙ = L ( q , q ˙ ) z + L ( q , q ˙ ) { C ^ ( χ ^ , q , q ˙ ) + G ^ ( χ ^ , q ) τ a p ( q , q ˙ ) τ ^ d = z + p ( q , q ˙ ) d d t p ( q , q ˙ ) = L ( q , q ˙ ) M ^ ( χ ^ , q ) q ¨ L ( q ) = X 1 M ^ 1 ( χ ^ , q ) p ( q ˙ ) = X 1 q ˙ τ ^ i n t = τ ^ d τ a τ ^ f r i
We can observe that the INDO uses the NDO-based reported by Saadatzi et al. [12] in terms of χ ^ . Therefore, the output of the calibration phase is used as an input for the estimation phase. This method uses the online parameter estimation to reduce BSIPs uncertainties. Subsequently, within the estimation phase, the NDO-based human–robot interaction torque uses the estimated parameters under the assumption that the BSIPs uncertainties have been reduced.
Let us now cite the following Theorem reported by Mohammadi et al. [34]. This Theorem is helpful since it presents sufficient conditions for asymptotic and exponential disturbance tracking when the robot shows fast-varying disturbances. Moreover, it evidences asymptotic and exponential convergence to estimate the NDO-based human–robot interaction torque estimation method.
Theorem 1
(Mohammadi et al. [34]). Consider the robotic manipulator subject to disturbances described by (1). The disturbance observer is given in (6) with the disturbance observer gain matrix L ( q ) defined in (7) and the disturbance observer auxiliary vector p ( q ˙ ) defined in (8). The disturbance tracking error τ d is globally uniformly ultimately bounded if:
  • The matrix X is invertible,
  • There exists a positive definite and symmetric matrix Γ such that
    X + X T X T M ^ ˙ ( q ) X Γ .
  • The rate of change of the lumped disturbance is bounded, i.e., κ > 0 such that τ ˙ d ( t ) κ for all t > 0 .
Under the above conditions and for all Δ τ d ( 0 ) R n , the tracking error converges with an exponential rate, equal to ( 1 q ) λ min ( Γ ) / 2 σ 2 X 2 to the ball with radius 2 κ σ 2 X 2 / q λ min ( Γ ) where 0 < q < 1 .
Where λ min ( · ) denotes the minimum eigenvalue of a matrix, the LLRR velocity vector lies in a bounded set, i.e., q ˙ D q ˙ = q ˙ R n : q ˙ q ˙ max , and M ( q ) v 2 ( q ) where v 2 = sup q D q σ 2 ( q ) .
Please refer to Mohammadi et al. [34] for a detailed proof of Theorem 1. In their paper, they show that the tracking error converges with an exponential rate, equal to ( 1 q ) λ min ( Γ ) / 2 σ 2 X 2 to the ball with radius 2 κ σ 2 X 2 / q λ min ( Γ ) where 0 < q < 1 for all Δ τ d ( 0 ) R n .
According to Theorem 1, the NDO design reduces to finding a constant invertible matrix X that satisfies the inequality (24). The following theorem, reported by Mohammadi et al. [34], states how (24) can be formulated as a linear matrix inequality (LMI) problem. It is important to mention that this theorem is useful to design the NDO for the estimation phase.
Theorem 2
(Mohammadi et al. [34]). Define the matrix Y = X 1 and assume that an upper bound of M ^ ˙ ( q ) is ζ. The inequality (24) holds if the following LMI is satisfied:
Y + Y T ζ I Y T Y Γ 1 0
Please refer to Mohammadi et al. [34] for a detailed proof of Theorem 2.

6. Simulation Results

This section will compare the IID and the INDO to define which one is the best approach for the LLRR. Therefore, we will evaluate if these methods can estimate the torque of the human lower limb joints even with BSIPs uncertainties. Thus, extending the methods reported by Saadatzi et al. [12].

6.1. Simulation Parameters

Table 2 presents the nominal LLRR parameters for simulation, considering a subject of height h = 1.75 m and weight m = 75 kg. To define the model parameters reported in this table, we used the CAD design to obtain its masses m i r and lengths i r . In order to define the lengths i s and weights m i s of the subject’s segments, models reported by Plagenhoef [26] and De Leva [27] were used. These models contain anatomical data for analyzing human motion, and we used the average value between male and female subjects. Segment weight m i was expressed as a percentage of total body weight m as reported by De Leva et al. [27]. Segment length i s were expressed as a percentage of total height h as reported by Plagenhoef et al. [26]. These parameters were used as the nominal values. However, we mentioned before that they are not the actual values when an amputee subject with protheses is analyzed. The backrest B r of the LLRR is not an actuated joint. However, it is imperative to report the B r of the robotic system for each test since it determines in which position the robotic system is. Figure 1 presents the conventions for the B r . A counterclockwise movement is the convention for a positive arc of movement, and B r is measured taking into account the horizontal plane. The ROM of the three joints of the lower limbs during the simulations should be within the values presented by Berryman et al. [40].
In this work, the viscous friction vec { v i } was assumed in a heuristic way, taking into account that the simulation showed a three-segment pendulum behavior when the control algorithm is deactivated. It was assumed in this way so that the model behaved as a pendulum that oscillated one time before stabilizing when it was dropped from the position q = ( 0 , 0 , 90 ) , as expected in the LLRR. The method consisted of recording a video of the actual LLRR in the laboratory and using Kinovea, free and open-source software used for the motion analysis of sports techniques and gestures, used mainly by sports coaches and athletes to explore, study or comment on performance. Additionally, a heuristic method was used to tune the gains of the PID controller, i.e., the outer loop feedback. These gains were configured to reduce the error and the oscillations. The sampling period was set to TS = 0.001   s .
Regarding the moment of inertia, we define that it would be calculated taking into account a bar centered on the CoM, and with length 2 b i , where b i is the distance from the joint i to the CoM of said segment i. According to Raymond et al. (2004) [41], the moment of inertia for a thin bar is I = 1 3 M L 2 . In this case, the thickness of the bar is not taken into account. In summary, for us, the inertia of each segment would be I c i , z z = 1 3 m i ( 2 b i ) 2 . Finally, the position of the CoM b i was configured as reported by De Leva et al. [27].

6.2. Control Algorithm

The control algorithm is necessary for the simulation. We require the robot to move in a certain way so that it allows us to excite the system and its dynamics to estimate the parameters or the human-robot torque. Therefore, we have used the Computed Torque Control (CTC) algorithm reported in previous works by Yepes et al. [42,43].

6.3. Calibration Phase Results

Let us take into account that the model for the LLRR has nine standard dynamic parameters, i.e., size of χ . Moreover, let us suppose an α = 1 , which defines the exponentially decreasing convergence ratio for the calibration phase. This value was selected taking into account that the convergence was obtained in less than 25 s. It was estimated that an initial phase of 30 s was appropriately short for a patient to accept being idle without executing any interaction torque. Therefore, the simulations sought to calibrate the convergence to take less than 30 s , achieving good results with 25 s . Moreover, a more significant α would reduce convergence time. However, this would turn into sensitivity within the calibration phase since a more considerable observer gain produces more fluctuations of the value.
In addition, let us define the identity matrix as I 9 , and (19) would turn into
H = γ I 9 γ 2 I 9 γ 2 I 9 1 2 K 1 > 0
Subsequently, we used the MATLAB® Robot Control Toolbox to solve this LMI-based analysis and, thus, design the calibration phase estimator. The solution with α = 1 for the LMI (26) was γ = 15.691 and K = 0.0212 I 9 , i.e., satisfying the LMI with H > 0
The trajectory used within the calibration phase must be selected with caution to improve the convergence rate of this phase and reduce noise sensitivity. The abovementioned trajectory should be a persistent exciting trajectory exciting trajectory and can be computed using some optimization criteria as reported by Khalil et al. [37]. Therefore, to ensure that W satisfies the persistent excitation condition, the robot’s trajectory was optimized for the calibration phase. This selection of the persistent exciting trajectory was conducted by calculating the movements of the robot whose points give a well-conditioned observation matrix [37]. The condition number κ 2 ( W ) of a matrix W allows evaluating its well or ill-conditioning. A problem with a low condition number may be defined as well-conditioned, while a problem with a high condition number is said to be ill-conditioned [44]. The κ 2 ( W ) is defined using the 2-norm, as
κ 2 ( W ) = σ max σ min 1
where σ max and σ min are defined as the maximum and minimum singular values of W . Therefore, the nonlinear optimization problem consisted of determining a trajectory that provides a κ 2 ( W ) that is close to one [37]. For this purpose, we used a Genetic Algorithm (GA) algorithm with the MATLAB® Robot Control Toolbox to optimize a trajectory whose points gave a well-conditioned observation matrix. The code minimizes the κ 2 ( W ) by generating multiple trajectories to ensure the persistently exciting condition, i.e., the κ 2 ( W ) is approximately one. The code uses GA for this purpose and creates trajectories using a combination of two sinusoid signals per joint for the q d . The frequencies were restricted to values lower than 0.2 Hz to take into account physical restrictions for the LLRR. Moreover, the maximum number of generations was configured as 10. The initial position was set to q0 = (0, −90, 90)° with no speed. Moreover, the amplitude for each sinusoid was set to 30°, i.e., these signals indicate the LLRR movement in degrees concerning the origin. The obtained persistent exciting trajectories are defined as
q 1 d = 30 sin ( 0.0430 t ) + 30 sin ( 0.2316 t ) q 2 d = 60 sin ( 0.0938 t ) + 30 sin ( 0.0934 t ) q 2 d = 30 sin ( 0.0594 t ) + 30 sin ( 0.1375 t ) .
These optimized trajectories generate a well-conditioned observation matrix, since κ 2 ( W ) is optimized from values with a 10 35 magnitude to κ 2 ( W ) = 30.05 , which is closer to 1. Therefore, our problem has a low condition number and is well-conditioned. Thus, the persistent excitation is satisfied.
Figure 4 depicts the parameter estimation convergence. Here, the method estimates χ ^ , and the convergence takes around 20 s. We can see that the method has slow dynamics, as expected with an α = 1 . Moreover, as desired, there is no chattering after converging. The nominal dynamic parameters χ n o m of our LLRR along with the new subject, i.e., with the simulation parameters reported in Table 2, are
χ n o m = 10.0418 148.1905 3.8831 3.2052 74.6331 0.5343 0.7208 0.6976 16.2432 T .
However, the initial conditions χ ^ 0 , using Table 2 and 20% uncertainties in the BSIPs are
χ ^ 0 = 12.83 189.216 5.0249 4.0885 95.2003 0.7621 1.0379 1.0045 23.3902 T .
After the calibration phase, the estimated parameters are
χ ^ = 10.0395 148.2095 3.8878 3.177 74.6555 0.5428 0.7223 0.7091 16.2991 T
and we can observe that (31) is very similar to (29), despite using 20% uncertainties in the BSIPs within the simulations, as suggested by Grun et al. [15]. The online parameter estimator reduces the percentage of uncertainty in the base inertial parameters. All nine parameters of χ n o m are estimated in the simulated 25 s window, and these values of χ ^ are to be used within the second phase of the IID and INDO. See Section 6 for the results.

6.4. Estimation Phase Results

Based on the parameters from Table 2 for the simulation, we can see that
M ^ ˙ ( q ) ζ = 30
Finally, with the identity matrix I 3 R 3 × 3 and (25) would turn into
Y + Y T 30 I 3 Y T Y Γ 1 0 .
Subsequently, we used the MATLAB® Robot Control Toolbox to solve this LMI-based analysis, obtaining a design for the NDO. The solution with ζ = 30 for (33) was Y = X = 0.0028 I and Γ = 0.9542 × 10 3 I . Hence, satisfying the LMI. Therefore, we successfully designed the NDO for the human–robot interaction torque estimation phase despite BSIPs uncertainties for our LLRR.
For the estimation phase, we used two exercises. The first simulated exercise was the squat, based on the requirements defined by the physiotherapists, see Section 2. The second exercise was a leg press. This exercise presents similar movements to a squat exercise, but in a supine position, i.e., with the subject facing up and their lower limbs elevated. These two exercises were simulated, and a physiotherapist validated the simulation trajectories.
Figure 5 presents a simulation of the squat exercise sequence, using an standing initial position with q = ( 90 , 0 , 90 ) and B r = 90 . Subsequently, the subject bends the knees to descend the body toward the floor with the heels on the floor, and the upper body remains aligned in the vertical plane [45]. It is important to note that our LLRR cannot let the subject descend their body toward the floor with their heels on the floor. Therefore, we simulated an “inverted” squat exercise, i.e., the subject does not lower their body towards the ground, but their feet rise towards their chest. Then, finally, the subject extends their knees to return to the standing position.
Figure 6 presents the q d for the hip, knee, and ankle. These movements are equivalent to the squat exercise. Moreover, we ensured that q ˙ d was not more significant than the one a human joint can execute. In addition, we can observe τ d y n , τ i n t , τ f r i , and τ a during the estimation phase. It is important to note the τ i n t had constant value since we simulated an isotonic exercise. During this simulation, the subject was relaxed during the first 5 s. Subsequently, we simulated a constant torque of τ i n t = ( 9.8 , 9.8 , 0 ) N   m , which is equivalent to the requirements presented in Section 2 suggested by the physiotherapists.
The simulation evaluated the response to sensor noise, i.e., the simulation used white Gaussian noise (WGN) with a signal-to-noise ratio SNR = 40 dB and 20% BSIPs uncertainties. We selected a WGN with SNR = 40 dB since it is the one used in the paper reported by Riani et al. [36].
Figure 7 reports the results for the simulation during the squat exercise. This simulation presents the outcome of the four estimation methods with SNR = 40 dB and 20% BSIPs uncertainties. Figure 7a compares the results for the IID and the Saadatzi et al. ID method. In this figure, we can observe that the IID works when there are BSIPs uncertainties. As expected, sensor noise affects both methods due to the computation of double derivatives to calculate velocity and acceleration. Therefore, our IID extends the technique reported in the literature by Saadatzi et al. [12]. Figure 7b contrasts the simulations for the INDO and the Saadatzi et al. NDO method. In this figure, we can denote that the INDO notably reduces the sensitivity to BSIPs uncertainties. Therefore, our INDO extends the technique reported in the literature by Saadatzi et al. [12]. Both NDO-based methods have significantly lower errors than the ID-based method. Moreover, they have a lower sensitivity to sensor noise due to the non-dependence of double derivatives. Finally, the IID and INDO work for a simulated LLRR despite the BSIPs uncertainties.
Table 3 reports the performance indices for the four human–robot interaction torque estimation methods during a squat exercise. In this table, we can denote that the INDO excels over the other three approaches. We can see that the INDO obtains the lowest Mean Absolute Error (MAE) and Mean Absolute Percentage Error (MAPE) for the hip and knee joints. Moreover, it has the lowest Root Mean Squared Error (RMS) and Root Mean Squared Percentage Error (RMSPE) for the same joints and the second lower value for all four indices in the ankle joint. Moreover, it obtains a value closer to one for the coefficient of correlation ( R 2 ) index, i.e., the desired value. Therefore, the overall performance of the INDO is superior in comparison to the other three methods in cases with BSIPs uncertainties.
Figure 8 depicts the leg press exercise sequence, which presents similar movements to a squat exercise but in a supine position, i.e., with the subject facing up and their lower limbs elevated. The leg press exercise was simulated by using an initial position with q = ( 40 , 0 , 90 ) and B r = 30 . Subsequently, the subject bends the knees to approach their chest, i.e., the subject begins the downward movement by slowly bending their hips and knees until their legs form a 90° angle. Finally, the subject extends their knee.
Figure 9 presents the q d for the hip, knee, and ankle. These movements are equivalent to the leg press exercise. Moreover, we can ensure that q ˙ d was not more significant than the one a human joint can execute. In addition, we can denote τ d y n , τ i n t , τ f r i , and τ a during the estimation phase. It is essential to note that the τ i n t was simulated with a constant value since we have used an isotonic exercise. The subject was relaxed during the first 5 s and a constant torque of τ i n t = ( 9.8 , 9.8 , 0 ) N m was executed by the subject from 5 s to 20 s. These torque magnitudes are equivalent to the requirements presented in Section 2 suggested by the physiotherapists. Finally, at 20 s, the subject stopped performing the interaction torque.
Figure 10 reports the results for the simulation during a leg press exercise. This simulation presents the results of the four estimation methods with SNR = 40 dB and 20% BSIPs uncertainties. Figure 10a compares the results for the IID and the Saadatzi et al. [12] ID method. The IID has significantly lower errors. The proposed method can accurately estimate the human–robot interaction torque despite BSIPs uncertainties. Both methods have a similar sensitivity to sensor noise. Figure 10b contrasts the simulations for the INDO and the Saadatzi et al. [12] NDO method. This figure depicts that the result obtained by the proposed method is superior to the approaches because it extends it by reducing the sensitivity to BSIPSs uncertainties. Moreover, both methods present an estimation with the lowest noise since their non-dependence on double derivatives. As a result, the errors obtained with the IID and INDO are significantly lower. Therefore, we can state that our proposed methods stand out for a proper torque estimation of the human–robot interaction torque even with sensor noise and BSIPs uncertainties. Thus, extending the work reported by Saadatzi et al. [12].
Table 4 reports the performance indices for the four human–robot interaction torque estimation methods during a leg press exercise. In this table, we can observe that the INDO excels over the other three approaches. We can see that the INDO obtains the lowest MAE, MAPE, RMS, and RMSPE for all three joints. Moreover, it obtained the closest value to one for the R 2 index, i.e., the desired value. Therefore, the INDO’s general performance outweighs in comparison of the other three methods.
As mentioned by Saadatzi et al. [12], their ID method is suitable for estimating the interaction torque when there is no noise in the sensors and there are no uncertainties in the parameters. It is expected to be noticeably affected by noise since it uses double derivatives to calculate velocity and acceleration. Additionally, it does not have a correction term like the NDO and only depends on the system’s dynamics. In contrast, in Figure 7a and Figure 10a, it can be seen that the IID notably reduces the sensitivity to BSIPs uncertainties. Therefore, the estimation with the IID outperforms when there are BSIPs uncertainties. However, the IID is affected by noisy sensors since the dependence of the doubles derivatives is not removed.
As demonstrated by Saadatzi et al. [12], their NDO method is good at estimating the interaction torque when there are no uncertainties in the parameters. Moreover, it is less sensitive to noise in the sensors than the ID. The Saadatzi et al. NDO method is a filtered version of the Saadatzi et al. ID method due to the nature of the observer and its non-dependence on double derivatives, only the first-order one. However, the Saadatzi et al. NDO method is just as sensitive to parameter uncertainties as the Saadatzi et al. ID method [12]. The INDO method notably reduces the sensitivity to BSIPs uncertainties, see Figure 7b and Figure 10b. In addition, the INDO inherits the lower sensitivity to noise from the Saadatzi et al. NDO method. Hence, an adequate estimation of the interaction torque is obtained with the INDO. However, it does not entirely reduce the noise effect to zero. Moreover, it may add a delay to the estimates. Defining the gain α of the NDO changes this sensitivity to noise and delay. The design for our INDO was made with the LMI, but an analytical and optimal design can be evaluated in further studies.
The proposed methods have no negative or positive effect on the sensitivity to sensors noise. Moreover, we can observe that the ID-based methods have more significant errors than the NDO-based methods. To summarize, we can state that, in this case, ID-based methods have more sensitivity to sensor noise than the NDO-based methods.
To sum up, we can see that the errors for the IID are significantly lower than the errors obtained with the Saadatzi et al. ID method under BSIPs uncertainties. Moreover, we can denote that the errors for the INDO are also significantly lower than the errors obtained with the Saadatzi et al. NDO method under BSIPs uncertainties. This significant error reduction with both proposed methods under BSIPs uncertainties, compared to the methods reported by Saadatzi et al. [12], is the most significant contribution of this paper. Therefore, we have proposed two human–robot interaction torque estimation observers for a simulated LLRR that works despite the BSIPs uncertainties. This proposal extends the work reported by Saadatzi et al. [12].

6.5. Requirements Evaluation

Table 5 presents an evaluation of the performance of the two proposed human–robot interaction torque estimation methods according to the requirements defined in Section 2. These requirements were defined to evaluate the performance of the human–robot torque estimation methods, to be used as biofeedback to favor isotonicity, despite the BSIPs uncertainties, through a simulation. In this table, we can observe that both proposed methods meet the vast majority of the requirements, i.e., accomplishing 16 (88.9%) and 13 (72.2%) of them. In contrast, both methods reported in the literature only accomplish 27.8% and 33.3% of the requirements when using simulations under noise and BSIPs uncertainties. Therefore, our proposed methods may be employed as a human–robot interaction torque estimation method for a simulated lower limb rehabilitation robot. Moreover, it is usable for isotonic exercises, as biofeedback, despite the BSIPs uncertainties. In addition, the INDO is superior to the ID-based method since it reduces the sensitivity to sensor noise. This result is reasonable because of the nature of the observer and because the NDO works without double derivates to compute q ¨ .
It is noticeable that the proposed IID and INDO behaved appropriately in the area where the methods were working. The good behavior of both methods was intensively verified using stimuli and simulations in the region in which the LLRR works.
Let us perform a stability analysis of suggested observers, considering that the IID and INDO are methods comprised of two phases.
Both methods have the same phase one, i.e., a calibration phase to reduce BSIPs uncertainties. The calibration phase for the IID and INDO inherits the stability from the online parameter estimation reported by Riani et al. [36]. We have shown that the calibration phase complies with the input trajectory design and the LMI design, e.g., the Genetic Algorithm optimization for the condition number ensured a persistent excitation trajectory a gave us an adequate input for the online parameter estimation. Thus, we obtained a trajectory whose points gave a well-conditioned observation matrix. Therefore, the calibration phase can use the same trajectory regardless of the patient’s BSIPs, and ensure that this estimator can converge during the calibration phase.
The second phase for the IID is based on a transformed ID. This method and transformation turn into a dynamical equation rather than an observer. Moreover, the stability properties are invariant under base changes. Therefore, this phase does not require a stability analysis.
The second phase for the INDO focuses on the convergence of the transformed NDO, i.e., the methods are used one after the other, and no feedback is used within them. Convergence was shown for the simulations, evidencing that it is very likely to work in the region of operation. This condition was preserved by ensuring the requirements for the transformed NDO design by using the LMI stated in Theorems 1 and 2.
To sum up, the IID and INDO stability is based on the fact that the convergence of each estimator separately is ensured. Moreover, they are executed one after the other. Therefore, the stability of each phase is maintained.
This work’s contribution consists of two methods that do not require the physiotherapist to make an exact measurement of the BSIPs of the patient’s limbs, but allow the use of approximate values. This contribution is possible since we based our methods on separating the parameter identification phase from the torque estimation phase and both from the controller itself. To do this, the physiotherapist may ask the patient to be in a state of rest, that is, without exerting an intentional torque. The identification system would take the data from the first iterations to adjust those measurements given by the physiotherapist, which were an initial value. It would adapt them to a functional value without having to measure them. Therefore, it would provide the LLRR with more “intelligence” as it would tune itself. Thus, it allows that in a second stage, once identifying the robot–human system, a known method is applied for torque estimation, which traditionally has the problem of being very dependent on an accurate model and therefore sensitive to BSIPs uncertainties. Still, we proposed two methods that reduce this uncertainty. It is essential to mention that the torque estimation is not necessary for the execution of the controller because, in this type of rehabilitation, it is the patient who imposes the torque on the robot and not the other way around.

7. Conclusions

We proposed the Identified ID-based method (IID) and the Identified NDO-based method (INDO). Both methods estimate the human–robot interaction torque for a simulated Lower Limb Rehabilitation Robot (LLRR), which is usable for isotonic exercises despite the uncertainties of the Body Segment Inertial Parameter (BSIPs). Both methods have two phases. (1) Phase one is the calibration phase, which allows an online parameter estimation to reduce BSIPs uncertainties. (2) Phase two is the estimation phase, where the human–robot interaction torque uses the estimated parameters under the assumption that the BSIPs uncertainties have been reduced. The first method combines a transformed Inverse Dynamics (ID), a friction model, and an online parameter estimation method. The second method uses the conjunction of a transformed Nonlinear Disturbance Observer (NDO) and a friction model fed by the output of an online parameter estimation method. Two sets of simulations were conducted using white Gaussian noise (WGN) with a signal-to-noise ratio SNR = 40 dB and 20% BSIPs uncertainties. The two proposed methods’ performance was compared to two of the best methods reported in the literature. Both proposed methods obtained the best Coefficient of Correlation, Mean Absolute Error, and Root Mean Squared Error compared to the benchmarks. Moreover, the IID and INDO met the vast majority of the requirements, i.e., accomplishing 72.2% and 88.9%, respectively. In contrast, both methods reported in the literature only accomplish 27.8% and 33.3% of the requirements when using simulations with under noise and BSIPs uncertainties. In addition, the INDO is superior to the IID since it reduces the sensitivity to sensor noise. In conclusion, this paper proposes two human–robot interaction torque estimation methods, which extend two methods reported in the literature and copes with BSIPs uncertainties without using additional sensors.
Future work includes a practical realization of proposed observers within a real environment with the actual LLRR Nukawa and conducting a study with healthy subjects and patients. This work requires additional approval from the ethics committee. Moreover, the validation would ideally involve hardware using a torque sensor on each joint to measure the actual interaction torque experimentally. Thus, we would be able to compare this measure with the torque estimation made via the two proposed methods. A challenge for implementation in the physical robot is the existence of static friction, which the methods do not consider. Therefore, the estimated torques may be different from the real ones during the start and stop of the robot.

Author Contributions

Conceptualization, J.C.Y., S.R., M.O., V.Z.P., J.A.M., A.A.-J. and M.J.B.; methodology, J.C.Y., S.R., M.O., V.Z.P., J.A.M., A.A.-J. and M.J.B.; validation, J.C.Y., S.R., M.O., V.Z.P., J.A.M., A.A.-J. and M.J.B.; formal analysis, J.C.Y., S.R., M.O., V.Z.P., J.A.M. and M.J.B.; investigation, J.C.Y., S.R., M.O., V.Z.P., J.A.M., A.A.-J. and M.J.B.; writing—original draft preparation, J.C.Y. and S.R.; writing—review and editing, J.C.Y., S.R., M.O., V.Z.P., J.A.M., A.A.-J. and M.J.B.; software, J.C.Y., S.R., M.O., V.Z.P., J.A.M. and M.J.B.; funding acquisition, J.C.Y., S.R., V.Z.P., J.A.M. and M.J.B. All authors have read and agreed to the published version of the manuscript.

Funding

This paper was framed under the project “KINA: Virtual Reality System for Lower Limb Rehabilitation of APM or IED Victims” funded by the call 808 of 2018 for projects from the Ministerio de Ciencia, Tecnologia e Innovacion de Colombia (Minciencias), project 121080864383. Moreover, the Ph.D. studies of Juan C. Yepes were funded by Minciencias and the Universidad Pontificia Bolivariana (UPB) under call 727 of 2015. In addition, the Australian Academy of Science (AAS) and Minciencias granted the funds for his Ph.D. internship within the Australia-Americas Ph.D. Research Internships Program 2019. Additionally, this work was supported by UNAM-PAPIIT IN 102121. Finally, this work was supported by Universidad de Medellín (UdeM).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgments

First and foremost, we have to thank all our colleagues for their support. It has been a challenging time due to the COVID-19 pandemic. We also thank Minciencias and the AAS for their grants. It would not have been possible to carry out this research and Ph.D. internship without their help. Finally, we thank the UNAM and UdeM for their support.

Conflicts of Interest

The authors declare that they have no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
AASAustralian Academy of Science
AIArtifficial Intelligence
APMAntipersonnel mines
BSIPsBody Segment Inertial Parameters
CADComputer Assisted Design
CoMCenter of Mass
COVID-19Coronavirus disease 2019
CTCComputed torque control
D-HDenavit Hartenberg
DOFDegrees of freedom
DPDorsi/plantar
ESMADMobile Anti-Disturbances Squadron
FEFlexion/extension
GAGenetic Algorithm
GRFGround Reaction Force
IEDImprovised explosive devices
IDInverse Dynamics
IIDIdentified Inverse Dynamics
INDOIdentified Nonlinear Disturbance Observer
KEKinetic Energy
KINAVirtual Reality System for Lower Limb Rehabilitation of APM or IED Victims
LegSysOld name for the lower limb rehabilitation robot
LLRRLower Limb Rehabilitation
LMILinear Matrix Inequality
MAEMean Absolute Error
MAPEMean Absolute Percentage Error
MincienciasMinisterio de Ciencias, Tecnología e Innovación de Colombia
NDONonlinear Disturbance Observer
NukawaCurrent name for the lower limb rehabilitation robot system
PEPotential Energy
R 2 Coefficient of Correlation
ROMRange of Motion
RMSERoot Mean Square Error
RMSPERoot Mean Square Percentage Error
sEMGSurface Electromyography
SDGSSustainable Development Goals
SNRSignal-to-noise Ratio
UdeMUniversidad de Medellín
UPBUniversidad Pontificia Bolivariana
UNAM        Universidad Nacional Autónoma de México
UXOUnexploded Ordnance
WGNWhite Gaussian Noise
WHOWorld Health Organization

Appendix A. Nukawa

We have developed a lower limb rehabilitation robot (LLRR) called “Nukawa”. This system has its antecedents in the LegSys system [42,46,47,48]. Figure A1 presents the current version of the LLRR Nukawa. The robotic system Nukawa is a product of technical requirements proposed by an interdisciplinary group formed by physiotherapists, industrial designers, and engineers. The design consists of two limbs, each one composed of a three-link mechanism and an electronic position and force control, i.e., each leg has 3 degrees of freedom (3-DOF). The design also has brushless DC motors, power drivers, and position sensors to perform a control strategy capable of generating multiple rehabilitation patterns.
The system would perform flexion/extension (FE) movements of the hip, FE movements of the knee, and dorsi/plantar (DP) flexion movements of the ankle [46]. Additionally, the joints are collinear to human joints. The knee is a polycentric joint. However, a collinear simplification was conducted as presented by Zoss et al. [49] involves a pure rotational joint in the sagittal plane for the knee joint. To adjust the system for each person, the length of each segment of the robotic system Nukawa is adjustable, i.e., the length between each joint can be adjusted. The system is designed for people between 1.44   m and 1.85   m tall, adjustable using a couple of telescopic mechanical systems. Besides, the robot was designed for people up to 85 k g weight. The brushless DC motors were selected to compensate for the weight of the robot’s links and the subject limbs. When needed, the motors may also generate a resistive opposition hip and knee FE and ankle DP. Nukawa aims to help the physiotherapist in the process of rehabilitation of musculoskeletal pathologies of the lower limb, e.g., subjects victims of APM, IED, or UXO.
Figure A1. Nukawa, a 3-DOF Lower-limb robotic rehabilitation for hip, knee, and ankle joints. (a) CAD model and (b) current version of the robot with both limbs.
Figure A1. Nukawa, a 3-DOF Lower-limb robotic rehabilitation for hip, knee, and ankle joints. (a) CAD model and (b) current version of the robot with both limbs.
Applsci 12 05529 g0a1

Appendix B. Physiotherapist Surveys

The idea was to ask people in the health area about the requirements of a human–robot interaction torque estimation method for an LLRR usable for isotonic exercises. Moreover, the survey purpose was to define the movements for the simulations, i.e., the isotonic exercises to be used within the simulations. We designed the surveys involving the terms that physiotherapists commonly use during their evaluation and rehabilitation of patients with lower limb injuries due to APM, IED, and UXO. Therefore, the questionnaire was created by a physiotherapist. The expert had experience in research, clinical rehabilitation, and therapeutic services.
It is worth noting that in the framework of this survey, we took the definition of isotonic exercises reported by Lopez et al. [50] and this terminology was explained to the physiotherapists before the survey. Let us introduce this definition as:
An Isotonic contraction is an exercise producing constant tension. This term is generally applied when the external resistance is constant, achieved with a variable resistance machine. It is essential to bear in mind that in normal muscle movement in humans, there are no muscle contractions in which the force remains the same throughout a workout. Moreover, the tension generated in the muscle will change as the lever arms change. Therefore, a machine with variable resistance is required throughout the ROM [50].
Table A1 presents all the questions included in the survey. It consisted of two multiple-choice multiple-answer questions, two multiple-choice single-answer questions, and one free-response question. We contacted the physiotherapists by virtual means, and they filled out the form in the cloud. A total of 26 physiotherapists answered the anonymous questionnaire. Figure A2 reports the results of the four multiple-choice questions. We asked these four questions to define the simulations’ ROM, speed, force, and maximum tolerable error.
Figure A2a depicts the results of asking them about the usage of isotonic exercises according to ROM recovery, i.e., the question focused on determining if these exercises are used before recovering the ROM, after recovering the ROM, or in both cases. In this pie chart, we can observe that the tendency is greater toward using isotonic exercises in all phases of rehabilitation, with a total of 17 (65.4%) subjects who agree. Therefore, we require to let the user configure the angle limits of the joints, i.e., the ROM to be used within the simulations.
Table A1. Questions included in the survey.
Table A1. Questions included in the survey.
Questions
1. How are performed Isotonic exercises for the rehabilitation of antipersonnel mines (APM), improvised explosive devices (IED), and unexploded ordnance (UXO)?
2. List, in order of importance, five isotonic exercises necessary for lower limb rehabilitation in people victims of APM, IED, and UXO.
3. Which segmental speeds are used in isotonic exercises to rehabilitate lower limbs in people victims of APM, IED, and UXO?
4. In a rehabilitation process for APM, IED, and UXO victims, what are the ranges of force, in pounds or kilograms, during the exercises?
5. What would be the maximum percentage of error allowed in a device for the automatic estimation of the force performed by the subject during the execution of a motor activity?
6. Optional question: considering the following case study of an amputee, write in the table an example of a protocol that you would use in the rehabilitation of the IED victim, including five (5) isotonic exercises for lower limb rehabilitation. Define the ranges of motion used in each exercise, the speeds (reps per minute), and the force range (in kg or pounds) performed at each joint. Case study of a person amputated by IED: male patient, 30 years old, weighing 75 k g , 175 m tall; with right transfemoral amputation (the distal third of the knee), a stump with 18 c m length from the perineum to the femur section. A user without vascular problems and with good soft tissue healing. The amputation was caused by exposure to an improvised explosive device during his activities as part of the ESMAD (Command of special operational units) during the control of disturbances and blockades in the rural area of Florencia Caquetá. As a consequence of the blast wave, the patient lost the hearing capacity of the right ear (sensory hearing loss), additionally to the amputation. He is currently undergoing gait rehabilitation and uses an Ottobock 3R80 hydraulic knee prosthesis as a device.
Figure A2b presents the results of the number of repetitions used within a set of isotonic exercises for the rehabilitation of subjects victims of APM, IED, and UXO. In this figure, we can observe that 21 (80,7%) subjects answered that a total of 1–25 repetitions per set should be conducted. Therefore, we can indirectly have a notion of the speed of the movements. This result implies that slow exercises are required for the simulations.
Figure A2c presents the results of the ranges of force during the isotonic exercises for lower limb rehabilitation of subjects victims of APM, IED, and UXO. It is important to note that physiotherapists usually express force and torque in kilograms or pounds. Therefore, we used these terms during the questionnaire, as suggested by the expert in rehabilitation. In this figure, we can observe that the physiotherapists tend to use the lower force ranges for the hip with 14 (53.8%) subjects, lower force ranges for the knee with 13 (50%) subjects and lower force ranges for the ankle with 16 (61.5%) subjects. According to the evidence in these results, it is concluded that the human–robot interaction torque estimation algorithm is required to work with the following ranges: 0   k g to 15   k g for the hip, 0   k g to 15   k g for the knee, and 0   k g to 10   k g for the ankle.
Figure A2d presents the results of asking the physiotherapists about the maximum percentage of error allowed in a device for the automatic estimation of the force carried out by the subject during the execution of the motor activity. In this figure, we can observe that the predominant response is that the maximum percentage of error allowed is in the range from 1% to 3%. Therefore, the requirement from a physiotherapeutic point of view is that the human–robot interaction torque estimation method should have a percentage error within this range.
In addition, to define the movements conducted by the LLRR during the simulations, we presented a free-response question to the physiotherapists. These movements were used to evaluate the performance of the human–robot interaction torque estimation method. The free-response question was optional. However, we mentioned that it is the one that contributes to our study most. For this question, we presented a clinical case to the physiotherapists. We asked them to write an example of a protocol that they would use to rehabilitate the person victim of IED, including five isotonic exercises for rehabilitation of lower limbs. We also asked them to define the ROM, speed, and ranges of force performed by each joint during each isotonic exercise. The case study of a subject with an amputation due to an IED is presented in Table A1.
The physiotherapists suggested a total of 44 exercises. However, 25 of them did not meet the exclusion criteria:
  • Exercises that could not be executed by our LLRR.
  • Exercises with ROM, speed, or forces without quantitative quantities.
  • Repeated exercises.
  • Exercises with incomplete information.
A total of 19 isotonic exercises were obtained after evaluating the exclusion criteria. Therefore, they were considered according to the inclusion criteria: (a) Exercises that involve the movement of the three joints so that it is not trivial, and (b) symmetric exercises, i.e., executing the same movement for each leg as a mirror. Finally, we obtained one exercise that met all exclusion criteria and inclusion criteria, named squat:
  • ROM—Hip 90° and Knee 90°;
  • Speed—20 repetitions in 4 min;
  • Force—2 lb ( 1   k g ) at hip, 2 lb ( 1   k g ) at knee, 0 lb ( 0   k g ) at ankle.
Therefore, the squat exercise should be used to evaluate the performance of the human–robot interaction torque estimation method. Moreover, the leg press exercise was also suggested by a physiotherapist as a second exercise. This exercise presents similar movements to a squat exercise, but in a supine position, i.e., with the subject facing up and their lower limbs elevated. These two exercises were simulated, and a physiotherapist validated the simulation trajectories.
Figure A2. Results of the survey reporting that isotonic exercises should be (a) used in these rehabilitation phases, (b) executed taking into account these speeds for the movements, (c) executed using these force ranges, and (d) executed with a human-robot interaction torque estimation with these maximum errors.
Figure A2. Results of the survey reporting that isotonic exercises should be (a) used in these rehabilitation phases, (b) executed taking into account these speeds for the movements, (c) executed using these force ranges, and (d) executed with a human-robot interaction torque estimation with these maximum errors.
Applsci 12 05529 g0a2

Appendix C. Dynamic Model

Defining s · = sin ( · ) , c · = cos ( · ) , c i = c q i , s i = s q i , we obtained the matrix
M ( q ) = m 11 ( q ) m 12 ( q ) m 13 ( q ) m 21 ( q ) m 22 ( q ) m 23 ( q ) m 31 ( q ) m 32 ( q ) m 33 ( q )
with
m 11 ( q ) = I c 1 , z z + I c 2 , z z + I c 3 , z z + 1 2 m 2 + 1 2 m 3 + 2 2 m 3 + b 1 2 m 1 + b 2 2 m 2 + b 3 2 m 3 + 2 1 b 3 m 3 c 23 + 2 1 2 m 3 c 2 + 2 1 b 2 m 2 c 2 + 2 2 b 3 m 3 c 3 ,
m 22 ( q ) = m 3 2 2 + 2 m 3 c 3 2 b 3 + m 2 b 2 2 + m 3 b 3 2 + I c 2 , z z + I c 3 , z z ,
m 33 ( q ) = m 3 b 3 2 + I c 3 , z z ,
m 12 ( q ) = m 21 ( q ) = m 3 2 2 + 2 m 3 c 3 2 b 3 + 1 m 3 c 2 2 + m 2 b 2 2 + 1 m 2 c 2 b 2 + m 3 b 3 2 + 1 m 3 c 23 b 3 + I c 2 , z z + I c 3 , z z ,
m 13 ( q ) = m 31 ( q ) = I c 3 , z z + b 3 2 m 3 + 1 b 3 m 3 c 23 + 2 b 3 m 3 c 3 ,
and
m 23 ( q ) = m 32 ( q ) = m 3 b 3 2 + 2 m 3 c 3 b 3 + I c 3 , z z
Therefore, the C ( q , q ˙ ) terms is
C ( q , q ˙ ) = C 1 ( q , q ˙ ) C 2 ( q , q ˙ ) C 3 ( q , q ˙ )
with
C 1 ( q , q ˙ ) = 1 2 m 3 q ˙ 2 2 s 2 1 b 2 m 2 q ˙ 2 2 s 2 2 b 3 m 3 q ˙ 3 2 s 3 1 b 3 m 3 q ˙ 2 2 c 2 s 3 1 b 3 m 3 q ˙ 2 2 c 3 s 2 1 b 3 m 3 q ˙ 3 2 c 2 s 3 1 b 3 m 3 q ˙ 3 2 c 3 s 2 2 1 2 m 3 q ˙ 1 q ˙ 2 s 2 2 1 b 2 m 2 q ˙ 1 q ˙ 2 s 2 2 2 b 3 m 3 q ˙ 1 q ˙ 3 s 3 2 2 b 3 m 3 q ˙ 2 q ˙ 3 s 3 2 1 b 3 m 3 q ˙ 1 q ˙ 2 c 2 s 3 2 1 b 3 m 3 q ˙ 1 q ˙ 2 c 3 s 2 2 1 b 3 m 3 q ˙ 1 q ˙ 3 c 2 s 3 2 1 b 3 m 3 q ˙ 1 q ˙ 3 c 3 s 2 2 1 b 3 m 3 q ˙ 2 q ˙ 3 c 2 s 3 2 1 b 3 m 3 q ˙ 2 q ˙ 3 c 3 s 2
C 2 ( q , q ˙ ) = 1 2 m 3 q ˙ 1 2 s 2 + 1 b 2 m 2 q ˙ 1 2 s 2 2 b 3 m 3 q ˙ 3 2 s 3 + 1 b 3 m 3 q ˙ 1 2 c 2 s 3 + 1 b 3 m 3 q ˙ 1 2 c 3 s 2 2 2 b 3 m 3 q ˙ 1 q ˙ 3 s 3 2 2 b 3 m 3 q ˙ 2 q ˙ 3 s 3
and
C 3 ( q , q ˙ ) = q ˙ 1 q ˙ 1 1 b 3 m 3 s 23 + 2 b 3 m 3 s 3 + 2 b 3 m 3 q ˙ 2 s 3 + q ˙ 2 2 b 3 m 3 q ˙ 1 s 3 + 2 b 3 m 3 q ˙ 2 s 3
Now, taking into account that the gravity term is defined as
G ( q ) = U ( q ) q T = g 1 ( q ) g 2 ( q ) g 3 ( q )
we obtained
g 1 ( q ) = g 0 m 3 2 c 12 + 1 c 1 + b 3 c 123 + g 0 m 2 b 2 c 12 + 1 c 1 + b 1 g 0 m 1 c 1
g 2 ( q ) = g 0 m 3 2 c 12 + b 3 c 123 + b 2 g 0 m 2 c 12
g 3 ( q ) = b 3 g 0 m 3 c 123
The friction may be expressed as
τ f r i = F v q ˙
where F v is the coefficient matrix of viscous friction, which is proportional to the velocity of joint motion.

Appendix D. Transformed Model

The introduction of dynamic coefficients χ i is a convenient regrouping of the dynamic parameters, i.e., a linear parametrization of dynamics. The robot dynamics depend in a nonlinear way on some of these parameters. Therefore, we assumed
χ 1 = I c 1 , z z + I c 2 , z z + I c 3 , z z + L 1 2 m 2 + L 1 2 m 3 + L 2 2 m 3 + b 1 2 m 1 + b 2 2 m 2 + b 3 2 m 3 χ 2 = L 1 g 0 m 2 + L 1 g 0 m 3 + b 1 g 0 m 1 χ 3 = m 3 L 2 2 + m 2 b 2 2 + m 3 b 3 2 + I c 2 , z z + I c 3 , z z χ 4 = L 1 L 2 m 3 + L 1 b 2 m 2 χ 5 = L 2 g 0 m 3 + b 2 g 0 m 2 χ 6 = m 3 b 3 2 + I c 3 , z z χ 7 = L 2 b 3 m 3 χ 8 = L 1 b 3 m 3 χ 9 = b 3 g 0 m 3
Therefore, we obtained a transformed M ( χ , q ) defined as
M ( χ , q ) = m 11 ( χ , q ) m 12 ( χ , q ) m 13 ( χ , q ) m 21 ( χ , q ) m 22 ( χ , q ) m 23 ( χ , q ) m 31 ( χ , q ) m 32 ( χ , q ) m 33 ( χ , q )
with
m 11 ( χ , q ) = χ 1 + 2 χ 4 c 2 + 2 χ 7 c 3 + 2 χ 8 c 23
m 22 ( χ , q ) = χ 3 + 2 χ 7 c 3
m 33 ( χ , q ) = χ 6
m 12 ( χ , q ) = m 21 ( χ , q ) = χ 3 + χ 4 c 2 + 2 χ 7 c 3 + χ 8 c 23
m 13 ( χ , q ) = m 31 ( χ , q ) = χ 6 + χ 7 c 3 + χ 8 c 23
m 23 ( χ , q ) = m 32 ( χ , q ) = χ 6 + χ 7 c 3
The Lagrangian L ( q , q ˙ ) of the mechanical system can be defined as
L ( q , q ˙ ) = T ( q , q ˙ ) U ( q ) = 1 2 i , j m i j ( q ) q ˙ i q ˙ j U ( q )
where T and U denote the total KE and PE of the system, respectively.
The Euler–Lagrange equations can be defined as
d d t L q ˙ k L q k = u k k = 1 , 2 , 3
where u k is the generalized force associated with the generalized coordinate q i . Therefore, we rewrite the dynamic equations as
j m k j ( q ) q ¨ j Inertial terms + i , j c k i j ( q ) q ˙ i q ˙ j Centrifugal ( i = j ) and Coriolis ( i j ) terms + U q k Gravity terms g k ( q ) = u k
where c k i j = c k j i are the Christoffel symbols of the first kind. Subsequently, we computed the Centrifugal and Coriolis term c ( q , q ˙ ) taking into account that the kth component of vector C is defined as
c k ( , q ˙ ) = q ˙ T C k ( q ) q ˙
where
C k ( q ) = 1 2 M k q + M k q T M q k
and M k is the k-th column of matrix M ( q ) . Subsequently, we obtained a transformed C ( χ , q , q ˙ ) as
C ( χ , q , q ˙ ) = C 1 ( χ , q , q ˙ ) C 2 ( χ , q , q ˙ ) C 3 ( χ , q , q ˙ )
with
C 1 ( χ , q , q ˙ ) = q ˙ 2 q ˙ 1 χ 8 s 23 + χ 4 s 2 + q ˙ 2 χ 8 s 23 + χ 4 s 2 + q ˙ 3 χ 8 s 23 + χ 7 s 3 q ˙ 3 q ˙ 1 χ 8 s 23 + χ 7 s 3 + q ˙ 2 χ 8 s 23 + χ 7 s 3 + q ˙ 3 χ 8 s 23 + χ 7 s 3 q ˙ 1 q ˙ 2 χ 8 s 23 + χ 4 s 2 + q ˙ 3 χ 8 s 23 + χ 7 s 3
C 2 ( χ , q , q ˙ ) = q ˙ 1 q ˙ 1 χ 8 s 23 + χ 4 s 2 χ 7 q ˙ 3 s 3 q ˙ 3 χ 7 q ˙ 1 s 3 + χ 7 q ˙ 2 s 3 + χ 7 q ˙ 3 s 3 χ 7 q ˙ 2 q ˙ 3 s 3
C 3 ( χ , q , q ˙ ) = q ˙ 2 χ 7 q ˙ 1 s 3 + χ 7 q ˙ 2 s 3 + q ˙ 1 q ˙ 1 χ 8 s 23 + χ 7 s 3 + χ 7 q ˙ 2 s 3
Now, using the set of dynamic parameters (A17) and taking into account that G ( χ , q ) is defined as
G ( χ , q ) = g 1 ( χ , q ) g 2 ( χ , q ) g 3 ( χ , q )
we obtain the elements for G ( χ , q ) as
g 1 ( χ , q ) = χ 2 c 1 + χ 5 c 12 + χ 9 c 123
g 2 ( χ , q ) = χ 5 c 12 + χ 9 c 123
g 3 ( χ , q ) = χ 9 c 123
We combined these standard parameters to get the vector χ of the base inertial parameters as given in [38]. This linear parametrization of robot dynamics is not unique in the chosen set of dynamic coefficients χ and the associated regression matrix W . The W is a n × b matrix, where n = 3 is the number of DOF of the LLRR and b = 9 is the number of the dynamic parameters. The order of these standard parameters was selected to guarantee that the regressor W ( q , q ˙ , q ¨ ) has a block upper triangular structure as reported by [30]. This triangular block structure may enable computing the parameter estimation by using a sequential procedure. In this case, the parameter estimation can be executed by an iterative method, taking measurements from the distal limb to the proximal. However, such a technique may have the drawback of accumulating any error due to ill-conditioning of the matrices involved step-by-step. It may then be worthwhile to use a global approach that imposes motions on all LLRR joints at once [30]. The calibration phase presented in Section 5.2 will show a multi-DOF approach for an online estimation to reduce BSIPS uncertainties.
χ = χ 1 χ 2 χ 3 χ 4 χ 5 χ 6 χ 7 χ 8 χ 9 T
so the dynamic model (1) can be written as
τ = W ( q , q ˙ , q ¨ ) χ
where W ( q , q ˙ , q ¨ ) is the regressor matrix, and its notation may be simplified as W .
W ( q , q ˙ , q ¨ ) = w 11 w 12 w 13 w 14 w 15 w 16 w 17 w 18 w 19 0 0 w 23 w 24 w 25 w 26 w 27 w 28 w 29 0 0 0 0 0 w 36 w 37 w 38 w 39
and the elements of the W ( q , q ˙ , q ¨ ) are defined as
w 11 ( q , q ˙ , q ¨ ) = q ¨ 1 w 12 ( q , q ˙ , q ¨ ) = c 1 w 13 ( q , q ˙ , q ¨ ) = q ¨ 2 w 14 ( q , q ˙ , q ¨ ) = s 2 q ˙ 2 2 2 q ˙ 1 s 2 q ˙ 2 + 2 q ¨ 1 c 2 + q ¨ 2 c 2 w 15 ( q , q ˙ , q ¨ ) = c 12 w 16 ( q , q ˙ , q ¨ ) = q ¨ 3 w 17 ( q , q ˙ , q ¨ ) = 2 q ¨ 1 c 3 q ˙ 3 2 s 3 + 2 q ¨ 2 c 3 + q ¨ 3 c 3 2 q ˙ 1 q ˙ 3 s 3 2 q ˙ 2 q ˙ 3 s 3 w 18 ( q , q ˙ , q ¨ ) = s 23 q ˙ 2 2 2 s 23 q ˙ 2 q ˙ 3 2 q ˙ 1 s 23 q ˙ 2 s 23 q ˙ 3 2 2 q ˙ 1 s 23 q ˙ 3 + 2 q ¨ 1 c 23 + q ¨ 2 c 23 + q ¨ 3 c 23 w 19 ( q , q ˙ , q ¨ ) = c 123 w 23 ( q , q ˙ , q ¨ ) = q ¨ 1 + q ¨ 2 w 24 ( q , q ˙ , q ¨ ) = s 2 q ˙ 1 2 + q ¨ 1 c 2 w 25 ( q , q ˙ , q ¨ ) = c 12 w 26 ( q , q ˙ , q ¨ ) = q ¨ 3 w 27 ( q , q ˙ , q ¨ ) = 2 q ¨ 1 c 3 q ˙ 3 2 s 3 + 2 q ¨ 2 c 3 + q ¨ 3 c 3 2 q ˙ 1 q ˙ 3 s 3 2 q ˙ 2 q ˙ 3 s 3 w 28 ( q , q ˙ , q ¨ ) = s 23 q ˙ 1 2 + q ¨ 1 c 23 w 29 ( q , q ˙ , q ¨ ) = c 123 w 36 ( q , q ˙ , q ¨ ) = q ¨ 1 + q ¨ 2 + q ¨ 3 w 37 ( q , q ˙ , q ¨ ) = s 3 q ˙ 1 2 + 2 s 3 q ˙ 1 q ˙ 2 + s 3 q ˙ 2 2 + q ¨ 1 c 3 + q ¨ 2 c 3 w 38 ( q , q ˙ , q ¨ ) = s 23 q ˙ 1 2 + q ¨ 1 c 23 w 39 ( q , q ˙ , q ¨ ) = c 123

References

  1. World Health Organization. World Report on Disability; Technical Report; World Health Organization: Geneva, Switzerland, 2011. [Google Scholar]
  2. United Nations Development Programme. UNDP Support to the Implementation of Sustainable Development Goal 1: Poverty Reduction; Technical Report; United Nations Development Programme: New York, NY, USA, 2016. [Google Scholar]
  3. Departamento Administrativo Nacional de Estadística. Defunción por Grupo de Edad, Sexo, Segœn Departamento, Municipio y área de Residencia; Departamento Administrativo Nacional de Estadística: Bogotá, Colombia, 2014.
  4. Eurostat. Disability Statistics—Prevalence and Demographics; Technical Report; European Union: Maastricht, The Netherlands, 2014.
  5. Organización Panamericana de la Salud. Aplicacion de la Clasificación Internacional del Funcionamiento, de la Discapacidad y de la Salud en Estudios de Prevalencia de Discapacidad en las Américas; Technical Report; Organización Panamericana de la Salud: Washington, DC, USA, 2012. [Google Scholar]
  6. Almeras, D.; Milosavljevic, V. Informe regional sobre la medición de la discapacidad. Una mirada a los procedimientos de medición de la discapacidad en América Latina y el Caribe; Comisión Económica para América Latina y el Caribe: Santiago, Chile, 2014; Volume 3860, p. 48. [Google Scholar]
  7. Departamento Nacional de Planeación. Bases del Plan Nacional de Desarrollo 2018–2022; Technical Report; Departamento Nacional de Planeación: Bogotá, Colombia, 2019.
  8. Descontamina Colombia. Víctimas de Minas Antipersonal y Municiones sin Explosionar; Descontamina Colombia: Bogotá, Colombia, 2019.
  9. Dollar, A.M.; Herr, H. Lower extremity exoskeletons and active orthoses: Challenges and state-of-the-art. IEEE Trans. Robot. 2008, 24, 144–158. [Google Scholar] [CrossRef]
  10. Choi, J.H.; Kwak, J.H.; An, J.; Oh, S. Force Sensorless Multi-functional Impedance Control for Rehabilitation Robot. IFAC-PapersOnLine 2017, 50, 12077–12082. [Google Scholar] [CrossRef]
  11. Hwang, B.; Jeon, D. A method to accurately estimate the muscular torques of human wearing exoskeletons by Torque sensors. Sensors 2015, 15, 8337–8357. [Google Scholar] [CrossRef] [PubMed]
  12. Saadatzi, M.; Long, D.C.; Celik, O. Comparison of Human-Robot Interaction Torque Estimation Methods in a Wrist Rehabilitation Exoskeleton. J. Intell. Robot. Syst. 2018, 10, 1–17. [Google Scholar] [CrossRef]
  13. Kong, K.; Jeon, D. Design and control of an exoskeleton for the elderly and patients. IEEE/ASME Trans. Mechatron. 2006, 11, 428–432. [Google Scholar] [CrossRef]
  14. Wang, X.; Chen, F.; Zha, F.; Qiu, S.; Li, M.; Deng, J. Towards Online Estimation of Human Joint Muscular Torque with a Lower Limb Exoskeleton Robot. Appl. Sci. 2018, 8, 1610. [Google Scholar] [CrossRef] [Green Version]
  15. Grun, M.; Konigorski, U. Observer Based Method for Joint Torque Estimation in Active Orthoses. IFAC Proc. Vol. 2012, 45, 199–204. [Google Scholar] [CrossRef]
  16. Hayashi, T.; Kawamoto, H.; Sankai, Y. Control method of robot suit HAL working as operator’s muscle using biological and dynamical information. In Proceedings of the 2005 IEEE/RSJ International Conference on Intelligent Robots and Systems, Edmonton, AB, Canada, 2–6 August 2005; Volume 2, pp. 3455–3460. [Google Scholar] [CrossRef]
  17. Anwar, T.; Al-Dmour, H. RBF based adaptive neuro-fuzzy inference system to torque estimation from EMG signal. In Proceedings of the 2017 IEEE Symposium Series on Computational Intelligence (SSCI), Honolulu, HI, USA, 27 November–1 December 2017; pp. 1–8. [Google Scholar]
  18. Heine, C.B.; Menegaldo, L.L. Numerical validation of a subject-specific parameter identification approach of a quadriceps femoris EMG-driven model. Med Eng. Phys. 2018, 53, 66–74. [Google Scholar] [CrossRef]
  19. Zhang, L.; Zhu, X.; Gutierrez Farewik, E.; Wang, R. Estimation of ankle dynamic joint torque by a neuromusculoskeletal solver-informed NN model. In Proceedings of the 2021 6th IEEE International Conference on Advanced Robotics and Mechatronics, ICARM 2021, Chongqing, China, 3–5 July 2021; pp. 75–80. [Google Scholar] [CrossRef]
  20. Aljuboury, A.S.; Hameed, A.H.; Ajel, A.R.; Humaidi, A.J.; Alkhayyat, A.; Mhdawi, A.K.A. Robust Adaptive Control of Knee Exoskeleton-Assistant System Based on Nonlinear Disturbance Observer. Actuators 2022, 11, 78. [Google Scholar] [CrossRef]
  21. Saadatzi, M.; Long, D.C.; Celik, O. Torque estimation in a wrist rehabilitation robot using a nonlinear disturbance observer. In Proceedings of the ASME 2015 Dynamic Systems and Control Conference, Columbus, OH, USA, 28 October 2015. [Google Scholar]
  22. Kong, K.; Moon, H.; Hwang, B.; Jeon, D.; Tomizuka, M. Robotic rehabilitation treatments: Realization of aquatic therapy effects in exoskeleton systems. In Proceedings of the 2009 IEEE International Conference on Robotics and Automation, Kobe, Japan, 12–17 May 2009; Volume 1, pp. 1923–1928. [Google Scholar] [CrossRef]
  23. Brahmi, B.; Saad, M.; Ochoa-Luna, C.; Rahman, M.H. Adaptive control of an exoskeleton robot with uncertainties on kinematics and dynamics. In Proceedings of the IEEE International Conference on Rehabilitation Robotics, London, UK, 17–20 July 2017; pp. 1369–1374. [Google Scholar] [CrossRef]
  24. Fahmy, A.A.; Abdel Ghany, A.M. Neuro-fuzzy inverse model control structure of robotic manipulators utilized for physiotherapy applications. Ain Shams Eng. J. 2013, 4, 805–829. [Google Scholar] [CrossRef] [Green Version]
  25. Winter, D.A. Biomechanics and Motor Control of Human Movement, 2nd ed.; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2009; p. 277. [Google Scholar] [CrossRef]
  26. Plagenhoef, S.; Evans, F.G.; Abdelnour, T. Anatomical Data for Analyzing Human Motion. Res. Q. Exerc. Sport 1983, 54, 169–178. [Google Scholar] [CrossRef]
  27. De Leva, P. Adjustments to zatsiorsky-seluyanov’s segment inertia parameters. J. Biomech. 1996, 29, 1223–1230. [Google Scholar] [CrossRef]
  28. Barkan, U.; Masayoshi, N.; Kazuyuki, H.; Michihiro, K.; Tatsuo, N. Proof of Concept for Robot-Aided Upper Limb Rehabilitation Using Disturbance Observers. IEEE Trans. Hum.-Mach. Syst. 2015, 45, 110–118. [Google Scholar]
  29. Mohammed, S.; Huo, W.; Huang, J.; Rifaï, H.; Amirat, Y. Nonlinear disturbance observer based sliding mode control of a human-driven knee joint orthosis. Robot. Auton. Syst. 2014, 75, 41–49. [Google Scholar] [CrossRef]
  30. Bruno, S.; Lorenzo, S.; Luigi, V.; Giuseppe, O. Kinematics BT-Robotics: Modelling, Planning and Control; Springer: London, UK, 2009; pp. 39–103. [Google Scholar]
  31. Fu, K.; Gonzalez, R.; Lee, C. Robotics: Control, Sensing, Vision and Intelligence; McGraw-Hill: New York, NY, USA, 1987; p. 1. [Google Scholar]
  32. Lewis, F.L.; Dawson, D.M.; Abdallah, C.T. Robot Manipulator Control Theory and Practice; CRC Press: Boca Raton, FL, USA, 2004. [Google Scholar]
  33. Denavit, J.; Hartenberg, R.S. A Kinematic Notation for Lower-Pair Mechanisms Based on Matrices. J. Appl. Mech. 1955, 22, 215–221. [Google Scholar] [CrossRef]
  34. Mohammadi, A.; Tavakoli, M.; Marquez, H.J.; Hashemzadeh, F. Nonlinear disturbance observer design for robotic manipulators. Control Eng. Pract. 2013, 21, 253–267. [Google Scholar] [CrossRef] [Green Version]
  35. Chen, W.H.; Ballance, D.J.; Gawthrop, P.J.; O’Reilly, J. A nonlinear disturbance observer for robotic manipulators. IEEE Trans. Ind. Electron. 2000, 47, 932–938. [Google Scholar] [CrossRef] [Green Version]
  36. Riani, A.; Madani, T.; Hadri, A.E.; Benallegue, A. Adaptive control based on an on-line parameter estimation of an upper limb exoskeleton. In Proceedings of the IEEE International Conference on Rehabilitation Robotics, London, UK, 17–20 July 2017; pp. 695–701. [Google Scholar] [CrossRef]
  37. Khalil, W.; Dombre, E. Modeling, Identification and Control of Robots; Butterworth-Heinemann: Oxford, UK, 2004; pp. 1–480. [Google Scholar] [CrossRef]
  38. Gautier, M.; Khalil, W. Direct Calculation of minimum set of inertial Parameters of Serial Robots. IEEE Trans. Robot. Autom. 1990, 6, 368–373. [Google Scholar] [CrossRef] [Green Version]
  39. Shimkin, N.; Feuer, A.F. Persistency of excitation in continuous-time systems. Syst. Control Lett. 1987, 9, 225–233. [Google Scholar] [CrossRef]
  40. Nancy Berryman, W.D.B.P. Joint Range of Motion and Muscle Length Testing, 1st ed.; Saunders: Philadelphia, PA, USA, 2002. [Google Scholar]
  41. Serway, R.A.; Jewett, J.W.; Peroomian, V. Physics for Scientists and Engineers; Saunders: Philadelphia, PA, USA, 2018. [Google Scholar]
  42. Yepes, J.C.; Saldarriaga, A.J.; Vélez, J.M.; Pérez, V.Z.; Betancur, M.J. A Hardware-in-the-loop Simulation Study of a Mechatronic System for Anterior Cruciate Ligament Injuries Rehabilitation. In Proceedings of the BIODEVICES 2017—10th International Conference on Biomedical Electronics and Devices, Porto, Portugal, 21–23 February 2017; Volume 1. [Google Scholar]
  43. Yepes, J.C. Surface Electromyography Signal Processing Algorithm and Movement Control Algorithm for Mechatronic-Assisted Rehabilitation of Anterior Cruciate Ligament Injuries. Ph.D. Thesis, Universidad Pontificia Bolivariana, Medellin, Colombia, 2018. [Google Scholar] [CrossRef]
  44. Landau, I.D.; Lozano, R.; M’Saad, M.; Karimi, A. Discrete-Time System Models for Control; Springer: London, UK, 2011; pp. 35–53. [Google Scholar] [CrossRef]
  45. Vakanski, A.; Jun, H.P.; Paul, D.; Baker, R. A Data Set of Human Body Movements for Physical Rehabilitation Exercises. Data 2018, 3, 2. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  46. Patiño, J.G.; Bravo, E.E.; Perez, J.J.; Perez, V. Lower Limb Rehabilitation System Controlled by Robotics, Electromyography Surface and Functional Electrical Stimulation; Pan American Health Care Exchanges, PAHCE: Medellin, Colombia, 2013; p. 6257. [Google Scholar] [CrossRef]
  47. Kirby, F.A. Simulación de los Algoritmos de Control de un Sistema de Rehabilitación de Miembro Inferiro (LegSys). Master’s Thesis, Universidad Pontificia Bolivariana, Medellin, Colombia, 2016. [Google Scholar]
  48. Bustamante, S.; Yepes, J.C.; Pérez, V.Z.; Correa, J.C.; Betancur, M.J. Online Simulation of Mechatronic Neural Interface Systems: Two Case-Studies. In Proceedings of the Biomedical Engineering Systems and Technologies: 9th International Joint Conference, BIOSTEC 2016, Rome, Italy, 21–23 February 2016; Revised Selected Papers. Springer International Publishing: Cham, Switzerland, 2017; pp. 255–275. [Google Scholar] [CrossRef]
  49. Zoss, A.B.; Kazerooni, H.; Chu, A. Biomechanical Design of the Berkeley Lower Extremity Exoskeletong (BLEEX). IEEE/ASME Trans. Mechatron. 2006, 11, 128–138. [Google Scholar] [CrossRef]
  50. Vaquero, A.F.; Chicharro, J.L. Fisiología del Ejercicio. In Fisiología del Ejercicio; Editorial Médica Panamericana: Bogotá, Colombia, 2006; p. 987. [Google Scholar]
Figure 1. LLRR configuration.
Figure 1. LLRR configuration.
Applsci 12 05529 g001
Figure 2. Block diagrams for the IID. (a) Calibration phase and (b) torque estimation phase.
Figure 2. Block diagrams for the IID. (a) Calibration phase and (b) torque estimation phase.
Applsci 12 05529 g002
Figure 3. Block diagrams for the INDO. (a) Calibration phase and (b) torque estimation phase.
Figure 3. Block diagrams for the INDO. (a) Calibration phase and (b) torque estimation phase.
Applsci 12 05529 g003
Figure 4. LLRR vector parameters error χ ˜ during the estimation phase.
Figure 4. LLRR vector parameters error χ ˜ during the estimation phase.
Applsci 12 05529 g004
Figure 5. Simulated squat sequence.
Figure 5. Simulated squat sequence.
Applsci 12 05529 g005
Figure 6. LLRR dynamics during a squat exercise.
Figure 6. LLRR dynamics during a squat exercise.
Applsci 12 05529 g006
Figure 7. Comparison of (a) IID and Saadatzi et al. [12] ID method, and (b) INDO and Saadatzi et al. [12] NDO method, during a squat exercise with SNR = 40 dB and 20% BSIPs uncertainties.
Figure 7. Comparison of (a) IID and Saadatzi et al. [12] ID method, and (b) INDO and Saadatzi et al. [12] NDO method, during a squat exercise with SNR = 40 dB and 20% BSIPs uncertainties.
Applsci 12 05529 g007
Figure 8. Simulated leg press sequence.
Figure 8. Simulated leg press sequence.
Applsci 12 05529 g008
Figure 9. LLRR dynamics during a leg press exercise.
Figure 9. LLRR dynamics during a leg press exercise.
Applsci 12 05529 g009
Figure 10. Comparison of (a) IID and Saadatzi et al. [12] ID method, and (b) INDO and Saadatzi et al. [12] NDO method, during a leg press exercise with SNR = 40 dB and 20% BSIPs uncertainties.
Figure 10. Comparison of (a) IID and Saadatzi et al. [12] ID method, and (b) INDO and Saadatzi et al. [12] NDO method, during a leg press exercise with SNR = 40 dB and 20% BSIPs uncertainties.
Applsci 12 05529 g010
Table 1. D-H Parameters for our 3-DOF LLRR.
Table 1. D-H Parameters for our 3-DOF LLRR.
i α i a i d i q i
10 1 0 q 1
20 2 0 q 2
30 3 0 q 3
Table 2. LLRR parameters for simulation including a subject of h = 1.75 m height and m = 75 kg weight.
Table 2. LLRR parameters for simulation including a subject of h = 1.75 m height and m = 75 kg weight.
ParameterDescriptionValueUnits
1 s Thigh length 24.05 % × h = 0.4209 m
2 s Shank length 24.85 % × h = 0.4349 m
3 s Foot length 13.15 % × h = 0.2301 m
1 r Robot’s link 1 length0.4209 m
2 r Robot’s link 2 length0.4349 m
3 r Robot’s link 3 length0.2301 m
1 Robot + Subject’s link 1 length0.4209 m
2 Robot + Subject’s link 2 length0.4349 m
3 Robot + Subject’s link 3 length0.2301 m
m 1 s Thigh weight 14.47 % × m = 10.8525 k g
m 2 s Shank weight 4.57 % × m = 3.4275 k g
m 3 s Foot weight 1.33 % × m = 0.9975 k g
m 1 r Robot’s link 1 weight19 k g
m 2 r Robot’s link 2 weight9 k g
m 3 r Robot’s link 3 weight11 k g
m 1 Robot + Subject’s weight of link 129.8525 k g
m 2 Robot + Subject’s weight of link 212.4275 k g
m 3 Robot + Subject’s weight of link 311.9975 k g
b 1 Robot + Subject’s CoM of link 1 38.53 % × 1 = 0.1622 m
b 2 Robot + Subject’s CoM of link 2 44.37 % × 2 = 0.193 m
b 3 Robot + Subject’s CoM of link 3 42.14 % × 3 = 0.1382 m
I c 1 , z z Robot + Subject’s Inertia of link 1 1 3 m 1 ( 2 b 1 ) 2 = 1.047 k g   m 2
I c 2 , z z Robot + Subject’s Inertia of link 2 1 3 m 2 ( 2 b 2 ) 2 = 0.6171 k g   m 2
I c 3 , z z Robot + Subject’s Inertia of link 3 1 3 m 3 ( 2 b 3 ) 2 = 0.3053 k g   m 2
v 1 Viscous Friction in hip joint100 N   m   s rad
v 2 Viscous Friction in knee joint100 N   m   s rad
v 3 Viscous Friction in ankle joint60 N   m   s rad
g 0 Gravity9.8 m s 2
K p Proportional gain150 % p u
T i Integral time0 s
T d Derivative time0.1 s
s a t 1 Motor 1 saturation768.458 N   m
s a t 2 Motor 2 saturation371.377 N   m
s a t 3 Motor 3 saturation102.689 N   m
Table 3. Human–robot interaction torque estimation methods performance indices during a squat exercise with SNR = 40   d B and 20% BSIPs uncertainties. Bold numbers highlight the best metric value.
Table 3. Human–robot interaction torque estimation methods performance indices during a squat exercise with SNR = 40   d B and 20% BSIPs uncertainties. Bold numbers highlight the best metric value.
MethodError (Hip, Knee, Ankle)
MAEMAPE (%)RMSERMSPE (%) R 2
Saadatzi et al. [12] ID(73.0, 22.2, 6.67)(1192, 362, ∼)(312, 128, 6.75)(5105, 2095, ∼)(6.64 × 10−5, 0.002, ∼)
IID(37.9, 15.9, 0.468)(618, 260, ∼)(243, 100, 0.691)(3973, 1634, ∼)(3.74 × 10−4, 0.002, ∼)
Saadatzi et al. [12] NDO(31.3, 4.99, 6.65)(512, 81.4, ∼)(33.7, 5.47, 6.77)(550, 89.3, ∼)(0.047, 0.855, ∼)
INDO(1.04, 0.953, 0.814) (16.9, 15.5, ∼) (1.34, 1.22,1.02) (21.8, 19.9, ∼) (0.931, 0.939, ∼)
Table 4. Human–robot interaction torque estimation methods performance indices during a leg press exercise with SNR = 40 dB and 20% BSIPs uncertainties. Bold numbers highlight the best metric value.
Table 4. Human–robot interaction torque estimation methods performance indices during a leg press exercise with SNR = 40 dB and 20% BSIPs uncertainties. Bold numbers highlight the best metric value.
MethodError (Hip, Knee, Ankle)
MAEMAPE (%)RMSERMSPE (%) R 2
Saadatzi et al. [12] ID(40.9, 21.3, 6.59)(696, 361, ∼)(153, 73.2, 17.2)(616, 1246, ∼)(0.004, 0.004, ∼)
IID(21.3, 9.90, 2.03)(362, 168, ∼)(117, 55.2, 11.8)(1998, 939, ∼)(0.002, 0.007, ∼)
Saadatzi et al. [12] NDO(21.2, 12.2, 4.60)(360, 207, ∼)(26.1, 12.2, 4.65)(444, 208, ∼)(0.226, 0.931, ∼)
INDO (0.718, 0.609, 0.521) (12.2, 10.4, ∼) (0.993, 0.803, 0.650) (16.9, 13.6, ∼) (0.959, 0.973, ∼)
Table 5. Evaluation of the performance of the human–robot interaction torque estimation methods, during simulations with SNR = 40 dB and 20% BSIPs uncertainties, according to the requirements defined in Section 2. The checkmark and cross-mark symbols represent whether the requirement is fulfilled or not, respectively.
Table 5. Evaluation of the performance of the human–robot interaction torque estimation methods, during simulations with SNR = 40 dB and 20% BSIPs uncertainties, according to the requirements defined in Section 2. The checkmark and cross-mark symbols represent whether the requirement is fulfilled or not, respectively.
ItemRequirementInspired by(Saadatzi et al., 2018) ID MethodIID(Saadatzi et al., 2018) NDO MethodINDO
1Non-dependence of additional sensors[11,13]
2Low phase lag in the estimation[15]
3Low sensor noise sensitivity[15]
4Low sensitivity to BSIPs uncertainties[12]
5The average percentage error must be lower
than 20 to 22%
when using accurate model parameters
[12]
6Small error band[10]
7Approximately 0.5 s of settling time or lower[10]
8Overshoot of approximately 25% or lower within the estimation[10]
9 R 2 greater or approximately 0.935 for the hip joint and 0.924 for the knee joint[11]
10%RMSE lower than 8.74% for
the hip joint and 10.26% for the knee joint
[11]
11A maximum of 5% error when moving
just one joint, i.e., the distal one
[28]
12Finite-time convergence[29]
13Should not require calibration each time
that the user wears the LLRR, i.e., it requires
a maximum of one calibration per user.
[13]
14It works in all ROM, and the limiting
angles of the joints must be configurable
Figure A2a
15It works with slow exercises, executing between 1–25 repetitions per setFigure A2b
16It works with the following ranges of forces: 0   k g
to 15   k g for the hip, 0   k g to 15   k g
for the knee, and  0   k g to 10   k g for the ankle.
Figure A2c
17It works having a maximum percentage of error in the range from
1% to 3%
Figure A2d
18It works within a squat exerciseCase Study
Total passed  5 (27.8%)13 (72.2%)6 (33.3%)16 (88.9%)
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Yepes, J.C.; Rúa, S.; Osorio, M.; Pérez, V.Z.; Moreno, J.A.; Al-Jumaily, A.; Betancur, M.J. Human-Robot Interaction Torque Estimation Methods for a Lower Limb Rehabilitation Robotic System with Uncertainties. Appl. Sci. 2022, 12, 5529. https://doi.org/10.3390/app12115529

AMA Style

Yepes JC, Rúa S, Osorio M, Pérez VZ, Moreno JA, Al-Jumaily A, Betancur MJ. Human-Robot Interaction Torque Estimation Methods for a Lower Limb Rehabilitation Robotic System with Uncertainties. Applied Sciences. 2022; 12(11):5529. https://doi.org/10.3390/app12115529

Chicago/Turabian Style

Yepes, Juan C., Santiago Rúa, Marisol Osorio, Vera Z. Pérez, Jaime A. Moreno, Adel Al-Jumaily, and Manuel J. Betancur. 2022. "Human-Robot Interaction Torque Estimation Methods for a Lower Limb Rehabilitation Robotic System with Uncertainties" Applied Sciences 12, no. 11: 5529. https://doi.org/10.3390/app12115529

APA Style

Yepes, J. C., Rúa, S., Osorio, M., Pérez, V. Z., Moreno, J. A., Al-Jumaily, A., & Betancur, M. J. (2022). Human-Robot Interaction Torque Estimation Methods for a Lower Limb Rehabilitation Robotic System with Uncertainties. Applied Sciences, 12(11), 5529. https://doi.org/10.3390/app12115529

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop