Multi-Domain Dynamic Modelling of a Low-Cost Upper Limb Rehabilitation Robot

Tracking patient progress through a course of robotic tele-rehabilitation requires constant position data logging and comparison, alongside periodic testing with no powered assistance. The test data must be compared with previous test attempts and an ideal baseline, for which a good understanding of the dynamics of the robot is required. The traditional dynamic modelling techniques for serial chain robotics, which involve forming and solving equations of motion, do not adequately describe the multi-domain phenomena that affect the movement of the rehabilitation robot. In this study, a multi-domain dynamic model for an upper limb rehabilitation robot is described. The model, built using a combination of MATLAB, SimScape, and SimScape Multibody, comprises the mechanical electro-mechanical and control domains. The performance of the model was validated against the performance of the robot when unloaded and when loaded with a human arm proxy. It is shown that this combination of software is appropriate for building a dynamic model of the robot and provides advantages over the traditional modelling approach. It is demonstrated that the responses of the model match the responses of the robot with acceptable accuracy, though the inability to model backlash was a limitation.


Introduction
There is increasing interest in the use of robots to assist the rehabilitation of stroke patients, partly due to the repetitive nature of post-stroke rehabilitation and the integration of stimulating video games [1] to increase patient motivation and partly due to a projected increase in strokes due to an ageing population, causing a strain on rehabilitation services [2]. Stroke can have severe consequences for the patient's quality of life, with reports estimating that up to 20% of patients require intensive rehabilitation programs [3].
Rehabilitation robots may be broadly classified into two categories: Class 1 or Class 2 [4,5]. Class 1 robots are high cost and are designed for clinical use with a therapist. Class 2 robots are low cost and intended for home use. Rehabilitation robots that are commercially available are Class 1 robots, and therefore high cost. Examples include the InMOTION (previously known as MIT Manus [6,7]) and the Armeo Power (previously known as the ARMin [8,9]). However, a large study has recently shown that the use of expensive robotics in a clinical setting provides little additional benefit to traditional rehabilitation alone and is not cost effective [10]. Essentially, if there is a choice between traditional physiotherapist-led poststroke rehabilitation or robot-based rehabilitation, there is little benefit to selecting robot-based rehabilitation. Both options require travel to a clinical setting and face-to-face contact with a physiotherapist.
As a response to the COVID-19 pandemic, many countries imposed national lockdowns, with non-essential travel banned and social mixing criminalised. Health services were put under considerable strain dealing with the effects of COVID-19, with many services limited or cancelled. There was a noted decrease in stroke patient admissions, particularly patients presenting with minor symptoms and transient ischaemic attack (TIA) in Italy, France, Germany [11], Canada [12], and Norway [13]. In the U.K., there was a 12% reduction in the number of admissions for stroke patients between 1 October 2019 and 30 April 2020 compared to the same period in the three previous years [14]. It was noted that quality of care did not decrease for admitted stroke patients, but there are no data on access to rehabilitation or long-term outcomes. It is postulated that hospital avoidance is the likely cause of the decrease in stroke patients rather than a decrease in the number of strokes.
There is little doubt that social distancing has had an effect on post-stroke rehabilitation. Indeed, in Canada, it was noted that access to rehabilitation care has been significantly reduced [12]. According to the Stroke Association, around half of stroke survivors in the U.K. have had therapy appointments cancelled or postponed [15]. These cancellations may have occurred due to an increase in protective measures, which require more turnaround time between patients, a reduced capacity in rehabilitation centres due to social distancing requirements, and staff absence due to self-isolating or infection. Further to this, 56% of stroke patients have not felt safe to attend scheduled appointments. This is likely due to fear of becoming infected with COVID-19, especially considering that stroke patients are at higher risk [16].
The COVID-19 pandemic has created a need for social distancing and a new paradigm of hospital avoidance due to the fear of infection. It was noted by [12] that telerehabilitation is an effective and well-accepted method for providing access to therapy, though it was considered that this requires family members and care givers to be given additional training. Virtual reality (VR) technologies and existing computer game systems such as the Nintendo Wii could be used to supplement rehabilitation [17], though it was noted that most homebased exercises require oversight from a physiotherapist. It is clear from the literature that there is a need for Class 2 rehabilitation robotics in the current climate of social distancing.

MyPAM Rehabilitation Robot
MyPAM, shown by Figure 1, is a low-cost planar robot designed for upper limb rehabilitation of stroke patients in the home. The robot workspace exists in the x-y plane, with gravity acting in the z-plane. The robot is a 2 DoF powered arm, which assists the patient to reach targets presented to them on a computer screen. It is important to monitor patient progress during the course of rehabilitation to ensure that the recovery program is appropriate. Tracking patient progress through a course of tele-rehabilitation requires continuous logging of the robot position data during exercise and periodic testing with no robotic assistance. These test data must be compared with previous test attempts and an ideal baseline, which requires a good understanding of the dynamics of the robot when connected to a human arm. Joint 0 is driven using a MAXON RE40 DC motor and has gear reduction of 40:1, composed of a planetary gear set attached to the motor with a gear reduction of 15:1 driving a spur gear pair with a gear reduction of 8:3. Joint 1 is driven using a MAXON DCX32L 24V DC motor and has a gear reduction of 70:1, composed of a planetary gear set attached to the motor with a gear reduction of 35:1 driving a bevel gear pair with a gear reduction of 2:1. The motors are powered using bespoke motor control boards, which handle current protection and disconnect the motors from power when the motor demand is 0V. Disconnecting the motors in this way prevents the braking effect caused by back EMF and allows the back-drivability of the robot, which is necessary for rehabilitation robotics. Joint position feedback is provided by quadrature encoders. MyPAM is controlled using a National Instruments MyRIO programmed using LabVIEW.
Whilst most upper limb rehabilitation robots use an admittance control or an impedance control approach, the high cost of force sensors required makes admittance control or impedance control inappropriate for a low-cost Class 2 robot such as MyPAM. As such, MyPAM uses position control as the control paradigm. A previous iteration of MyPAM (known as HCARR [5]) using a similar position control approach was tested in a clinical trial. In this trial, seventeen participants used the HCAAR for eight weeks, and statistically significant improvements were observed in kinematic and clinical outcomes.

Dynamic Modelling
The traditional approach to modelling involves forming dynamic equations (equations of motion), typically using Newtonian or Lagrangian formulations [18]. This becomes cumbersome and error prone as more components and more degrees of freedom are modelled. There has been a growing trend towards simulation-based dynamic modelling due to the increased capability of physics-based simulation software and the increased availability of computing power. Furthermore, efficient recursive algorithms have been developed, increasing the performance of simulation software [19].
The model presented in this paper incorporates three domains (mechanical, electromechanical, and control) and is constructed using a combination of software packages from MATLAB. The mechanical components were modelled using SimScape Multibody, which allows three-dimensional dynamic modelling and simulation of mechanical systems. The motors in each of the powered joints of MyPAM were modelled using SimScape, which allows modelling and simulation of electrical and electro-mechanical systems. Simulink and MATLAB were used to apply control to the model. SimScape and SimScape Multibody utilise a physical network approach, whereby block connections communicate information about power [20] and are presented within the Simulink environment. The SimScape Multibody MultiPhysics library was used to provide an interface between SimScape and SimScape Multibody, such that the motor models may be used to drive the joints.

Human Arm Proxy
A human arm proxy capable of replicating the seven degrees of freedom of a human arm was used for this work. The human arm proxy was previously developed for the validation of a different upper limb rehabilitation robot [21]. The human arm proxy emulates posterior/anterior translation and superior/inferior translation at the shoulder using a flexible steel rod secured at a distance away from the shoulder joint. Extension/flexion of the upper arm is achieved using a rotational joint at the shoulder. Abduction/adduction is achieved using hinge at the shoulder. Extension/flexion of the lower arm is achieved by a hinge at the elbow, with the addition of a spring contributing to the stiffness of the elbow joint. External/internal rotation and pronation/supination are achieved with rotational joints, with stiffness of these rotations facilitated by the addition of friction clutches. The human arm proxy is shown by Figure 2, with four of the seven degrees of freedom illustrated in the front view and the remaining three degrees of freedom shown in the lateral view.

Aim
The aim of this study is to create a multi-domain dynamic model capable of characterising movement to better understand patient progress.

Computational Modelling and Simulation-SimScape Multibody and the Mechanical Domain
The mechanical components of MyPAM were identified and replicated in SimScape Multibody to create a three-dimensional dynamic model in the mechanical domain. Joint 0 was constrained to rotate about the global origin. Similarly, a 3D dynamic model was created for the human arm proxy. The shoulder of the human arm proxy was constrained to rotate about an appropriate point in the global workspace, and the hand was connected to the MyPAM model with a 6-DoF rotational joint. The modelling procedure for each model followed the workflow:

1.
Define the global reference point, coordinate system, and simulation settings by placing the world frame block in parallel with the solver configuration block and the mechanism configuration block; 2.
In turn, place and configure a solid body block for a component or place and configure a joint; 3.
Connect the block input/output ports, ensuring that rigid body transforms are used where appropriate to translate or rotate frames.
All solid bodies were configured with dimensions and assigned a material density, from which the inertia and mass were automatically calculated. Coordinate frames were assigned at suitable locations for each solid body, which created the input/output connection ports required for components and joints to be connected. Joints 0 and 1 of the MyPAM model were not assigned internal mechanics (friction) because these were instead applied to the motor model for each joint. The two joints were configured to receive a torque input, which was provided by the motor model. All joints in the human arm proxy model were configured with internal mechanics. Finally, Joint 1 and the end-effector of MyPAM were configured to the output position in the global workspace, which was externally logged to MATLAB as the output of each test. The SimScape Multibody model for the unloaded MyPAM and the SimScape Multibody model for MyPAM connected to the human arm proxy are shown by Figure 3.

Computational Modelling and Simulation-SimScape and the Electro-Mechanical Domain
A multi-domain motor model for each powered joint of MyPAM was created in SimScape using a combination of electrical and mechanical blocks. For each powered joint, the motor demand from the control model was provided to a controllable voltage source block. The output of the controllable voltage source block was provided to a DC motor block.
The DC motor block in each motor model was configured with the parameters defined in the datasheet for each joint, which are shown by Table 1. The output of the motor was passed through a gear ratio block configured with the correct gear ratio. A friction block was placed in parallel across the motor and gear blocks to model the friction across the joint. The friction block in each motor model was configured as the friction parameters found in a previous experiment, which are shown by Table 2. Note that the breakaway friction velocity was set to the default value, which is close to zero. Finally, a conversion block from the SimScape Multibody MultiPhysics library was placed so that the torque output of the motor modelled in SimScape may be provided to the relevant joint of MyPAM modelled in SimScape Multibody. The motor model is shown by Figure 4.

Computational Modelling and Simulation-MATLAB, Simulink, and the Control Domain
The control scheme was implemented in MATLAB and Simulink and followed the block diagram shown by Figure 5.  A discretised minimum jerk trajectory between the start and end point of the reaching movement was initially generated using MATLAB to produce a set of Cartesian position vectors for the end-effector, which were used as an input to the whole Simulink model. A MATLAB function block was used in Simulink to perform the inverse kinematics, converting the Cartesian position vectors into joint position demand vectors used for control. Joint position feedback from the MyPAM SimScape Multibody models was used to evaluate the joint position error using standard Simulink mathematics blocks. Simulink PID blocks were used to generate a motor control signal u (in Volts) for each joint using the respective joint position errors. MyPAM uses PI position control only, as stated in Section 1.1. In the unloaded case for both the robot and the simulation, the gains were P = 1, I = 0.01. In the loaded case for both the robot and the simulation, the gains were P = 1.5, I = 0.01.

Testing and Validation Methodology
A series of tests was performed to validate the performance of the model against the performance of MyPAM. The first pair of tests compared the performance of the MyPAM model against the performance of the MyPAM in the unloaded condition, with no external loading applied. In the first of these tests, a desired trajectory was applied only in the x-direction, and in the second test, a desired trajectory was applied only in the y-direction. The second pair of tests compared the performance of the MyPAM model against the performance of the MyPAM in the loaded condition, with the human arm proxy model connected. The same desired trajectories were applied only in the x-direction and only in the y-direction for each test, respectively. A summary of the testing is provided by Table 3. Figure 6 shows MyPAM with the human arm proxy and a simulation of the MyPAM and human arm proxy models.

Unloaded MyPAM: Tests 1 and 2
The graph shown by Figure 7 shows the simulated response and the mean real response of MyPAM when subjected to a minimum jerk trajectory in the x-direction, with the mean x-and y-positions against time shown by Figure 8. It may be observed that both responses show a characteristic curve, caused by the small demand of the motor at Joint 0 and the large demand of the motor at Joint 1. The graph shown by Figure 9 shows the simulated response and the real responses of MyPAM when subjected to a minimum jerk trajectory in the x-direction, with the x-and y-positions across all repeats against time shown by Figure 10.   The graph given by Figure 11 shows the simulated response and the mean real response of MyPAM when subjected to a minimum jerk trajectory in the y-direction, with the mean x-and y-positions against time shown by Figure 12. The graph given by Figure 13 shows the simulated response and the real responses of MyPAM when subjected to a minimum jerk trajectory in the y-direction, with the x-and y-positions across all repeats against time shown by Figure 14.

Loaded MyPAM: Tests 3 and 4
The graph given by Figure 15 shows the simulated response and the mean real response of MyPAM when subjected to an x-direction minimum jerk trajectory, with the mean x-and y-positions against time shown by Figure 16. It may be observed that the curved response in Figure 7 is present. The graph given by Figure 17 shows the simulated response and the real responses of MyPAM when subjected to an x-direction minimum jerk trajectory, with the x-and y-positions across all repeats against time shown by Figure 18.  The graph given by Figure 19 shows the simulated response and the mean real response of MyPAM when subjected to a y-direction minimum jerk trajectory, with the mean x-and y-positions against time shown by Figure 20. The graph given by Figure 21 shows the simulated response and the real responses of MyPAM when subjected to a y-direction minimum jerk trajectory, with the x-and y-positions across all repeats against time shown by Figure 22.

Discussion
In both the unloaded and loaded cases, the x-direction trajectory tracking shown by Figures 7 and 15 (Tests 1 and 3) is much smoother than the y-direction trajectory tracking shown by Figures 11 and 19 (Tests 2 and 4) for both the simulated response and the robot response. This is because in this area of the robot workspace, x-direction movement is achieved mainly by rotation of Joint 1, with very little movement required by Joint 0. As expected, the effects of inertia are less pronounced when the majority of the movement occurs at Joint 1 because the majority of the mass of the robot is located between Joint 0 and Joint 1. Similarly, the y-direction trajectory tracking was poorer because most of the movement is achieved moving Joint 0, where the effects of inertia are more pronounced.
In the unloaded case, the x-direction trajectory tracking (Test 1), seen in Figure 8, shows good x-position agreement between the simulated response and the robot response, with both closely following the trajectory demand. It should be noted, however, that the robot response shows a small amount of overshoot that the simulated response does not show. There is a similar good x-position agreement between the simulated response and the robot response in the loaded case (Test 3) seen in Figure 16, though there is a greater difference between the robot response and the simulated response than is present in the unloaded case. This is likely caused by imperfect modelling of the human arm proxy. The human arm proxy has a greater number of joints than MyPAM, making it considerably more difficult to correctly account for the effects of friction.
The y-position agreement between the robot response and the simulated response for the unloaded x-direction trajectory tracking (Test 1) shown by Figure 8 is satisfactory since both lie within the 10 mm dead-zone allowed between the position of the end-effector and the target, though it is noted that the robot response shows greater overshoot and steady-state error than the simulated response. This may be accounted for by the absence of backlash modelling in the MyPAM model. There is a greater y-position agreement between the robot response and the simulated response for the loaded x-direction trajectory tracking (Test 3) shown by Figure 16, though the simulated response shows a small overshoot at around 1 s that the robot does not. This is again likely caused by imperfect modelling of friction in the human arm proxy, which added less damping to the simulated response than the human arm proxy added to the robot response.
In the unloaded case, the y-direction trajectory tracking (Test 2), seen in Figure 11, shows poorer x-position agreement between the simulated response and the robot response than may be observed in the response to the x-direction trajectory (Test 1). This is likely due to the lack of backlash modelling, which is relatively large at Joint 1. The effect of backlash is pronounced in the robot response because tracking a y-direction trajectory in this area of the robot workspace requires small movements of Joint 1 where the backlash to the x-position demand ratio is high. Despite this, it was observed that the steady-state response in the x-position, seen in Figure 12, is closely matched between the robot response and the simulated response. In the loaded case, the x-position agreement between the simulated response and the robot response while tracking the y-direction trajectory (Test 4), seen in Figure 20, is poorer than was observed in the unloaded case. This is again likely caused by imperfect modelling of friction in the joints of the human arm proxy model.
In the unloaded case, the y-direction trajectory tracking (Test 2), seen in Figure 11, shows reasonable y-position agreement between the simulated response and the robot response, though it was observed that the model tracks the trajectory demand better than the robot. This likely occurred due to the backlash at Joint 1, which caused greater error in the x-position with corresponding inertial effects on the robot. The loaded case (Test 4) shows a better y-position agreement between the simulated response and the robot response, seen in Figure 19, due to the damping effects of the friction in the human arm proxy.
It has been noted that SimScape Multibody was poor at modelling instability. This became apparent when attempting to tune the PID gains for the motor controllers using the Ziegler-Nichols method on the unloaded MyPAM model, where it was not possible to introduce instability even with an extremely large derivative gain. Indeed, the model response with large PID gains followed the trajectory perfectly. This means that the PID tuning strategy required a balance of rough tuning with the models to identify sensible gains, followed by finer tuning and validation on the robot. The set of software cannot be relied upon for tuning the control system.
SimScape Multibody does not provide a function for modelling backlash. This was a problem because the backlash in Joint 1 of MyPAM is relatively large. It would be preferred to model the backlash because it has a significant effect on the dynamics, particularly in circumstances where Joint 0 must make frequent changes in direction to track the trajectory. The effects of the backlash on MyPAM can be clearly seen by the jagged motion of the real MyPAM in Figure 10. In the same figure, the jagged motion, which would be caused by the effects of the backlash, is absent from the response of the simulated MyPAM, which instead shows a smoother response.
The agreement between the robot response and the model response is acceptable in the primary direction of travel in all cases (i.e., for motion primarily in the x-direction, there is good agreement in the x-position against time). The agreement between the robot response and the model response perpendicular to the primary direction of travel is less good, though it is considered acceptable because this particular region of the workspace was selected for testing as it is where the robot has the most difficulty. It is apparent that the performance of the model more closely matches the performance of the robot when loaded than when unloaded, which is a useful outcome from a rehabilitation perspective since the main purpose of the model is to provide a baseline against which patient data may be compared.

Conclusions
The combination of MATLAB, SimScape, and SimScape Multibody provides an appropriate tool for multi-domain modelling of MyPAM. While the response of the simulation is closer to the response of the robot in the loaded case than in the unloaded case, the model has sufficient accuracy to allow the dynamics of the robot to be accounted for when analysing patient movement data. Importantly, it can be seen that the shape and direction of the response curves produced by the models adequately match the shape and direction of the response curves produced by MyPAM. Modelling the robot using this combination of software allowed the creation of a model with a greater deal of both fidelity and com-plexity than traditional mathematical modelling alone would allow, though it was noted that the absence of backlash from the model is a limitation. The failure of the models to produce instability in response to a large control input means that this set of software is not appropriate for tuning the control system.