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

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 dB 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.


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 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τ int despite BSIPs uncertainties. To the best of the authors' knowledge, none of the previous research has examined how a humanrobot 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. Sections 3 and 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.

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:
The average percentage error must be lower than 20 to 22% when using accurate model parameters [12]. 6.
Coefficient of Correlation (R 2 ) greater or approximately 0.935 for the hip joint and 0.924 for the knee joint [11]. 10. Root Mean Squared Percent Error (RMSPE) lower than 8.74% for the hip joint and 10.26% for the knee joint [11]. 11. A maximum of 5% error when moving just one joint, i.e., the distal one [28]. 12. Finite-time convergence [29]. 13. 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]. 14. It works in all ROM, and the limiting angles of the joints must be configurable. 15. It works with slow exercises, executing between 1-25 repetitions per set. 16. 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. 17. It works having a maximum percentage of error in the range from 1% to 3%. 18. It works within a squat exercise.

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 where q ∈ R 3×1 ,q ∈ R 3×1 , andq ∈ 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 τ int ∈ R 3×1 are the disturbance torques, e.g., the torque exerted by the subject. Appendix C presents detailed information of the dynamics.

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]: It is important to note that this ID-based method uses double derivatives to computë q, which significantly increases the sensitivity of the model to sensor noise. Moreover, this method uses the estimatesM(q),Ĉ(q,q), andĜ(q) of the actual M(q), C(q,q), and G(q) dynamic matrices, respectively. Taking into account that and 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 matricesM(q),Ĉ(q,q), andĜ(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: where L(q,q) is the observer gain matrix,τ d is an estimation of the disturbances in a lumped term,M(q),Ĉ(q,q), andĜ(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 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 ofq, as shown below in the disturbance observer auxiliary vector: 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 NDObased 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τ int when there are uncertainties. Therefore, this paper extends the ID-based and NDO-based methods to cope with BSIPs uncertainties.

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.

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 However, the LLRR's dynamics depend in a nonlinear way on some of these parameters, e.g., through the combination I ci,zz + m i r 2 xi . 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)χ (10) 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 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 and 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 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) =  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.

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 are, respectively, the estimates of the vectors χ and Γ.
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 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: 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.

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 τ int = 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 τ int = 0, i.e., the subject should is exerting torque on the robot and the goal of the phase is to estimateτ int . Here, we can observe that the IID is defined as 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 τ int = 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 τ int = 0, i.e., the subject should is exerting torque on the robot and the goal of the phase is to estimateτ int .
Here, we can observe that the INDO is defined as 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, 2.
There exists a positive definite and symmetric matrix Γ such that 3.
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 ∈ Dq = {q ∈ R n : q ≤ q max }, and M(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: Please refer to Mohammadi et al. [34] for a detailed proof of Theorem 2.

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]. 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 ir and lengths ir . In order to define the lengths is and weights m is 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 is 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].

Simulation Parameters
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 2b 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 ML 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 ci,zz = 1 3 m i (2b i ) 2 . Finally, the position of the CoM b i was configured as reported by De Leva et al. [27].

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].

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
Subsequently, we used the MATLAB ® Robot Control Toolbox ™ to solve this LMIbased analysis and, thus, design the calibration phase estimator. The solution with α = 1 for the LMI (26) was γ = 15.691 and K = 0.0212I 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 illconditioning. 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 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 q 0 = (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 1d = 30 sin(0.0430t) + 30 sin(0.2316t) q 2d = 60 sin(0.0938t) + 30 sin(0.0934t) q 2d = 30 sin(0.0594t) + 30 sin(0.1375t) .
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 χ nom of our LLRR along with the new subject, i.e., with the simulation parameters reported in Table 2 However, the initial conditionsχ 0 , using 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 χ nom 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.

Estimation Phase Results
Based on the parameters from Table 2 for the simulation, we can see that Finally, with the identity matrix I 3 ∈ R 3×3 and (25) would turn into Subsequently, we used the MATLAB ® Robot Control Toolbox ™ to solve this LMIbased analysis, obtaining a design for the NDO. The solution with ζ = 30 for (33) was Y = X = 0.0028I 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 thatq d was not more significant than the one a human joint can execute. In addition, we can observe τ dyn , τ int , τ f ri , and τ a during the estimation phase. It is important to note the τ int 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 τ int = (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 8. Simulated leg press sequence. 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 thatq d was not more significant than the one a human joint can execute. In addition, we can denote τ dyn , τ int , τ f ri , and τ a during the estimation phase. It is essential to note that the τ int 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 τ int = (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. . LLRR dynamics during a leg press exercise. 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 Figures 7a and 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 nondependence 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 Figures 7b and 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]. 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 computeq.

Requirements Evaluation
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. 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. Non-dependence of additional sensors [11,13] 2 Low phase lag in the estimation [15] 3 Low sensor noise sensitivity [15] 4 Low sensitivity to BSIPs uncertainties [12] 5 The average percentage error must be lower than 20 to 22% when using accurate model parameters [12] 6 Small error band [10] 7 Approximately 0.5 s of settling time or lower [10] 8 Overshoot 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] 11 A maximum of 5% error when moving just one joint, i.e., the distal one [28] 12 Finite-time convergence [29] 13 Should not require calibration each time that the user wears the LLRR, i.e., it requires a maximum of one calibration per user. [13] 14 It works in all ROM, and the limiting angles of the joints must be configurable Figure A2a 15 It works with slow exercises, executing between 1-25 repetitions per set Figure A2b 16 It works with 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. Figure A2c 17 It works having a maximum percentage of error in the range from 1% to 3% Figure A2d 18 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.

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 humanrobot 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.  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:

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 threelink 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 kg 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.

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 multiplechoice 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.  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 kg, 175 m tall; with right transfemoral amputation (the distal third of the knee), a stump with 18 cm 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 kg to 15 kg for the hip, 0 kg to 15 kg for the knee, and 0 kg to 10 kg 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: 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. (c) (d) 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.

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 = Ic 1,zz + Ic 2,zz + Ic 3,zz + 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 + Ic 2,zz + Ic 3,zz where T and U denote the total KE and PE of the system, respectively. The Euler-Lagrange equations can be defined as d dt where u k is the generalized force associated with the generalized coordinate q i . Therefore, we rewrite the dynamic equations as where c kij = c kji 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