Impedance Control for Robotic Rehabilitation: A Robust Markovian Approach

The human-robot interaction has played an important role in rehabilitation robotics and impedance control has been used in the regulation of interaction forces between the robot actuator and human limbs. Series elastic actuators (SEAs) have been an efficient solution in the design of this kind of robotic application. Standard implementations of impedance control with SEAs require an internal force control loop for guaranteeing the desired impedance output. However, nonlinearities and uncertainties hamper such a guarantee of an accurate force level in this human-robot interaction. This paper addresses the dependence of the impedance control performance on the force control and proposes a control approach that improves the force control robustness. A unified model of the human-robot system that considers the ankle impedance by a second-order dynamics subject to uncertainties in the stiffness, damping, and inertia parameters has been developed. Fixed, resistive, and passive operation modes of the robotics system were defined, where transition probabilities among the modes were modeled through a Markov chain. A robust regulator for Markovian jump linear systems was used in the design of the force control. Experimental results show the approach improves the impedance control performance. For comparison purposes, a standard H∞ force controller based on the fixed operation mode has also been designed. The Markovian control approach outperformed the H∞ control when all operation modes were taken into account.


INTRODUCTION
Physical therapy represents a well-accepted procedure for improvements in the recovery of human motor function and promotion of higher performance in Activities of Daily Life (ADLs) (Krebs et al., 2008). Mainly when people have been affected by injuries such as stroke (Hatano, 1976) and Multiple Sclerosis (MS) (Cattaneo et al., 2002). Robotic-assisted therapy is a promising field for the development of rehabilitation tasks. Among the advantages offered by robotic devices are uniformity in the repetition of long-time routines, reliable records of measured variables, and motivation for the patient's participation using interactive environments like serious games (Lum et al., 2002;Chang and Kim, 2013;Gonçalves et al., 2014). However, the fact that these robots interact with humans during therapeutic movement, they require a high degree of security and reliability. Therefore, rehabilitation robots should identify activities performed by the patient to reach only pre-defined training objectives, whose principle is known as Assist-as-Needed Paradigm (Radomski and Trombly, 2013). Impedance control has been used in the implementation of this kind of paradigm in rehabilitation robotic systems. It was initially proposed for manipulator robots to obtain a safe physical interaction with the environment (Hogan, 1985). Such a controller aims at establishing a dynamic relationship between the force and velocity of an actuator. Series elastic actuators (SEA) provide a simple and efficient solution for the implementation of impedance controllers (Calanca et al., 2016). However, the impedance control of SEAs requires an explicit force control loop whose performance is sensitive to uncertainties and time-varying human dynamics. Hogan (1988, 1989) addressed the problem of interaction control and analyzed the energy exchange between a robotic system and its environment, defined as passive. They also determined stability criteria for a coupled system. Some years later this analysis was brought into the context of humanrobot interaction considering three new issues: (1) Human dynamics is now the environment; (2) This dynamics can exhibit passive and active behaviors; and, (3) The stability criterion has been reformulated as complementary stability (Buerger and Hogan, 2007). Such new concepts emphasize the way the human dynamics is taken into account. For example, Vallery et al. (2008) analyzed limits of coupled stability and performance of an SEA actuator for rehabilitation applications considering the human dynamics a spring and damping system. Kong et al. (2009) analyzed those limits considering the human dynamics only as a mass (inertia). In Oh and Kong (2017), the human dynamics was considered as stiffness, damping, and inertia. Li et al. (2017) proposed an adaptive control scheme for SEA-driven robots which consider two operation modes in the adaptation process, namely robot-in-charge and human-in-charge. Similar approaches were proposed in Yu et al. (2015) and Pan et al. (2017), however, the authors did not take into account the fact the human-robot system can be modeled by different operation modes related to abrupt changes in the dynamic behavior.
This paper reports on the implementation of an impedance controller in a robotic platform for ankle rehabilitation based on a Markovian approach. The platform uses an SEA and enables plantarflexion and dorsiflexion movements. The following three operation modes that may occur in the robot-human interaction were defined: (1) fixed mode, in which the platform is mechanically fixed; (2) resistive mode, in which the user makes efforts against the platform movement; and (3) passive mode, in which the user does not make any effort against the platform movement. Such operation modes are modeled as states of a Markov chain. Based on this modeling, a recursive robust regulator for discrete Markov jump linear systems (RR-DMJLS) proposed in Cerri and Terra (2017) is designed to regulate the force control. It guarantees mean square stability and optimal performance for this class of stochastic system (Jutinico et al., 2017). In order to check the effectiveness of the approach, we performed a comparative study with a standard H ∞ force controller proposed in Pérez-Ibarra et al. (2017) which is designed based only on the higher-impedance operation mode (fixed). Although this control approach provides robust stability for the whole system, including all operation modes, its performance was outperformed by RR-DMJLS. We present actual results based on force-and impedance-control for both controllers.
The paper is organized as follows: Section 2 introduces the SEA-based robotic platform and its respective dynamic model; Section 3 describes the design of the robust controllers; Section 4 reports experimental results of a comparative study between the controllers; finally, Section 5 provides the conclusions and some final remarks for future work.

SYSTEM DESCRIPTION AND MODELING
The SEA-based Robotic Platform for Ankle Rehabilitation (SRPAR), Figure 1, is a device for robot-assisted training. It works under two conditions: guidance of a physiotherapist during pre-established dorsi/plantarflexion movements of the ankle and active participation of the patient by using serious games. Both conditions benefit individuals who have suffered a stroke. Also, the platform is a tool to evaluate the ankle force and range of movement (Gonçalves et al., 2013).
The platform uses a DC motor (Maxon Motor RE 40, graphite brushes, 150 Watt DC motor) linked to a ballscrew through a belt and pulleys with 2:1 reduction ratio. A recirculating ballscrew nut converts the rotational motion of the screw in a linear motion. A pair of steel springs is attached to the nut and to a movable piece. When the motor is driven, the nut moves forward or backward compressing the springs. The movable part is connected to a kinematic chain which converts linear force into torque which is transferred to the user. We estimate the nut's position by measuring the motor rotation with a magneto-resistant incremental encoder. Finally, we estimate the spring force by measuring the spring deformation. A logarithmic sliding potentiometer is attached to the movable piece located between the springs. The potentiometer's cursor moves along with the piece generating a voltage proportional to the spring deformation.

Transfer Function Model
In order to describe the dynamic behavior of the system, we define three sub-systems in the SRPAR: motor-transmission system, series elastic element, and human-load system ( Figure 1B). A set of assumptions is made to simplify the dynamic modeling, resulting in a linear model for the SRPAR. Although these simplifications can result in a limited model, mainly due to the presence of nonlinearities, the use of robust controllers to deal with uncertain parameters can improve the performance of the system.
Concerning the motor-transmission sub-system, although some studies deal with only the effect of the inertia in the motor model (Vallery et al., 2008;Calanca et al., 2014), in this paper we also consider the influence of the motor damping. The effects of the inertia and damping parameters of the pulley and ballscrew are also considered as in Yu et al. (2013) and dos Santos et al. (2015). In addition, the nonlinear effects of friction, backslash, and efficiency losses in the motor-transmission system are minimized by controlling the motor velocity instead of directly controlling the motor torque or current, as discussed in Wyeth (2006).
Regarding the series elastic configuration, in Petit et al. (2015) is presented a generic model for robotic systems with variable stiffness. They describe the output torque of the actuator by a nonlinear function that depends on the joint deflection and mechanical stiffness variation of the springs. In this paper, the springs are modeled as constant stiffness and operate in a linear region. Since the human limb is attached to the four-bar mechanism of the SRPAR, rotation of the ankle is transformed into a linear movement in the same direction of the spring force by a nonlinear Jacobian. However, due to a small range of allowed movements related to mechanical constraints of the SRPAR, this transformation is approximated by a linear relationship (see Equation 6). The gravity effects on the human foot are also neglected since the distance between the foots center of gravity and the rotation axis of the platform is small.
Different approaches have been proposed to model the dynamics of the human-load sub-system. For example, Tagliamonte and Accoto (2014) used a set of second-order plants to represent human dynamics in order to implement an impedance controller where passivity concepts are explored. In Lee et al. (2016), similar second-order models are used to measure the ankle mechanical impedance in a unified way. In this paper, we model the human dynamics as a set of second-order linear plants subject to parameter uncertainties.
In the following, we present a transfer function between the commanded motor velocity ω d m and the spring force F s .

Dynamics of the Motor-Transmission System
Motor-transmission can be modeled as a second-order control system, where X w denotes the displacement of ball screw nut, M meq and B meq are respectively the equivalent inertia and damping of the system as defined in Equation (2), and F meq is the output force as defined in Equation (3), (2) In Equations (2) and (3), M t and B t are the mass and damping of the ball screw nut, J and C are torsional inertia and damping where subscripts m and p stand for motor and pulley, respectively. N p is the pulley ratio, ρ = 2π l is a rotational-to-linear factor where l is the ball screw lead. K t is the motor constant, i m is the motor current determined by the inner velocity control of the motor, and K p and K i are the proportional and integral gains of the controller; ω m and ω d m are actual and desired motor velocities, respectively.
where M l , B l , and K l are respectively equivalent inertia, damping and stiffness of the human-load system, defined by:

Dynamics of the Series Elastic Actuator
Taking the Laplace transform of Equations (1) and (7), and solving for X w and X l , we obtain: (9) We make use of the Hooke's Law to define output spring force F s : From Equations (9) and (10), we have: where Z meq = M meq s + B meq and Z l = M l s + B l + Kl s are the mechanical impedances of the motor-transmission and human-load system, respectively. Finally, the spring force F s (s) is expressed as a function of the desired motor velocity ω d m and human torque τ h , as: thus, the system dynamics is defined by the transfer functions G n (s) and G h (s), given by: where K PI = K p + K i s . The nominal SEA model is obtained fixing the output load making X l = 0 (Robinson et al., 1999;dos Santos et al., 2015). Thus making Z l → ∞ in Equation (13), we obtain a transfer function G(s) that only considers platform parameters: G(s) = G n (s) Z l →∞ = ρN p K s K t K PI K PI M meq s 2 + K PI B meq s + K s .

State Space Model
Consider again the system shown in Figure 1B. In order to simplify the model, the inner velocity loop control allows us to model the motor as a pure velocity source; therefore the torque of the motor τ m is given by: and, in consequence, From Equations (1) and (17), and taking into account that the angular position of the motor φ m is described in function of the displacement X w by φ m = ρN p X w , we obtain: (18) From Equations (7) and (10), we obtain the following expression: Using Equations (6), (8) where F s is the spring force and, φ m and φ l are angular positions for motor and load, respectively. The system control input is the motor angular velocity w m and τ h is the human torque which is considered as an input disturbance.

Experimental Validation
In order to validate the proposed theoretical model, we identify the frequency response function (FRF) of the force modeled in Equation (13). We consider in the human-load interaction, three specific operation modes: (1) a fixed mode, in which the platform is fixed in a neutral position, i.e., φ l ≈ 0 and Z l → ∞; (2) a resistive mode, where the human being makes effort against the platform in order to hold it in the neutral position; and (3) a passive mode, when the user leaves the platform leads the movement. For all modes, we apply a desired motor velocity given by a chirp signal with an amplitude of 209.4 rad/s (2,000 rpm) and sweeping frequencies between 0 and 20 Hz (Figure 2). Due to variations in the activation level of the muscles acting on the ankle joint, abrupt changes in ankle mechanical impedance are expected. In addition, operation modes defined are properly matched with the phases of human gait pattern. Thus, the fixed mode can be associated with the mid-stance phase; the resistive mode with the initial contact, terminal stance and double support of the stance phase; and the passive mode with the swing phase. Table 1 shows the platform specifications, as well as the human parameters used for each mode and their corresponding lower limits. Parameters for Mode 1 were chosen to satisfy Z l → ∞, in order to obtain similar results to those presented in Robinson et al. (1999). For Modes 2 and 3, they were chosen from Lee et al. (2016) taking into account the actual human impedance during stance and swing phases of walking, respectively.

CONTROL APPROACHES
The block diagram of the control system for the SRPAR is presented in Figure 3. The platform uses an EPOS driver that performs two internal control loops for motor velocity and current. They are based on PI controllers with Kp and Ki gains shown in (see Table 1). Platform sensors provide data to a signal conditioning block where SRPAR-human system variables are computed. The sampling frequency used in this process is 1 kHz.
In order to ensure an appropriate interaction between human and platform, we implement a standard impedance control configuration for SEAs that uses the following internal force control law: where F d is the desired force computed from a sequence of desired trajectories for load angular position φ d l and velocity ω d l . It is determined by the virtual stiffness K v and damping B v .
We aim to improve the performance and to guarantee the stability of the system for different modes of human activities, where changes among them are modeled as random jumps. In this sense, we design a recursive robust regulator for discrete-time Markov jump linear systems (RR-DMJLS). Disturbances and uncertainties due to human-robot interactions, mainly because FIGURE 2 | Experimental identification. Frequency response measurements of theoretical (dashed) and experimental (solid) transfer functions between input ω d m and output F s . Graphs show responses for passive (blue), resistive (red) and fixed (black) modes.
of human parameters (inertia, stiffness, and damping), are considered in this control approach. For comparison purposes, we also design an H ∞ force controller that does not take into account these different modes of human activities. We synthesize a fixed-gain controller using a similar approach presented in Mehling and O'Malley (2014) and dos Santos et al. (2015). Both controllers are presented in the following.

Recursive RR-DMJLS Force Control Design
In this section, we design a robust force controller for the SRPAR ( Figure 4A). We present first nominal and uncertain representations for different operation modes of the system. Then, we model the system as a discrete-time Markovian jump linear system (DMJLS). In order to guarantee stability and robust performance, we use the recursive RR-DMJLS algorithm developed in Cerri and Terra (2017).

Nominal Model
Consider the model presented in Equation (20), for each operation mode θ ∈ : = {1, ..., s}, where s is the number of nominal models, F a θ ∈ R nxn , B a θ ∈ R nxm , and G a θ ∈ R nxm are nominal parameter matrices, x a ∈ R n is the state vector, w m ∈ R m is the input control and τ h ∈ R m is the input disturbance. We discretized this model by using F a θ ,k = I + F a θ T s , To eliminate the steady state error, we augment the model by including an integral action: where r k is a force reference signal and C a = [0 − 1 0 0]. F θ ,k and B θ ,k are nominal matrices for three nominal models according to Equation (20), which are based on parameters presented in Table 1: Frontiers in Neurorobotics | www.frontiersin.org

Uncertain Model
Consider the system presented in Section 3.1.1 subject to parametric uncertainties given by: Uncertain matrices δF θ ,k and δB θ ,k are modeled by: where H θ ,k ∈ R nxk is a nonzero matrix, E F θ ,k ∈ R lxn and E B θ ,k ∈ R lxm are known matrices, θ ,k is an arbitrary matrix that satisfies || θ ,k || ≤ 1. In order to identify matrices H θ ,k , E F θ ,k and E B θ ,k , we analyze frequency responses of the nominal and uncertain models described in Equations (24) and (28), where τ h k = 0, by: x n (z) = (zI − F θ ,k ) −1 (B θ ,k u(z) + Br θ ,k r(z)), with z = e jωT s ≈ 1+jωT s /2 1−jωT s /2 . For each operation mode, we compute a transfer function W a,θ in order to satisfy: where σ G n θ and σ G un θ are singular values of the nominal and uncertain models, respectively, see Figure 5. Notice that upper bounds defined by singular values of W a,θ are effective for frequencies until 1.7 Hz for resistive mode and 3.8 Hz for passive mode.
Finally, matrices H θ ,k , E F θ ,k , and E B θ ,k are found using a genetic algorithm considering Equation (32)

Probability Matrix
Time transitions among modes defined in Section 2.3 for the system presented in Equation (28) can be modeled by a Markov chain {θ k } N −1 k=0 , where θ k is called the jump parameter and belongs to a finite set : = {1, ..., s}. Transitions among these different modes are determined by a probability matrix for state transitions P = [p ij,k ] ∈ R s×s where each input must satisfies the following constraints: Different intervals for load velocity and position can describe different behaviors of the human-robot system. Namely, if the user is passive (1), resistive (2) or if the platform is fixed (3), ω l and φ l provide the jump parameter where, for each k = 0, ..., N − 1, where w f l k is the low-pass filtered signal from load velocity absolute value, with f c = 0.1 Hz and T s = 1 ms, given by: Since w l k is related to the frequency response of the system, the jump parameter θ k is also related to the robustness of the regulator. In fact, robust stability and optimal performance is only guaranteed in the interval of frequencies in which relative errors are represented by upper bounds ||σ W a,i (jω) || (Figure 5).
FIGURE 5 | Relative errors and upper bounds. Relative errors between nominal and uncertain models with respect to inputs u k (green) and r k (blue). Graphs also show upper bounds for inputs u k (red) and r k (black).

RR-DMJLS Algorithm
Equation (28) is rewritten as a discrete-time Markov jump linear system subject to parametric uncertainties, by: for all k = 0, ..., N − 1, where F θ k ,k and B θ k ,k are nominal parameter matrices for each mode given by Equation (25), x k is the state vector, u k is an input control, and δF θ k ,k and δB θ k ,k are uncertain matrices modeled as: according to parameters presented in Equations (33) and (34). Consider the robust control problem to regulate the DMJLS subject to parametric uncertainties defined in Equation (39). The solution for this problem is achieved through the min-max optimization problem: for each k = N − 1, ..., 0 and θ k ∈ , where J µ K is the uncertain penalized-regularized quadratic functional: where F δ θ k ,k and B δ θ k ,k are defined in Equation (28); θ k ,k+1 = s j=1 P j,k+1 p ij,k ; P θ k ,k is a positive definite matrix; Q θ k ,k and R θ k ,k are semi-definite weighting matrices. The solution to the optimization problem expressed in Equations (41) and (42) that guarantees the optimal state-control sequence {(x * µ,k+1 , u * µ,k )} N−1 k=0 for a fixed instant k and state θ k , is given by the following Robust Regulator for Discrete Markov Jump Linear Systems (RR-DMJLS): Robust Regulator for DMJLS (Cerri and Terra, 2017). Initial Conditions: Set x 0 , θ 0 , P, P i (N ) ≻ 0, ∀i ∈ : = {1, ..., s}.
Step 1: (Backward). Calculate for all k = N − 1, . . . , 0: Step 2: (Forward). Obtain for each k = 0, . . . N − 1: Equation (43) uses the following auxiliary matrices: In this formulation, µ > 0 is a penalty parameter responsible to guarantee the robustness of the RR-DMJLS. In fact, when µ → +∞ then W i,k → 0. In consequence, the DMJLS closed-loop response is given by: which provides the robust optimal response (x * k+1 , u * k ). Details of the necessary and sufficient conditions for existence of the mean square stabilizing solution and robustness of this regulator can be found in Cerri and Terra (2017).
Let µ = 9.998 × 10 6 and λ i,k = 1 × 10 17 in order to satisfy Equations (44) and (45); weighting matrices R 1,k = R 2,k = R 3,k = 1, Q 1,k = Q 2,k = Q 3,k = I 5 and P 1 (N ) = P 2 (N ) = P 3 (N ) = 1 × 10 10 × I 5 ; and the probability matrix P defined in Equation (36). By using the robust regulator presented in Equation (43), we obtain the following control law: where K a,θ k is the gain to the states x a k and K int,θ k is the gain to the state x int k . Table 2 shows the control gains obtained for three Markovian modes. We do not consider uncertainties in the third term of E F i,k , as a consequence, the algorithm decouples the state variable φ m guaranteeing the controllability of the system (see K 3 ).

H ∞ Force Control Design
Consider the nominal model: We use a mixed-sensitivity shaping approach S = (1 + G n K c ) −1 to ensure tracking performance and disturbance rejection, and K c S to limit the control signal (Skogestad and Postlethwaite, 2007). The H ∞ problem is defined as: where w e (s) and w u (s) are respectively performance and control weights, K c is a stabilizing controller that bounds N by an attenuation level γ , andσ (N) is given by the usual Euclidean vector norm:σ (N) = |w e S| 2 + |w u K c S| 2 .

EXPERIMENTAL RESULTS
In order to evaluate the stability and performance of the proposed control approaches, time and frequency response tests are performed with a healthy user for three different operation modes presented in Section 2.3. We show a comparative study of force controllers and impedance controllers, according to Figure 4. This study was approved and carried out in accordance with the recommendations of the Ethics Committee of the Federal University of São Carlos (Number 26054813.1.0000.5504).

Force Controllers
In this section, robust force controllers presented in Section 3, RR-DMJLS and H ∞ , are compared. Two performance criteria are used for this purpose. We consider the rise time t r and a normalized mean error between spring and desired forces defined as where N = 8/T s is the number of samples. Figure 6 shows time responses of the system using the RR-DMJLS force controller. In these experiments, three different levels (100, 0, and −100 N) of desired forces, F d , are set during a total time of T = 8 s. In tests performed for fixed and resistive modes, shown in Figures 6A,B, the Markovian modes remained constant. Notice in the test shown in Figure 6C, the Markov chain changes between Modes 2 and 3. The RR-DMJLS responses are similar for three modes available, maintaining an appropriate tracking despite natural differences between human and robot dynamics. Notice that t r are similar to the respective desired force levels and the mean errors e qc are 14.24, 11.97, and 11.9%, for tests shown in Figures 6A-C, respectively. Figure 7 shows the force control response of the RR-DMJLS while tracking a sinusoidal force reference. In this experiment, the user is asked to be resistive to the movement of the platform during the first fifteen seconds, and then to remain passive during the next fifteen seconds. We observe that force control maintains a similar response despite the abrupt change in the human dynamics (resistive → passive). Notice that jumps between Markovian states reflect time transitions between different operation modes of the system. However, for a condition where the operation modes change more frequently, these transitions would be detected with a certain delay. We hypothesize that this delay is related to the jump parameter identification method considered (Equation 37). It could be improved by estimating alternative variables. For example, the stiffness and damping parameters of the human being. Figure 8 presents time responses of the system using the H ∞ force control approach. In this case, also are applied three different levels of desired forces F d in T = 8 s, for each operation mode: fixed, resistive, and passive. The rise time t r is lower for fixed and resistive modes in comparison with the passive mode. The mean errors obtained for fixed, resistive and passive modes were 5.89, 8.23, and 23.29%, respectively. Comparing mean  errors between both controllers, we obtain better performance for the H ∞ force controller for the fixed mode. For the passive mode, the RR-DMJLS provides better performance. However, when we design the H ∞ control based on the passive mode, it does not stabilize the system when it is operating in the fixed mode. An advantage of the RR-DMJLS is the uniformity of performance obtained for the whole system. Figure 9 shows frequency responses of the closed-loop system using the H ∞ and recursive RR-DMJLS force controllers. For passive and resistive modes, we apply a desired force signal given by a chirp signal with amplitude of 100 N, sweeping frequencies between 0 and 8 Hz. Some indexes for both controllers are compared. For the H ∞ force controller, shown in Figure 9A, we have: maximum magnitude during passive mode, |G fp | max = −1.6 dB, and resistive mode, |G fr | max = 0.86 dB; bandwidth for passive mode, Bw fp −3dB = 0.2 Hz, and bandwidth for resistive mode, Bw fr −3dB = 1.39 Hz; phase at cut-off frequency for passive mode, α fp = −43 • , and for resistive mode, α fr = −79 • . For the frequency response of the RR-DMJLS shown in Figure 9B, we obtain: |G fp | max = 0.9 dB; |G fr | max = 0.32 dB; Bw fp −3dB = 1.2 Hz; Bw fr −3dB = 1.5 Hz; α fp = −110 • ; and α fr = −75 • . Notice that the bandwidth of the closed-loop system is greater for the RR-DMJLS controller.

Impedance Control
To obtain the RR-DMJLS force control response during impedance controlled movements, we perform two experiments in which the system must track a kinematic reference. In both cases, the user was asked to be resistive during first ten seconds and then to be passive during the next ten seconds. Figure 10A shows the force and impedance control for a pure stiffness control configuration and Figure 10B a stiffnessdamping control configuration. Notice that the torque tracking performance is similar for two modes in spite of the transition between them. Regarding the passive case, the position tracking is better for the pure-stiffness configuration. However, since this configuration has no damping parameter, the velocity tracking is worse. We can include a damping coefficient. However, it can decrease the performance of the position tracking. Thus, there exists a compromise between these control objectives that must be considered by the assistance strategy. Regarding the Markovian states, notice that increasing B v it reduces the velocity of the system and enforces the system to remain in the same Markovian state, Figure 10B. Figure 11 shows time responses using the impedance control with the inner H ∞ force control loop, for K = 15 N·m/rad, following a sinusoidal trajectory for ankle angular position. In  this case, the user is asked to remain passive during the first ten seconds, and then to be resistive during the next ten seconds. Results show that the torque tracking performance is worse when the human being is in passive mode.
In order to quantify the performance achieved by the proposed controllers, we compute the real stiffness and damping parameters of the system. For this purpose, we calculate the torque τ K v generated by the virtual stiffness K v and the torque τ B v generated by the virtual damping B v with where φ e = φ d l − φ l is the angular position error and ω e = ω d l − ω l is the angular velocity error. We compute the root mean square (RMS) errors of the measured variables in the experiments shown in Figures 10, 11. The RMS of the angular error φ e is given by: where N = 10/T s is the number of samples in each test (passive or resistive). In a similar way, we calculate RMS values for angular velocity error ω e , torque generated by virtual stiffness τ K v , and damping τ B v . Notice that, when B v = 0 then τ B v = 0 and τ K v = τ plat . Actual stiffness K r and damping B r are given by: see Table 3. The error between the virtual stiffness and the actual stiffness are calculated, e K v = |(K v − K r )/K v |100% and the error between the virtual damping and the actual damping, We compare results of the impedance control with K v = 15 Nm/rad and B v = 0 Nms/rad. Impedance control with RR-DMJLS presents higher stiffness accuracy in both passive and resistive modes, 15.28 Nm/rad and 16.13 Nm/rad, respectively. With H ∞ controller, it provides 16.44 Nm/rad in passive mode and 7 Nm/rad in resistive mode. Stiffness errors e K v of the RR-DMJLS force control are smaller than the H ∞ force control. This difference can be seen in the passive mode (θ = 3), where the e K v error of the H ∞ is approximately seven times greater than the RR-DMJLS. In the test with K v = 15N·m/ rad and B v = 5N·m·s/ rad, the performance is preserved. Figure 12 shows the frequency response for the impedance which is measured between output torque and angular position error. For passive and resistive modes, we apply a desired angular position trajectory given by a chirp signal with amplitude 0.2 rad, sweeping frequencies between 0 and 8 Hz. In these experiments, impedance controller is defined as pure stiffness configuration; therefore, magnitude of the impedance should be almost constant. Notice that it is guaranteed only until 1 Hz. Figure 12A shows the behavior of the RR-DMJLS-based impedance control. Figure 12B shows the frequency response for the H ∞ -based impedance control. Notice for this controller that the performance decreases when the system operates in the passive mode. We can see that the RR-DMJLS is a more resilient impedance control if compared with H ∞ -based control for different operating modes and desired impedances.

DISCUSSION AND CONCLUSIONS
The model-based robust force control approach was evaluated in a robotic platform for ankle rehabilitation, improving FIGURE 11 | Impedance control with inner H ∞ force control. Graphs show desired (blue) and measured (red) values for the platform torque, τ plat (top), and angular position of the ankle joint (bottom). its impedance control performance. First, the robot-human dynamics was modeled considering the ankle impedance a second-order system with inertia, stiffness, and damping parameters. Since we do not have exact knowledge of the human parameters, our system is subject to parametric uncertainties and we have defined three operation modes related to the humanrobot activity whose time transitions were modeled via a Markov chain. For comparison purposes, an H ∞ force controller was also evaluated. Although it can guarantee coupled stability, the force control performance decreases when the system is in the passive operation mode. We have also designed and implemented an RR-DMJLS for dealing with abrupt changes and system uncertainties by guaranteeing robust mean square stability of the system. Experimental results show RR-DMJLS outperformed the H ∞ force controller.

Related Work
The impedance control configuration used is based on Hogan (1985) and it is aimed at the regulation of the dynamical behavior in the interaction port by variables that do not depend on the environment. The actuator together with the controller are modeled as an impedance, Z r , with velocity inputs (angular) and force outputs (torque). The environment is considered an admittance, Y e , in the interaction port. Hogan (1988, 1989) presented sufficient conditions for the determination of stability of two coupled systems and explained how two physically coupled systems with Z r and Y e with passive port functions can guarantee stability. These concepts have been useful for the implementation of interaction controls for almost three decades. The stability of two coupled systems is given by Z r and Y e eigenvalues and the performance is evaluated by through the impedance Z r . Buerger and Hogan (2007) described a methodology in which an interaction control is designed for a robot module used for rehabilitation purposes. They considered an environment with restricted uncertain characteristics, therefore the admittance is rewritten as Y e (s) = Y n (s) + W(s) (s). The authors also used a second-order dynamics to model the stiffness, damping, and inertia of the human parameters. Complementary stability for interacting systems was defined, where stability is determined by an environment subject to uncertainties. Therefore, a coupled stability problem is considered a robust stability problem.
Regarding the human modeling, the dynamic properties of the lower limbs and muscular activities vary considerably among subjects. This is relevant since SRPAR has been designed for users that suffer diseases that affect the human motor control system, e.g., stroke and other conditions that cause hemiplegia. Typically, such diseases change stiffness and damping in the ankle and knee joints, hence producing spasticity or hypertonia (Lin et al., 2006;Chow et al., 2012). Therefore, the development of a control strategy that guarantees a safe interaction between patient and platform, mainly in virtue of uncertainties related to the human being, is fundamental. Li et al. (2017) and Pan et al. (2017) proposed adaptive control schemes for SEA-driven robots. They considered two operation modes in the adaptation process, namely robot-in-charge and human-in-charge, which are close related to the passive and resistive operation modes, respectively, proposed in this paper. However, the control adaptation is based on changes in the desired position input of the SEA controller and estimation of coordinate accelerations through nonlinear filtering.
In human-robot interaction control systems, the efficiency of the force actuator operation deserves special attention. Although SEAs are characterized by a low output impedance, an important requirement for improving such efficiency is the achievement of a precise and proportional output torque with respect to the desired input. Pratt (2002), Au et al. (2006), Kong et al. (2009), Mehling andO'Malley (2014), and dos Santos et al. (2015) developed force controllers for ankle actuators using SEA. In this paper, we proposed a force control methodology that can deal with system uncertainties and guarantee robust mean square stability. Similar performance was obtained in different tests performed. Accuracies of 98.14% for resistive mode and 92.47% for passive mode were obtained in the pure stiffness configuration. In the stiffness-damping configuration, with K v = 15 and B v = 5, the accuracy obtained in the resistive case was of 97.47% for stiffness and 97.2% for damping, and in the passive case was of 93.94% for stiffness and 94% for damping.
On the contrary, using a fixed-gain control approach based on H ∞ synthesis, the performance was not similar among operation modes. We showed that this strategy can guarantee coupled stability; nevertheless, force control performance decreases when the system is in the passive operation mode. This is reflected in the impedance control accuracy for the pure stiffness configuration, falling from 90.4% in the resistive mode to 46.67% in the passive mode.

Shortcomings and Possible Improvements
In order to control the SRPAR, we proposed a methodology to force control based on RR-DMJLS. It considers a discretetime Markovian model with three states associated with the operating modes of the system. The model was augmented for eliminating the steady state error through the inclusion of an integral action. We also found appropriate uncertain matrices considering frequency responses for relative errors between perturbed and nominal models.
Regarding the observation of the operation modes, in Figure 10A there exists a delay in the jump identification from resistive to passive mode, and in Figure 10B the proposed jump identification method was not even able to identify the jump between modes. This behavior is directly related to virtual damping B v selected since load angular velocity w f l k is lower than 2 rad/s. A possible solution to the problem is the estimation of the virtual stiffness and damping of the human being for the definition of the bounds of the identification method.
Based on our observations of the system behavior, we have defined a probability matrix P that models those transitions among different modes. We hypothesize the probabilities may vary in function of the user's physiological conditions, therefore, our probability matrix can be considered partially or completely uncertain. In a future study, we aim at using our methodology subject to uncertain transition probabilities with the unknown Markov chain proposed in Bortolin and Terra (2016).
Other approaches may improve our methodology. For example, a disturbance observer, as proposed in Kong et al. (2009), could compensate the effect of human torque τ h in Equations (14) and (20). Optimal robust filter for DMJLS (Ishihara et al., 2015) and extended robust Kalman filter proposed in Inoue et al. (2016) could better estimate the states of the SRPAR.

AUTHOR CONTRIBUTIONS
AJ, JJ, FE, and JP conceived research, performed the experiments and the data analysis, drafted the manuscript, and coded the controllers. MT and AS participated in the design of the controllers and in the draft of the manuscript. All authors read and approved the final manuscript.