Modeling and Emulating a Physiotherapist's Role in Robot‐Assisted Rehabilitation

In home‐based rehabilitation, one possible approach is haptic teleoperation in which a hospital‐based therapist is haptically linked and tele‐presented to a home‐based patient to effectively simulate traditional in‐hospital therapies over a distance. In this context, this article proposes a learn‐and‐replay (LAR) paradigm that consists of two phases: a therapist‐in‐loop (interactive) phase where the therapist interacts through the haptic teleoperation loop with the patient to perform the cooperative therapy task, and a therapist‐out‐of‐loop (standalone) phase where the therapist's task is played by the patient‐side robot in future repetitions. During the interactive phase, the therapist demonstrates impedance during cooperating with the patient. During the standalone phase, the patient‐side robot is automatically controlled to mimic the therapist's demonstrated impedance which is learned in the interactive phase. The direct force reflection (DFR) architecture is utilized as the control method for the bilateral telerehabilitation system. Case studies involving 1‐degree‐of‐freedom and 2‐degree‐of‐freedom cooperative manipulation tasks are tested for proof of concept. The results show that the impedance of the therapist's arm can be replicated by the patient‐side robot for both tasks and proposed LAR telerehabilitation paradigm that assists the therapist in the rehabilitation procedure to take care of other tasks or attend to other patients.

DOI: 10.1002/aisy.201900181 In home-based rehabilitation, one possible approach is haptic teleoperation in which a hospital-based therapist is haptically linked and tele-presented to a home-based patient to effectively simulate traditional in-hospital therapies over a distance. In this context, this article proposes a learn-and-replay (LAR) paradigm that consists of two phases: a therapist-in-loop (interactive) phase where the therapist interacts through the haptic teleoperation loop with the patient to perform the cooperative therapy task, and a therapist-out-of-loop (standalone) phase where the therapist's task is played by the patient-side robot in future repetitions. During the interactive phase, the therapist demonstrates impedance during cooperating with the patient. During the standalone phase, the patientside robot is automatically controlled to mimic the therapist's demonstrated impedance which is learned in the interactive phase. The direct force reflection (DFR) architecture is utilized as the control method for the bilateral telerehabilitation system. Case studies involving 1-degree-of-freedom and 2-degree-offreedom cooperative manipulation tasks are tested for proof of concept. The results show that the impedance of the therapist's arm can be replicated by the patient-side robot for both tasks and proposed LAR telerehabilitation paradigm that assists the therapist in the rehabilitation procedure to take care of other tasks or attend to other patients.
impedance displayed by the therapist using therapist-side robot (master robot) and then emulate the therapist's behaviour for the task.
In this article, we propose a novel paradigm called learn and replay (LAR) to realize direct bilateral telerehabilitation that encompasses two distinct phases to achieve the time-sharing of a therapist. During the first phase, the therapist interacts directly with a patient through bilateral teleoperation (Figure 1a) to complete a cooperative task. This phase is called the therapistin-loop (interactive) phase. Throughout the interaction, the therapist's task-specific impedance is measured through the therapist-side robot. During the next stage, the therapist is no longer in the rehabilitation loop, thus giving rise to the therapistout-of-loop (standalone) phase ( Figure 1b). In the standalone phase, the therapist's measured arm impedance is displayed to the patient by the patient-side robot via an impedance control loop so that the cooperative rehabilitation task can be carried out in the absence of the therapist.
The mechanical impedance of the human arm (which we aim to teach to our robotic system) is often measured via system identification methods. Either position or force perturbations are applied by a robot to the hand, and the resulting force and motion responses are analyzed to fit typically to a second-order impedance model. For example, a planar robot can be used to impose step position disturbances to the human hand, enabling the calculation of the endpoint stiffness of the arm in the Cartesian plane. [11] Later, dynamic components, i.e., damping and inertia, can be added to the impedance model. [12][13][14] In rehabilitation, the mechanical arm impedance can potentially serve as a quantifiable index to measure patient recovery. Both position and force perturbations can be applied to the arm to identify a second-order impedance model. [15] When force perturbations are applied, the subject is required to maintain a certain posture. While this approach maximizes the model precision and consistency, it is geared toward assessing patient recovery and not therapists. It will be impractical for the interactive phase of our proposed LAR telerehabilitation paradigm, as to restrain the therapist to a certain posture during the entire session of therapy does not facilitate therapy. The endpoint impedance of the human arm can be measured in a real-world task by applying short, impulsive force perturbations to the arm during the execution of the task. [16][17][18] Force perturbations with 3N amplitude and 100 ms duration can be exerted on the hand holding the weld gun. [18] Segments of data immediately following the onset of perturbations with 200 ms of duration were used to identify a second-order impedance model. A modified version of this approach is implemented in this article.
We believe that the proposed paradigm is worth researching because of the following three reasons: 1) It provides the rehabilitation robot with a desired therapy-oriented behaviors by demonstrating the task rather than explicitly programming it through machine commands. This feature is useful in clinical settings where the therapist normally does not have the knowledge to reprogram/reconfigure the robot. 2) It also facilitates time-sharing the same therapist across multiple home-based patients engaged in task-oriented therapy. The therapist can engage the next patient in the interactive phase, whereas the previous patient starts standalone phase exercises. The different therapist arm impedances related to different patients and different therapy tasks will be identified and emulated by different patient-side robots. 3) Finally, the LAR  Figure 1. Illustration of the proposed LAR paradigm with a) interactive phase, where the hospital-based therapist interacts with the home-based patient via a teleoperation system, and b) standalone phase where the therapist's behaviour is emulated by the patient-side robot.
www.advancedsciencenews.com www.advintellsyst.com paradigm is also particularly useful for sophisticated impedance-based rehabilitation tasks that are highly relevant to activities of daily living (ADL); this paradigm can also accomplish simpler trajectory-following tasks. This article is organized as follows: in Section 2, we present a general control architecture involving up to 6 degrees-of-freedom (DOF) for LAR therapy. In Section 3, we introduce the method used for human arm impedance measurement in the interactive phase without interrupting the normal flow of the therapy task. In Section 4, we present the impedance control implementation for the patient-side robot during the standalone phase. In Section 5, we present two case studies involving a 1-DOF task and a 2-DOF task. Finally, Section 6 presents the concluding remarks.

Therapist-Side and patient-Side Robot Kinematics
Consider a bilateral teleoperation system operating in n DOF (n up to 6), comprising a therapist-side robot (placed on the therapist side), a patient-side robot (placed on the patient side), a communication system, and a control architecture. We can denote the joint angles of the therapist-side robot with θ m ¼ ½ θ 1m θ 2m : : : θ nm T and the joint angles of the patientside robot with θ s ¼ ½θ 1s θ 2s : : : θ ns T . The position of the endeffector of the therapist-side and patient-side robots in Cartesian space can be denoted by p m and p s respectively, each being a n Â 1 vector. The velocity in joint space and Cartesian space are linked by the n Â n Jacobian matrices J m and J s for the therapist-side and patient-side robotṡ Joint torques applied by the motors are denoted by τ m ¼ ½τ 1m τ 2m : : : τ nm T for the therapist-side robot and τ s ¼ ½τ 1s τ 2s : : : τ ns T for the patient-side robot. We can therefore relate static joint torques to the end-effector Cartesian forces, f m for the therapist-side robot and f s for the patient-side robot, by Similarly, Cartesian forces applied by human operators onto the robots, f th for the forces applied to the therapist-side robot by the therapist and f p for the forces applied to the patient-side robot by the patient, can be mapped to their corresponding joint torques τ th ¼ ½τ th1 τ th2 : : : τ thn T and τ p ¼ ½τ p1 τ p2 : : : τ pn T by

Therapist-Side and Patient-Side Robot Dynamics
Both the therapist-side and the patient-side robots can be modeled with the following motion equations. The therapist-side robot dynamics can be modeled by where and M m denotes the inertia matrix for the therapist-side robot, C m denotes the Coriolis and centrifugal matrix, G m denotes the therapist-side robot's gravity vector, and f rm denotes the joint friction torque vector. τ th is preceded by a negative sign because in this article we decide to consider forces/torques applied by the robot to be positive, and that the forces/torques applied by the environment to the robot to be negative. This choice is arbitrary. The patient-side robot dynamics can be modeled by and M s denotes the inertia matrix for the patient-side robot, C s denotes the patient-side robot's Coriolis and centrifugal matrix, G s denotes the patient-side robot's gravity vector, and f rs denotes the joint friction torque vector.

Telerehabilitation System Controller
A DFR architecture will be used for control of the bilateral telerehabilitation system. In the DFR architecture, the patient-side robot follows the position of the therapist-side robot, whereas the therapist-side robot displays to the human operator the interaction forces measured by a force sensor at the patient-side robot's end-effector. This architecture achieves perfect force tracking and assumes that the environment is modeled by a linear spring. For further details of this teleoperation control method, readers can refer to ref. [19]. A detailed schematic of the bilateral teleoperation system is shown in Figure 2. Note that uppercase letters are used to denote the Laplace transforms of the corresponding time-domain position, velocity, force, and torque variables. The matrix Q is introduced to transform the therapist-side robot joint angles into reference joint angles for the corresponding patient-side joints. K s refers to the proportional-derivative position controller for the patient-side robot Note that the inputs to the controller matrix K s are velocities, not positions, as shown later in Equation (18).
K f refers to the force feedback gain matrix for the therapistside robot Z m and Z s denote the linear impedance matrices in the joint domain of the therapist-side and patient-side robot that we approximate from the nonlinear robot dynamics. Z th and Z p refer to the impedance matrices in the Cartesian domain of the therapist's arm and the patient's arm, respectively. As shown in Figure 2, the teleoperation system is divided into five subsystems: therapist, therapist-side robot, control and communication, patient-side robot, and patient. For the therapist and patient, we have where F * th and F * p denote the therapist's and patient's input forces (generated by the muscles with command sent from the central nervous system). For the therapist-side robot and the patient-side robot we have As the therapist and patient work in Cartesian space but the robots work in joint spaces, Jacobian matrices are needed as interfaces between the therapist and therapist-side robot, as well as between the patient and the patient-side robot, to convert Cartesian forces to joint domain torques (based on Equation (5) and (6)) and to convert joint velocities to Cartesian velocities (based on Equation (1) and (2)). As for the controller, we have for the force feedback control on the therapist-side robot, and for the patient-side robot's position control. By combining Equation (15)-(18), the overall system dynamics in the frequency domain can be derived Equation (19) and (20) can be manipulated into the following 2-port network hybrid matrix form where www.advancedsciencenews.com www.advintellsyst.com and I nÂn is the n Â n identity matrix. The 2-port network hybrid matrix will be used later in Section 4 to help derive the desired impedance matrices for the impedance controller.

Identification of Human Arm Impedance
In this section, we first have an overview of the human arm impedance identification techniques used in the literature. Then, based on the nature of the therapy tasks and the requirement to not disrupt the normal flow of teleoperation during the interactive phase, we present the human arm impedance strategy that this work utilizes.

The Proposed Human Arm Impedance Strategy
The approach seen in ref. [18] used very brief force perturbations to minimize the impact on task execution. Such a scheme will be inconvenient if used in our interactive phase as the motions introduced on the therapist's arm following such disturbances may confuse the patient on the other side of the teleoperation and hamper the execution of the cooperative task. The strategy used in this article is similar to this approach with the difference that no extra force perturbations are added. Essentially, the very input forces imported by the patient that occur naturally during the telerehabilitation task are regarded as force disturbances, and a relatively short duration of data is used to identify a secondorder passive impedance model for the human arm. A challenge we face is that, unlike [18] , we do not have an a priori force perturbation sequence and therefore determining the onset of force perturbations is a challenge. As a work-around, let us first analyze the energy absorbed by the therapist's arm during a therapy task over timespan T in one Cartesian direction, where the scalar f th denotes the force/torque applied by the therapist's arm on the therapistside robot in the Cartesian direction under consideration, and the scalar p m denotes the position of the therapist-side robot's end-effector position measured in that Cartesian direction (therefore also the therapist's arm position because the therapist holds onto the therapist-side robot's end-effector). Although the therapy task may involve movements in multiple Cartesian directions, we will consider only one direction for the time being, and duplicate the process to all of the Cartesian directions involved. As shown in Figure 3a for a typical therapy task in one specific Cartesian direction, each rising edge of the absorbed energy corresponds to a task-relevant force onset in that direction. The energy absorption provides us with a good criterion to determine the onset of force perturbations, as it clearly distinguishes when the arm is moving voluntarily (during initial adjustments of the arm position, for example) from when the arm is moving involuntarily (when the arm is knocked away by the force perturbation). Therefore, the onset of a force perturbation in a certain Cartesian direction is determined by the moment whenĖ absorbed (which is the time derivative of the absorbed power, E absorbed in Equation (26)) associated with that Cartesian direction becomes positive (i.e., when energy absorption by the arm begins). A typical perturbation onset determination result is shown in Figure 3b. As can be seen, it would be difficult to determine the perturbation from the force signal alone as it shows no clear distinction between voluntary movement and involuntary movement.
After the perturbation onset is determined to be at time t p for each perturbation for one of the Cartesian DOFs involved, we will try to determine the following impedance model for that Cartesian DOF with respect to force and position data in the time window ½t p , t w M ip ðtÞ þ B iṗ ðtÞ þ K i pðtÞ ¼ Àf ðtÞ, t ∈ ½t p , t w (27) where pðtÞ ¼ p m ðtÞ À p m ðt p Þ, f ðtÞ ¼ F th ðtÞ À F th ðt p Þ, and t w is a selected time window end. The subscript i for the massspring-damper model coefficients indicates that the model is identified in the ith DOF of the Cartesian space. Note also that we subtract the force and position/velocity/acceleration readings at time t p to consider only the force and position/velocity/acceleration changes arising from the disturbance. Classical linear least-squares regression will be used to identify the 1-DOF impedance model in Equation (27)   Determination of t w in Equation (27) depends on the desired duration of the data window T w as t w ¼ t p þ T w . In ref. [18], T w was chosen to be 200 ms as a compromise between the need to use as little data as possible (to accommodate the 100-150 ms window in which the human cannot react voluntarily to the abrupt motion and therefore the arm impedance does not change [20] ) and the model identification calculation that requires sufficient data points. In that work, 200 ms was found to provide all-positive impedance values (corresponding to positive M i , B i , K i in Equation (27)) for over 90% of the perturbations. The need to identify all-positive impedance parameters for the arm in each involved Cartesian DOF is also present in our case because we desire to only obtain the passive impedance of the human arm, as later during the standalone phase, implementing an active impedance may endanger the patient's safety. Furthermore, due to the uncertain nature of the patient-applied (rather than robot-generated) perturbation, we determine a time window length while taking into consideration the goodness of model fit that we measure using the variance accounted for (VAF) test statistic for each perturbation from t p to t p þ T w , with T w incrementing from 100 to 1000 ms at steps of 10 ms Here,f th ðtÞ is the force applied to the therapist-side robot by the therapist's arm as predicted by the identified arm impedance model for a given position trajectory, and F th ðtÞ is the actual force applied to the therapist-side robot by the therapist. The upper bound for T w search is chosen to be 1s as a compromise between the need for a longer window length for impedance model identification and a shorter window for passivity concerns.
Finally, by combining the need to identify all-positive parameters for the impedance model in Equation (27), the need to use data from a short period of time immediately after the perturbation onset, and the need for a good model fit, the strategy used to determine T w associated with the arm impedance model on each involved Cartesian DOF in our work is formulated as follows: for one interactive session, for each involved Cartesian DOF, find the minimum T w ∈ ½100 ms, 1000 ms at 10 ms steps such that the number of all positive-valued identified impedances, Nb pos is maximized subject to the constraint that the average VAF value for the all-positive impedance identification results should be above 95. The algorithm applied to each involved Cartesian DOF is summarized in the flowchart in Figure 4. The T w that leads to the maximum Nb pos is then retained and the impedance model parameters identified with this T w will be used for the next steps.
After determining the 1-DOF model for each Cartesian DOF involved, the complete n-DOF Cartesian arm impedance model can be put together as where pðtÞ ¼ p th ðtÞ À p th ðt p Þ and f ðtÞ ¼ F th ðtÞ À F th ðt p Þ for the same reasons. The matrices M, B, and K can be written as A decoupled n-DOF impedance model with 3 Â n parameters is used here because otherwise there would have been 3 Â n 2 parameters to be identified for a complete n-DOF impedance model given less than 1s of data (1000 data points since for our system, the control and data recording operate at 1 kHz) www.advancedsciencenews.com www.advintellsyst.com via linear least-squares regression, which is very difficult. Therefore, we have decided to simplify the problem by adopting the approach used in ref. [18] in which decoupled models are used. It is important to note that this simplification means the learned impedance is fixed along only the coordinate axes of the reference frame and more complex demonstrations of impedance (where the off-diagonal terms would be nonzero) will not be learned. In addition, the impedance parameters are treated as time-invariant, although for more complicated upper-limb interactions this may not be the case. In light of these caveats, the proposed approach is appropriate for learning simpler interactions, such as those to be discussed in Section 5.

Impedance Control with Disturbance Observer
In this section, we first discuss the impedance controller implemented on the patient-side robot during the standalone phase, once we have the desired impedance parameters M d , B d , and K d . Then we discuss the derivation of the desired impedance to be implemented in the standalone phase.

Impedance Controller Design
Impedance control of a robot can be achieved with model-based approaches or model-free approaches. Model-based approaches such as the one introduced in ref. [21] require the exact knowledge of robot dynamics including joint friction. In our application with an industrial robot in the case studies (a 7-DOF SIA5F robot from Yaskawa Motoman, Miamisburg, OH, USA, as shown in Figure 5b as the patient-side robot and a 2-DOF planar rehabilitation robot from Quanser, Inc., Markham, Ontario, Canada, as shown in Figure 5c as the therapist-side user interface), whereas most dynamic terms can be calculated or estimated, the joint friction is hard to obtain. Unlike the therapistside haptic device which is designed to have low friction, there is considerable friction in the patient-side robot joint that concerns us. The simple yet widely used coulomb þ viscous friction model performed poorly because, first, the linear model cannot capture the highly nonlinear friction in reality. [22] The model's dependency on velocity also makes static friction compensation ineffective, and in our system, the stiction is very large (in the order of 9 Nm in the joint space). Elaborate nonlinear models such as the ones mentioned in refs. [23] or [24] can capture various nonlinear phenomena such as the Stribeck curve, stiction, and presliding displacement, but it is difficult to identify or tune their parameters due to the complexity. In ref. [25], a simpler approach is presented to determine a LuGre model, which is a "dynamic model of friction that describes the behavior of elastic bristles as a function of relative velocity between two surfaces." However, the approach remains largely empirical and requires very refined encoder resolution (especially for measuring the presliding state). Also, the dynamic nature of the LuGre model makes it a potential source of numerical instability in real-time friction compensation implementations. Given the considerable difficulties associated with determining the robot joint friction precisely, as well as the uncertainties in estimating robot link inertias, we decided to implement an online version of the impedance control based on a nonlinear disturbance observer, which has been used to estimate and compensate for the joint frictions and payload of a robotic device in a trajectory following task. [26] For more details of the nonlinear disturbance observer, readers can refer to ref. [26]. The design of the nonlinear disturbance observer-based impedance control is discussed in the following paragraphs.
Consider the general robot dynamics model for the patientside robot shown in Equation (9). By assuming thatM s ðθ s Þ andN s ðθ s ,θ s Þ are the estimates of the actual M s ðθ s Þ and N s ðθ s ,θ s Þ, and that ΔM s and ΔN s are the associated uncertainties N s ðθ s ,θ s Þ ¼N s ðθ s ,θ s Þ þ ΔN (34) Figure 5. a) Direction manipulation: the screwdriving task without using any robots. b) Configuration of the patient-side (patient-side) Yaskawa Motoman SIA5F robot. c) Configuration of the therapist-side (therapist side-side) rehabilitation robot.
www.advancedsciencenews.com www.advintellsyst.com we can define the lumped disturbance vector τ d as in which dynamic uncertainties and joint frictions are lumped into one single disturbance vector. From Equation (9), the patient-side robot dynamics can be rewritten intô For the sake of brevity, the equations describing the nonlinear disturbance observer can be found in the Appendix. After getting an estimation of the τ d term in Equation (36) with the disturbance observer, we can implement a classic impedance control approach in ref. [21] to achieve the following desired dynamics on the patient-side robot where M d , B d , and K d are the desired impedance matrices and p ̤ sd ,ṗ sd , and p sd denote the desired robot acceleration, velocity, and position vectors.
Consider the control law where a q is the reference joint acceleration. By combining Equation (6), (36), and (38), we have θ ̤ s ¼ a q . Let a p be the reference acceleration in Cartesian space, and let or equivalently Together with the kinematic relationship linking joint space acceleration with Cartesian space acceleration we have from Equation (37) and (39) thaẗ Furthermore, take Equation (42) into (40), and then combining with Equation (38) and the disturbance observer (Equation (A.1)) we have Next, we discuss how the desired impedance is found.

Derivation of Desired Impedance
The desired impedance matrices M d , B d , and K d come from the identified therapist arm impedance Z th as distorted by the teleoperation system dynamics. This is because during the interactive phase, it is the therapist arm impedance after the distortion introduced by the teleoperation system that completes cooperative tasks with the patient. From Equation (21), we can derive that the impedance Z d displayed to the patient via the patient-side robot during the interactive phase defined by where Z th is expressed using the identified M, B, and K matrices identified in Section 3 as and the H ij are defined from Equation (22)- (25). With respect to Figure 2, this equates to removing the control and communications channel block and all of the therapist and therapist-side robot components, and replacing them with Equation (43) with inputs from the observer estimateτ d and the robot joint velocity dataθ s .

Case Studies
In this section, the proposed LAR telerehabilitation paradigm is applied to two tasks: a 1-DOF task involving a cooperative screwdriver task and a 2-DOF task involving a collaborative peg-in-thehole task. The 2-DOF task is further supported by including a simulated patient in the trials. We select the screw-driving task because it represents an ADL with synergistic movements similar to opening a doorknob, another common ADL. Peg-in-thehole tasks have seen use in rehabilitation assessment regimes, such as the Nine-Hole Peg Test. [27,28]

Therapy Task Description and Robot Configuration
Consider a task in which a screw is driven by the patient into a surface held in position by the therapist as shown in Figure 5a. For this screwdriver task to be done in the proposed telerehabilitation context, the patient will be tasked to drive the woodscrew into the wooden plate attached to the end-effector of a patientside robot (7-DOF SIA5F robot from Yaskawa Motoman, Miamisburg, OH, USA, as shown in Figure 5b) that is teleoperated by the therapist from a therapist-side user interface (2-DOF planar rehabilitation robot from Quanser, Inc., Markham, Ontario, Canada, as shown in Figure 5c). Thus, in the interactive phase, while the patient uses a screwdriver to drive the screw into the wooden plate fixed to the patient-side robot, the therapist firmly holds the therapist-side haptic device in position.
Although both the therapist-side and patient-side are multi-DOF robots, they have been configured to accommodate the aforementioned task, which naturally involves only a 1-DOF motion in the Cartesian space. Consider motor 0 of the therapistside robot and the sixth joint of the patient-side robot corresponding to the joint angles θ m and θ s in Figure 5, respectively. The patient-side position θ s is made to follow the therapist-side position θ m , whereas interaction forces at the patient-side side are reflected to the therapist-side side. Because of this 1 to 1 correspondence, Q in Figure 2 becomes Q ¼ 1. The second joint of the therapist-side robot is passively (physically) clamped in its home position with the corresponding motor (motor 1) turned off, whereas the other six joints of the patient-side robot are actively held in position via high-gain PID position control. In Figure 5, both robots are at their home positions (θ m ¼ θ s ¼ 0). Cartesian frames are attached to the end-effectors of the therapist-side and the patient-side, as shown in Figure 5. Note that for small θ m and θ s , the motions of the two robots can be approximated to be along a Cartesian axis. The task is therefore in the Y direction. Through teleoperation, the therapist tries to resist the pushing forces of the patient applied in the Y direction by displaying a stiff impedance to the therapist-side robot in that direction, so that the patient can complete the screwdriver task. The link lengths for the therapist-side and patient-side robot are l m and l s , respectively, measuring from the rotating axes to the centre of the handle for the therapist-side robot (Figure 5c) and to the screw location on the wooden plate for the patientside robot (Figure 5b). As a result, the Jacobian matrices in Equation (1) and (2) have become scalars and are respectively. p m and p s are therefore measured in the Y direction. 1-DOF interaction forces applied onto the robots by the therapist or the patient: f th and f p in Equation (5) and (6) are directly measured with two ATI Gamma NET force/torque transducers (Apex, NC, USA) attached to the end-effectors of the two robots at 1 kHz sampling rate. The communication channel is implemented using the Winsock application programming interface over the Ethernet using the UDP protocol at 1 kHz sampling rate-the same rate as the one used in the robot control loops (for reading encoders and issuing torque commands) of both robots. Furthermore, teleoperation system parameters Z m in Equation (19) can be modeled as an inertia, Z m ¼ M m s, as we choose to ignore the friction associated with the actuated link of this haptic device. The patient-side robot is modeled as an inertia and a damper because it has significant damping and friction needing to be modeled, leading to Z s ¼ M s s þ B s in Equation (20).

Experimental Results
During the interactive phase, we use a feedback gain of K f ¼ 0.5 for DFR teleoperation because of the following two reasons: 1) it provides the therapist with a good perception of the perturbations from the patient side; 2) the force feedback is not too strong for the therapist to handle with ease. Numerical values for other teleoperation system parameters are either measured or chosen, as shown in Table 1. Note that the values for M m , M s , and B s are obtained using system identification method similar to what was done in ref. [15]. l m and l s are obtained by directly measuring the concerned robot links. K f , K v , and K s are directly specified in the controller software. A total of three interactive sessions were carried out, each lasting around 60 s containing 20-26 perturbations. The experimental trials were carried out during the same day with one healthy person (the first author of the article) acting as the therapist and another healthy person (the second author of the article) acting as the patient. The arm impedance identification results are shown in Table 2, including M, B, K, the percentage of all-positive identification results Nb pos over the total number of considered perturbations in each session, and the chosen data window length T w for model identification.
As can be seen from Table 2, with our impedance model identification algorithm shown in Figure 4 and 90% of the perturbations yield all-positive identification results with an average data window length of about 600 ms. We use the averaged impedance parameters: M ¼ 1.92 kg, B ¼ 55.54 Ns m À1 , and K ¼ 1454.52 N m À1 to to make up the Z th term defined in Equation (46). Combining with the definition in (45), the displayed impedance during interactive phase Z d can be numerically approximated into the following form for impedance control implementation 14 Ns m À1 , and K d ¼ 3467.36 N m À1 following an optimization operation.
In the standalone phase, the patient-side robot is programmed under the impedance control law specified in Equation (43) with the parameters of M d , B d , and K d found in interactive phase discussed in the previous subsection. For this 1-DOF screwdriving task under consideration, p ̤ sd ,ṗ sd , and p sd in Equation (43) are all 0. For the disturbance observer, we need to first find Y using (Equation (A.8)) from the parameters ζ, β, and σ 2 . As the dynamics are in 1-DOF, M s does not vary over time. Therefore, ζ, representing an upper bound of˙M s , is zero; σ 2 , being an upper bound ofM s , can be taken as the same value www.advancedsciencenews.com www.advintellsyst.com ofM s . However, the value of M s , or the theoretical value of the patient-side robot inertia, provides a good candidate forM s (therefore σ 2 , in practice it is tuned together with β for best observer convergence while being able to reject noise and instability). Through experiments, we takeM s ¼ σ 2 ¼ 0.15 and β ¼ 1050 for the disturbance observer. The patient (imitated by a healthy person) is then able to complete the screwdriver task with the patient-side robot alone in this phase. The impedance control implementation is validated by comparing the actual patient-side robot position against the simulated patient-side robot position based on the desired impedance and the measured interaction forces. Based on the therapist's arm impedance that was retrieved during the interactive phase, patient-side robot positions upon perturbation are simulated. This simulation is the ideal movement, based on the therapist's arm impedance, that we hope the patient-side robot will do in response to measured interaction forces. However, factors such as the patient-side robot's inertia and joint friction may impact the actual performance of the robot in replicating these movements. The graph shows the comparison between the simulated (ideal) movement of the patient-side robot versus its actual movement.
The comparison result is shown in Figure 6 where a good match is shown. Note that the effect of the strong stiction inherent in our industrial manipulator arm is still visible, but it is mitigated (position error no more than 1 mm during the task) due to the disturbance observer approach. We can also see from Figure 6 that the maximum displacement in the Y direction of the patient-side arm is around 12 mm, which is also the case in the interactive phase, showing that the desired behaviour has been successfully mapped from the therapist to the patient-side robot during the standalone phase.

Therapy Task Description and Robot Configuration
The 2-DOF task consists of the patient trying to put an aluminum mechanical part representing a 1D "hole" onto a peg held by the therapist, as shown in Figure 7a. In the proposed LAR paradigm, the patient will insert the hole onto the peg attached to the endeffector of the patient-side robot shown in Figure 7b, which is teleoperated by the therapist from the therapist-side interface, as shown in Figure 7c. While "hole-onto-the-peg" insertion seems to be a more appropriate term to name the considered task, as the peg and the hole play interchangeable roles, we still name the task "peg-in-the-hole" insertion. The peg-in-the-hole insertion task is a challenging manipulation task, involving both position and force control. In our case, the task requires the patient to first align the hole with the peg by wiggling it primarily in the Y direction, as shown in Figure 8a. This procedure is made possible by the curved opening of the hole and the curvature of the peg tip. This step is also shown in Figure 7a in the direction manipulation case, where no robots are involved. Once the hole and the peg are lined up (illustrated in Figure 8b), the patient pushes the hole onto the peg by exerting a force in the X direction (shown in Figure 8). The X and Y Cartesian directions are defined  The peg-in-the-hole task without using any robots. b) Configuration of the patient-side (patient-side) Yaskawa Motoman SIA5F robot as the patient-side. c) Configuration of the therapist-side (therapist-side) rehabilitation robot as the therapist-side.
www.advancedsciencenews.com www.advintellsyst.com for the patient-side robot in Figure 7b. Note that the X-Y frame origin is placed at the centre of the peg tool when the patient-side robot is at home position, and it does not move with the peg tool. We name this frame the patient-side base frame. For the therapist-side robot, X and Y Cartesian directions are the same as the patient-side robot and the frame origin is attached to the handle when the therapist-side robot is at its home position.
Although the therapist-side frame origin is also defined using the end-effector, it does not move with the handle either once it is defined. We name this X-Y frame related to the therapistside robot, the therapist-side base frame. Note that the therapistside robot is not at its home position in Figure 7c. The therapist-side robot is at its home position when θ 1m and θ 2m shown in Figure 7c are both zero. Axes parallel to the X-and Y-axis of the therapist-side robot but attached to the robot centre axis shown as X 0 and Y 0 in Figure 7c, which will be used later for details on the joint angle definition for the therapistside robot. Similar to the 1-DOF task screwdriver task, the therapist holds firmly the therapist-side haptic device in both X and Y directions, whereas the patient tries to complete the task during the interactive phase. We make this choice because it is desirable for the patient to learn to behave compliantly. To complete the peg-in-the-hole task, both operators cannot be both compliant or both rigid at the same time (or else the hole risks not being able to be lined up with and inserted into the peg).
In this task, we face the difficulty of teleoperating two robots with different kinematics, workspaces and DOFs, because of our choice of the therapist-side and patient-side robots. We solve this problem by making the robot with more DOFs, and a larger workspace has a similar effective geometry as the robot with fewer DOFs and a smaller workspace. To do so, we make the 7-DOF patient-side robot take on an effective 2-DOF geometry. The first, fourth, and sixth joints of the patient-side robot are arranged according to inverse kinematics such that the distance between the first joint axis and the sixth joint axis, which is l 1s in Figure 7b, takes on a desired value and becomes the length of the first "effective link" of the patient-side robot.
During teleoperation, high gain PID control is applied to joints 2-5 of the patient-side robot to maintain this geometry. The second "effective link" of the patient-side robot is the same as the actuated link presented in the 1-DOF case study (comprising the last link of the patient-side robot, the force sensor and peg tool) with link length l 2s , as shown in Figure 7b. For the therapistside robot, both robot joints are actuated and the link lengths are l 1m and l 2m , respectively, as shown in Figure 7c. Note that l 2m is the same as l m for the 1-DOF screwdriver task. We define that at home position, the first effective link and the second effective link of the patient-side robot are perpendicular to each other (shown in Figure 7b), whereas the second effective link points to the positive X direction. The same can be said about the therapist-side robot: at home position, the first link and the second link are perpendicular to one another, whereas the second link points to the positive X direction.
In terms of joint angles, we use the first and sixth joint angles of the patient-side robot minus their values when the patient-side robot is at its home position as the new joint angles: θ 1s and θ 2s . Therefore, at the home position shown in Figure 7b, θ 1s ¼ θ 2s ¼ 0. For the therapist-side robot, we define the first joint angle to be the angle formed by rotating the Y 0 axis counterclockwise around the robot centre joint to be parallel to the first robot link (l 1m ), and the second joint angle to be the angle formed by rotating the X 0 axis counterclockwise to be parallel to the second robot link (l 2m ). Note that in the configuration shown in Figure 7c, θ 1m has a negative value and θ 2m has a positive value. The joint angle configurations for the therapist-side and the patient-side are shown in Figure 9.   We made l 1s ¼ l 1m ¼ 0.254 m. As we saw in the 1-DOF task, l 2s ¼ l 2m ¼ 0.2667m (listed in Table 1). Therefore, we have made the therapist-side and patient-side robots have similar effective geometries. During teleoperation, we simply make the patientside robot follow the therapist-side robot by using the current values of θ 1m and θ 2m À θ 1m as the reference positions for θ 1s and θ 2s . This feature can be seen from from Figure 9, by defining the angles between the X-axis and the second link of both robots (l 2m and l 2s ) as q m and q s , respectively. As q m and q s should be equal to each other, we have Therefore, As θ 1s takes θ 1m as its reference, we can replace θ 1s in Equation (52) by θ 1m and obtain as the reference angle position for θ 2s . Accordingly, matrix Q in Figure 1 take the following form To recapitulate mathematically the above description of the robot geometries, if we express Cartesian positions p m and p s as ½x m y m T and ½x s y s , respectively, where x and y represent the Cartesian coordinate on the X-and Y-axis, then p m and p s can be expressed as Àl 1m sinðθ 1m Þ þ l 2m cosðθ 2m Þ À l 2m l 1m cosðθ 1m Þ þ l 2m sinðθ 2m Þ À l 1m # (55) Àl 1s sinðθ 1s Þ þ l 2s cosðθ 1s þ θ 2s Þ À l 2s l 1s cosðθ 1s Þ þ l 2s sinðθ 1s þ θ 2s Þ À l 1s # Accordingly, Jacobian matrices J m and J s in Equation (1) and (2) become Àl 1m cosðθ 1m Þ Àl 2m sinðθ 2m Þ Àl 1m sinðθ 1m Þ l 2m cosðθ 2m Þ # (57) Àl 1s cosðθ 1s Þ À l 2s sinðθ 1s þ θ 2s Þ Àl 2s sinðθ 1s þ θ 2s Þ Àl 1s sinðθ 1s Þ þ l 2s cosðθ 1s þ θ 2s Þ l 2s cosðθ 1s þ θ 2s Þ # (58)

Experimental Results
During the interactive phase, a series of 11 peg-in-the-hole insertion tasks was completed via teleoperation. The experimental trials were carried out during the same day with one healthy person (first author of the article) acting as the therapist and another healthy person (second author of the article) acting as the patient. The force feedback gains for the therapist-side robot in Equation (12) and the position control gains for the patient-side robot in Equation (11) were chosen to be While the tasks were being completed, the position and force data on the therapist's side were recorded. The T w search algorithm shown in Figure 4 yields 960 ms for the X direction and 170 ms for the Y direction. By using the technique described in Section 3, we obtain the following identification results for the therapist's arm impedance in Equation (30) where the subscripts 1 and 2 represent the X and Y directions, respectively.
We can see that the arm demonstrates a stiffer impedance along the X direction, which agrees with the findings of ref. [15]. We can write Z th as Z th ¼ Now we need to find the the impedance displayed to the patient during the interactive phase Z d in Equation (45). By applying a persistent perturbation to the therapist-side robot, we can identify the following linearized form of Equation ( In the standalone phase, the patient-side robot is programmed to demonstrate the desired impedance (Equation (37)) via the controller (Equation (43)). With trial-and-error, the best disturbance observer was obtained with the following parameters for (Equation (A.8)): The patient (imitated by a healthy person) is then able to complete the peg-in-the-hole task with the patient-side robot alone in this phase. We repeated the peg-in-the-hole insertion task in the absence of the therapist for several of times. The task was completed successfully with the patient-side robot alone. If the robot impedance was not controlled properly, the robot could yield to pressure from the patient and the task would fail. The impedance control implementation is validated by comparing the actual robot position in both X and Y directions with what the target impedance model predicts based on the measured forces. A good match in both directions is shown in Figure 10. The desired behaviour has therefore been successfully mapped from the therapist to the patient-side robot during the standalone phase for the 2-DOF peg-in-the-hole insertion task.

2-DOF Task with Simulated Patient
Following the same procedure as the previous peg-in-the-hole trial and using the same gains, two healthy persons (different people compared from the previous trial) took the role of the therapist and the patient. However, this time we simulated a patient with arm tremors using a Transcutaneous Electrical Nerve Stimulation (TENS) device. Two electrodes were attached to the acting patient's forearm and were set to provide stimulation at roughly 3 Hz. The stimulation of the upper arm and wrist at low frequency induces a behavior similar to poststroke patients. [29] The T w search algorithm yields 170 ms for the X direction and 240 ms for the Y direction. The therapist's arm impedance is It can be seen that the arm is stiffer in the X direction but is less stiff when compared to the therapist impedance from the previous trial. This feature can be attributed to a different person acting as the therapist.
As the Z m and Z s parameters come from robot identification, they remain unchanged. After the optimization approach, the following desired impedance parameters are determined In the standalone phase, the disturbance observer parameters are found using trial-and-error ζ ¼ 3 Nm s rad À1 β ¼ 50 Nm s rad À1 (69) As shown in Figure 11, the Motoman Robot, put under impedance control, can match the movement of a simulated patient www.advancedsciencenews.com www.advintellsyst.com with tremors in the standalone phase. The TENS device produced wrist flexion and extension, resulting in more movement in the Y direction.

Conclusions and Future Work
In this study, we proposed and demonstrated a novel telerehabilitation approach: the LAR over two collaborative tasks: a screwdriver task and a peg-in-the-hole insertion task. During the interactive phase, the therapist supported the patient in completing the task and simultaneously the impedance of his arm was measured by the therapist-side haptic device. It is worth noting that the interactive phase should last long enough in time to allow a sufficient number of task repetitions to take place for impedance model identification. In the case studies of this article, we found that a few minutes' of interaction provided enough data for impedance identification. The measured impedance was then processed taking into account the teleoperation system dynamics to obtain the desired impedance parameters used for impedance control implementation during the standalone phase. In the standalone phase, the impedance control was successfully implemented on the patient-side robot and the therapist's role in completing the tasks was successfully replicated by the patient-side robot. We showed the feasibility of the proposed LAR telerehabilitation paradigm and its potential in time-sharing a therapist and thus partly automating the rehabilitation therapy. One important caveat regarding the formulation of the learned impedance behavior used in this work is that the impedance is regarded as time-invariant and fixed along the coordinate axes. We acknowledge this as a potential limitation to our approach, given the knowledge that therapeutic assistance is often provided in a time-varying and direction-dependent manner with respect to the task. Another consideration is to explore different methods of selecting the desired impedance matrices, instead of only considering the mean of the fitted values which is prone to outliers. This could include taking the median of the values, or fitting the parameters in a coordinate space that better encompasses the behavior of the responses such as the natural frequency of the impedance.
In the future, new multi-DOF tasks that are closely related to current rehabilitation practices for upper-limb rehabilitation can be designed. The tasks need to be closely related to ADL and therapist evaluation of the paradigm will be focused on. Doing so will likely require reformulating the method of learning impedance to consider tasks as dynamic and trajectory-based. Clinical patient-oriented studies will be carried out to study the usefulness of the proposed paradigm in real clinical settings.