Vibrations of cable-suspended rehabilitation robots

Abstract

Rehabilitation robots help the treatment of diseases by performing cyclic exercises for a long period of time. These exercises must perform movements of the patient’s limbs; thus, the robots are required to be flexible and safe. Among rehabilitation robots, cable robots are widely used due to their unique properties, such as being lightweight and the possibility of being equipped with magnetic hooks to improve both safety and ease of use. However, the elasticity and flexibility of cables result in vibrations of the payload and hooks. In this paper, the forced vibrations due to rehabilitation exercises are studied. Since the previous studies of the authors showed a weak coupling between longitudinal and transverse vibrations, a two-cable planar model for the study of transverse vibrations is developed. The model makes it possible to study the forced transverse vibrations due to both cable motion and robot motion. Stiffness and damping of the patient’s arm are considered. Results show that the cable system exhibits a simple linear behavior when excited by robot motion and a non-linear behavior when excited by cable motion.

1. Introduction

Various movement-related disabilities are experienced each year by millions of individuals in the world [2] frequently as a result of sensory impairments, traumatic brain injuries, and musculoskeletal and neurological disorders [8, 20]. The benefits of rehabilitation and medical robotics are considered in many review articles, and meta-analyses [4, 7, 33, 37] demonstrating as robotic rehabilitation devices (RRDs) are well-suited to assist patients, based on their ability to carry out simple and repetitive tasks with consistency.

Cable-driven rehabilitation robots (CDRR) are commonly used in the rehabilitation field (e.g., Carex [21], CUBE [6], NereBot [14], CADEL [19, 35]), offering several promising features such as low inertia, high payload-to-weight ratio, and large workspace [5, 13, 16]. These advantages are largely due to the position of the actuators which are usually fixed to the ground, strongly reducing the mass to be moved: pulleys are used to make possible variations in the orientation of the cables, while the payload is connected to the cables via hooks.

Although CDRRs have a lot of promising features, there are also some limitations and deficiencies due to the intrinsic properties of cables, which result in unidirectional power transmission, vibrations, and maintenance. These limitations increase the complexity of kinematic and dynamic modeling of CDRRs [30]. The elasticity or flexibility of cables of CDRRs causes undesirable vibrations, which may generate position and orientation errors and compromise patient comfort. The importance of magnetic hooks is mostly related to operational flexibility. In fact, the orthosis (i.e., the payload of the rehabilitation robot) can be detached and can be fitted to a patient while another one is performing an exercise. Moreover, thanks to the hooks, different orthoses can be installed, both passive [28] and active [26].

Since there is a strong interaction between the patient and the robotic platform, it is fundamental to consider the characteristics of the patient’s dynamics [32]. The literature provides different values of the human arm dynamic properties, in particular for the stiffness $k_h$ and damping factor $c_h$ , depending on the field of application and therefore on the investigated population. The diversity in the fields of application leads to different postures assumed by the individuals which imply different values of the parameters.

Dolan et al. [9] focused on characterizing arm dynamics in order to properly design human-robot interaction systems. In particular, the authors provided experimental measurements of the damping and mass parameters, both in the unloaded and loaded cases. Dyck et Tavakoli [12] proposed a method to measure the dynamic impedance of the human arm to assess the motor function and muscle tone of impaired individuals. Since the sensors may increase the costs and complexity of these systems, the authors proposed a sensorless approach based on a virtual sensor derived from the robot’s kinematic and dynamic models. Woo and Lee [34] proposed a dynamic model of the human arm to describe the dynamic behavior of a haptic system interacting with the patient. The model, represented by a second-order mass-damper-spring system, shows that it is possible to exploit the patient’s damping to maintain passivity in the haptic system. To estimate the impedance parameters, the authors measured the force-displacement data when the patient was excited by small external disturbances during the interaction with the haptic device.

Human dynamic models are also frequent in vehicle dynamics, where the human body influences the dynamic behavior of the vehicle by means of both voluntary reactions and passive responses to vehicle oscillations. In particular, Doria et al. [11] presented different lumped-mass models to explain the measured behavior of the rider in experimental tests. The authors developed modular testing equipment composed of a hydraulic shaker and a motorcycle mockup to replicate the oscillation of a motorcycle, in particular roll oscillations. To identify the human mechanical parameters, the authors fitted the motion and torque frequency response functions with the developed models. Lastly, Höhne [15] expanded existing biomechanical pilot models to properly model the human dynamics and reduce the possibility of aircraft-pilot coupling problems. In particular, the author adopted data from the literature and a maximum likelihood cost function to identify the model parameters.

To the best of the author’s knowledge, no previous work on CDRRs has discussed the forced vibrations of cable robots equipped with hooks. Hence, this paper aims at filling this research gap. Another motivation behind this paper is to carry out the analysis of forced vibrations of CDRR due to rehabilitation exercises, considering also the stiffness and damping of the patient’s arm. Since massless cables are assumed, cable sagging is negligible. This assumption is acceptable for applications that do not involve large workspaces [17, 25].

Starting from the cable system of the Maribot, which is a typical cable-suspended rehabilitation robot, in Section 2, possible simplifications of the dynamic model are analyzed. The forced vibrations due to robot arm motion during rehabilitation exercises are studied in Section 3. In Section 4, the model is extended to simulate forced vibrations due to cable motion. In both cases, a harmonic input is considered, but in Section 5, the analysis is extended considering a polynomial input. Finally, in Section 6, results are discussed, and future applications of the cable models are illustrated.

2. Models for the analysis of vibrations of cable robots with hooks

The Maribot [27, 29] (Fig. 1a) is a 5-degrees-of-freedom (DOF) CDRR composed of a rigid planar jointed arm robot (2-DOF) and of a payload driven by cables (3-DOF). The three cables are controlled by DC motors fixed to the robot links, and they support an orthosis used to hold up the patient’s arm during rehabilitation exercises (Fig. 1b). Cable length and orientation are determined by the presence of pulleys. To improve the ease of use, each cable is connected to a magnetic hook near the orthosis, which makes the instantaneous release of the cable possible when a safety-limit force is overcome. Magnetic hooks are used in other CDRR like [3, 21, 23, 28].

Figure 1. Maribot rehabilitation robot: (a) Maribot; (b) orthosis and cable system of the robot.

To analyze the vibrations of this CDRR, first a single cable analysis was performed [38]. The mathematical model considered the elastic deformation of the cables and the inertia of the pulleys and had 4 DOFs. Figure 2 depicts the model and its DOFs: pulley rotation $\theta _1$ , rotation $\alpha _3$ of the upper cable, vertical translation of the orthosis $y_G$ , and vertical translation of the hook $y_3$ .

Figure 2. Single cable model: (a) longitudinal vibration; (b) transverse vibration.

Using the parameters of the single cable reported in Table I [38], the analysis of this model showed that:

• Transverse vibrations of the hook (caused by rotation $\alpha _3$ ) are completely uncoupled with longitudinal vibrations.

• Longitudinal vibrations are very little affected by pulley moment of inertia.

• The first longitudinal mode of vibration is very little affected by the hook mass.

• The natural frequencies of the longitudinal modes of vibration are higher than the natural frequency of the transverse modes of vibration.

Table I. Parameters of free vibration mathematical model.

Some experimental tests were carried out to find the natural frequencies of the system [38]. The measured results confirmed the calculated results.

Then, accordingly to the flowchart in Fig. 3, the model was extended considering two and three cables [36], but, in agreement with the results obtained with the single model, the moments of inertia of the pulleys were neglected.

Figure 3. Flowchart of the evolution of the mathematical model of the cable rehabilitation robot.

The planar model, which is depicted in Fig. 4, has 5 DOF: rotation $\theta _1$ of the left link of the four-bar linkage composed by the base, the two left cables ( $L_3+L_5$ ), the orthosis and the two right cables ( $L_4+L_6)$ ; rotations $\alpha _3$ and $\alpha _4$ of the cables that connect the base with the hooks; elongation of the left cable $s_9$ ; elongation of the right cable $s_{10}$ . Simulations showed the presence of two longitudinal modes of vibration and of three transverse modes of vibration. The most important results were:

• The transverse modes include the hook modes, in which the hooks move sideways whereas the orthosis essentially stands still, and a four-bar mode in which the left cable, the orthosis, and the right cable move as a four-bar linkage.

• The two hook modes have essentially the same natural frequency, which is much higher than the one of the four-bar mode.

• There is a weak or null coupling between the longitudinal and the transverse modes of vibration, the degree of coupling depends on the initial configuration of the system (if the system is a parallelogram the coupling is null).

• The natural frequencies of the longitudinal modes are higher than the natural frequencies of the transverse mode.

The last results actually confirmed and extended the results obtained by means of the single cable model. Then, some assumptions were made and the planar model of [36] was extended to study the free vibrations of a 3D symmetric system.

Figure 4. Planar model of the four-bar linkage with elongations and rotations of links.

The aim of the present research is the study of the forced vibrations of the suspended system caused by rehabilitation exercises. Rehabilitation exercises are performed cyclically varying both the vertical and horizontal position of the orthosis. The vertical motion is obtained by changing the length of the cables, whereas the horizontal motion is performed by the joints of the planar robot. Both motions can excite the vibrations of the suspended system, that is, both the cables and the orthosis. In typical rehabilitation exercises, the motions of the cables and of the robot are periodic with very low fundamental frequency (0.01–0.1 Hz) [28]. According to the results obtained with the single cable model and the planar model, these low-frequency motions cannot excite the longitudinal modes of the suspended system but can excite the transverse modes of vibrations.

For this reason, longitudinal compliance of the cables was neglected, and from the 5 DOF planar model of Fig. 4, a 3-DOF planar model of Figs. 5 and 9 was developed. Moreover, in the new model, the effect of patient arm stiffness and damping was added. The model was derived with the Langrange’s approach, and small modifications in the formulas of velocity make it possible to take into account both excitation due to robot motion (Section 3) and excitation due to cable motion (Section 4).

Figure 5. Model of transverse forced vibration due to robot arm motion. The input motion is represented in red and the DOFs $\theta _3$ , $\alpha _3$ , and $\alpha _4$ in blue.

The rehabilitation robot holds the patient’s arm during the exercise. However, the arm must not be simply considered as a lumped mass: in fact, the arm is attached to the patient’s body, which is then supported by the ground. As a result, the mass $M$ of the orthosis (with the arm) is connected to the ground also via a linear spring with stiffness $k_h$ and a damper with damping coefficient $c_h$ .

To define the values of $k_h$ and $c_h$ , many tests and models have been proposed during the years (Table II). In [11], the stiffness and damping of the arm are calculated considering the series between the upper arm and the forearm stiffness of a single upper limb, and the series between the upper arm and the forearm damping of a single upper limb. In [15], the stiffness of the arm is calculated considering the series between the translatory and the rotatory stiffness of the arm (the last one defined with a mean arm length of 0.7 m). A similar approach was adopted to evaluate arm damping. Ref. [34] presents sets of values of arm stiffness and damping of three subjects for different positions of the hand in the workspace and intensities of interaction force. Among those sets, the one with an interaction force set to zero and the location that most resembles the one depicted in Fig. 1 was considered. Therefore, the stiffness and damping of the patient’s arm are set equal to the mean value between the stiffness and damping of the arm of the three subjects, respectively.

Table II. Stiffness and damping of the patient’s arm.

The measured values vary widely, due to the different anatomical characteristics of the tested subjects and postures [18, 22], with an average value of $k_h$ = 78.66 N/m and $c_{h}$ = 7.40 Ns/m. In the model presented in this work, the stiffness and damping along the $x$ -axis ( $k_x$ and $c_x$ ) and along the $y$ -axis ( $k_y$ and $c_y$ ) are needed.

It is worth noticing that the values of Table II are found by experiments in which the movements of the arm are mostly restricted to the arm plane (i.e., the plane containing the arm and forearm axes) when the arm is horizontal. Very little literature has investigated arm stiffness in directions outside of such a plane [1, 31]. Indeed, the Maribot (Fig. 1a) is designed in such a way that motions in 3D space can be achieved. In this sense, the values of $k_h$ and $c_h$ of Table II should be adopted only for motions performed in the arm plane. However, [31] highlights that the impedance of the arm in the 3D space is characterized by components in $x$ - $y$ - $z$ -axis which are on the same order of magnitude. As a result, the values of Table II are considered as an overall representation of the arm suspended by the orthosis.

Finally, the Maribot is a 3-cable robot, but the behavior of a 3-cable robot can be studied by considering solely two cables [36]. As a result, $k_h$ and $c_h$ used in the following simulations will be considered as $2/3$ of the values of Table II, since only $2/3$ of the effect of the arm on the orthosis is supported by two cables.

3. Transverse vibration due to robot arm motion

3.1. Model description

When the robot moves in the horizontal plane with motion law $x_b(t)$ , the suspended orthosis behaves like a pendulum. To describe the transverse vibrations, the planar model with two cables of Fig. 5 was developed. The planar model is composed of two upper cables ( $L_{30}$ and $L_{40}$ ) which connect the pulleys to the magnetic hooks (mass $m_3$ and $m_4$ ). The negligible influence of the pulleys makes it possible to neglect the effect of the horizontal cables $L_1$ and $L_2$ that connect the motors to the pulleys, since they do not influence transverse vibrations. A payload equivalent to 2/3 of the expected load on the orthosis (mass $M$ ) is supported by two lower cables with length $L_{50}$ and $L_{60}$ . Similarly, the values of $k_x, k_y$ and $c_x, c_y$ adopted in the simulations are 2/3 of the values in Table II.

In the present analysis, the system is considered symmetric with $L_{30}=L_{40}$ , $L_{50}=L_{60}$ and $m_3=m_4$ . In the reference configuration, the rotations of the payload and the cables are equal to $\varphi _2 = 0$ and $\varphi _3 = \varphi _4 = \frac{3}{2}\pi$ . The quadrilateral $OABC$ is a parallelogram. The vibrations of the system about the reference configuration are described by the three coordinates depicted in Fig. 5 ( $\theta _3$ , $\alpha _3$ and $\alpha _4$ ). $\alpha _3$ and $\alpha _4$ are the rotations of the upper cables due to the lateral motion of the hooks. Since small oscillations are considered, these simple relations between the rotations of the upper and lower cables hold:

(1) \begin{equation} L_{30}\alpha _3 = L_{50}\alpha _5 \end{equation}

(2) \begin{equation} L_{40}\alpha _4 = L_{60}\alpha _6 \end{equation}

$\theta _3$ is the rotation of the quadrilateral $OABC$ with respect to the reference configuration. Rotations $\theta _2$ and $\theta _4$ are always related to $\theta _3$ by the loop equations of quadrilateral $OABC$ .

(3) \begin{equation} \left \{{\begin{array}{*{20}c}{OA \cos \left ( \varphi _3 + \theta _3 \right ) + b\cos \left ( \varphi _2 + \theta _2 \right ) - CB \cos (\varphi _4 + \theta _4) - c = 0} \\[5pt]{OA \sin \left ( \varphi _3 + \theta _3 \right ) + b\sin \left ( \varphi _2 + \theta _2 \right ) - CB \sin (\varphi _4 + \theta _4) = 0} \\ \end{array}} \right. \end{equation}

where $b$ is the length of the payload and $c$ is the distance between cable pulleys.

If $\alpha _3$ and $\alpha _4$ are small, the variations in the lengths of sides $OA$ and $CB$ due to the lateral motion of the hooks are small. Length $OA$ is given by:

(4) \begin{equation} OA = L_{30}cos(\alpha _3)+L_{50}cos(\alpha _5) \end{equation}

The Taylor’s expansion shows that the lateral motion of the hook has a second-order effect on length $OA$ .

(5) \begin{align} OA = L_{30}\left ( 1-\frac{\alpha _3^2}{2} \right ) + L_{50} \left ( 1-\frac{\alpha _5^2}{2} \right ) \end{align}

(6) \begin{align} OA = (L_{30}+L_{50})\left ( 1- \frac{L_{30}}{2L_{50}}\alpha _3^2 \right ) \end{align}

A similar equation holds true for $CB$ :

(7) \begin{equation} CB = \left(L_{30} + L_{50}\right) \left( 1- \frac{L_{30}}{2L_{50}}\alpha _4^2 \right) \end{equation}

If the symmetry conditions are introduced in the loop equations, the following equations hold:

(8) \begin{equation} \theta _4=\theta _3 \end{equation}

(9) \begin{equation} \theta _2 = \frac{(L_{30} + L_{50})(\alpha _4^2 L_{30}^2 - \alpha _3^2 L_{30}^2)}{2bL_{30}L_{50}} \end{equation}

It is worth noticing that the payload rotation is a second-order term.

The equations of motion of the cable system are developed with the Lagrange’s approach.

The system kinetic energy ( $E_k$ ) is equal to:

(10) \begin{equation} E_k = \frac{1}{2}m_3\left ({\dot x_3 ^2 + \dot y_3 ^2 } \right ) + \frac{1}{2}m_4\left ({\dot x_4 ^2 + \dot y_4 ^2 } \right ) + \frac{1}{2}M\left ({\dot x_G ^2 + \dot y_G ^2 } \right ) \end{equation}

The velocities of masses are calculated using first-order approximations:

(11) \begin{equation} \left \{{\begin{array}{*{20}c}{\dot x_3 = L_{30} \dot \alpha _3 + L_{30} \dot \theta _3 + \dot x_b } \\[4pt]{\dot y_3 = 0} \\ \end{array}} \right. \end{equation}

(12) \begin{equation} \left \{{\begin{array}{*{20}c}{\dot x_4 = L_{30} \dot \alpha _4 + L_{30} \dot \theta _3 + \dot x_b } \\[4pt]{\dot y_4 = 0} \\ \end{array}} \right. \end{equation}

(13) \begin{equation} \left \{{\begin{array}{*{20}c}{\dot x_G = \left ({L_{30} + L_{50} } \right )\dot \theta _3 + \dot x_b } \\[4pt]{\dot y_G = 0} \\ \end{array}} \right. \end{equation}

The springs and the dampers that represent the stiffness and the damping properties of the human arm have the first end-point fixed to the body and the second end-point fixed to the center of the orthosis (point G). The elastic potential energy ( $E_{p,el}$ ) is given by:

(14) \begin{equation} E_{p,el} = \frac{1}{2}k_x \left ({x_G - \frac{b}{2}} \right )^2 + \frac{1}{2}k_y \left ( y_G + (L_{30} + L_{50}) \right )^2 \end{equation}

where $k_x$ and $k_y$ are the horizontal and vertical arm stiffness; the gravity potential energy ( $E_{p,g}$ ) is given by:

(15) \begin{equation} E_{p,g} = m_3gy_3 + m_4gy_4 + Mgy_G \end{equation}

where the coordinates of the masses are given by:

(16) \begin{equation} y_3 = L_{30} sin \left ({3\pi }/{2} + \theta _3 + \alpha _3 \right ) = - L_{30} cos(\theta _3 + \alpha _3) \end{equation}

(17) \begin{equation} y_4 = L_{30} sin \left ({3\pi }/{2} + \theta _3 + \alpha _4 \right ) = - L_{30} cos(\theta _3 + \alpha _4) \end{equation}

(18) \begin{equation} \begin{cases} x_G = OA cos \left ({3\pi }/{2} + \theta _3 \right ) + x_b + \dfrac{b}{2}cos \left ( \theta _2 \right )= (L_{30} + L_{50}) \left ( 1 - \dfrac{L_{30}}{2L_{50}} \alpha _{3}^{2} \right ) sin \theta _3 + x_b + \dfrac{b}{2} \\[15pt] y_G = OA sin \left ({3\pi }/{2} + \theta _3 \right ) + \dfrac{b}{2} sin\left ( \theta _2 \right ) = - (L_{30} + L_{50}) \left ( 1- \dfrac{L_{30}}{2L_{50}}\alpha _3^2 \right ) cos \theta _3 + \dfrac{\left ( L_{30} + L_{50} \right )\left (\alpha _4^2 - \alpha _3^2 \right ) L_{30}}{4L_{50}} \end{cases} \end{equation}

If second-order Taylor’s expansions of the trigonometric functions are adopted, the coordinates needed for the calculations of the potential energy are given by the following equations:

(19) \begin{equation} y_3 = -L_{30} \left ( 1 - \frac{(\theta _3 + \alpha _3)^2}{2} \right ) \end{equation}

(20) \begin{equation} y_4 = -L_{30} \left ( 1- \frac{(\theta _3 + \alpha _4)^2}{2} \right ) \end{equation}

(21) \begin{equation} \begin{cases} x_G = \left ( L_{30} + L_{50} \right )\theta _3 + x_b + \dfrac{b}{2}\\[10pt] y_G = - \left ( L_{30} + L_{50} \right ) \left ( 1 - \dfrac{\theta _3^2}{2} - \dfrac{L_{30}}{2L_{50}} \alpha _3^2 \right ) + \dfrac{\left (L_{30} + L_{50} \right ) \left ( \alpha _4^2 - \alpha _3^2 \right ) L_{30}}{4L_{50}} \end{cases} \end{equation}

where in Equation 21, the term proportional to $\alpha _3^2\theta _3^2$ has been neglected since it is a fourth-order term.

Using the Lagrange’s approach, the equations of forced damped vibrations in matrix form are as follows:

(22) \begin{equation} {\bf M}_{\bf s} \cdot \ddot{\mathbf{q}} +{\bf C}_{\bf s} \cdot \dot{\mathbf{q}} +{\bf K}_{\bf s} \cdot \mathbf{q} = \mathbf{T} \end{equation}

where ${\bf M}_{\bf s}$ is the mass matrix and ${\bf C}_{\bf s}$ is the damping matrix, whereas $\mathbf{T}$ is the forcing torque and $\mathbf{q}=\{\theta _3,\alpha _3,\alpha _4\}^T$ . Since the stiffness matrix ${\bf K}_{\bf s}$ includes both the elastic and the gravitational terms, two matrices $K_e$ and $K_g$ are introduced to represent the elastic and gravitational parts, respectively. The mass, damping and the two stiffness matrices and the forcing torque of the model have the following form:

(23) \begin{equation} \begin{array}{*{20}c}{{\bf M}_{\bf s} = \left [{\begin{array}{*{20}c}{L_{30} ^2 \left ({m_3 + m_4 } \right ) + M\left ({L_{30} + L_{50} } \right )^2 } & \quad {L_{30} ^2 m_3 } & \quad {L_{30} ^2 m_4 } \\[5pt]{L_{30} ^2 m_3 } & \quad {L_{30} ^2 m_3 } & \quad 0 \\[5pt]{L_{30} ^2 m_4 } & \quad 0 & \quad {L_{30} ^2 m_4 } \\[5pt] \end{array}} \right ]} & \quad {{\bf C}_{\bf s} = \left [{\begin{array}{*{20}c}{c_x\left ( L_{30} + L_{50} \right )^2} & \quad 0 & \quad 0 \\[5pt] 0 & \quad 0 & \quad 0 \\[5pt] 0 & \quad 0 & \quad 0 \\[5pt] \end{array}} \right ]} \\[5pt] \end{array} \end{equation}

(24) \begin{equation}{\bf K}_{\bf e} = \left [{\begin{array}{*{20}c}{k_x\left ({L_{30} + L_{50} } \right )^2} & \quad 0 & \quad 0 \\[5pt] 0 & \quad 0 & \quad 0 \\[5pt] 0 & \quad 0 & \quad 0 \\[5pt] \end{array}} \right ] \end{equation}

(25) \begin{equation}{\bf K}_{\bf g} = \left [{\begin{array}{*{20}c}{g\left [{L_{30} \left ({m_3 + m_4 } \right ) + M\left ({L_{30} + L_{50} } \right )} \right ] } & \quad {m_3 gL_{30} } & \quad {m_4 gL_{30} } \\[5pt]{m_3 gL_{30} } & \quad {\dfrac{{\left ({2L_{50} m_3 + M\left ({L_{30} + L_{50} } \right )} \right )L_{30} g}}{{2L_{50} }}} & \quad 0 \\[5pt]{m_4 gL_{30} } & \quad 0 & \quad {\dfrac{{\left ({2L_{50} m_4 + M\left ({L_{30} + L_{50} } \right )} \right )L_{30} g}}{{2L_{50} }}} \\[5pt] \end{array}} \right ] \end{equation}

(26) \begin{equation}{\bf T} = - \left \{{\begin{array}{*{20}c}{\left [{L_{30} \left ({m_3 + m_4 } \right ) + M\left ({L_{30} + L_{50} } \right )} \right ]\ddot x_b + c_x\left ({L_{30} + L_{50} } \right )\dot x_b + k_x\left ({L_{30} + L_{50} } \right )x_b } \\[5pt]{L_{30} m_3 \ddot x_b } \\[5pt]{L_{30} m_4 \ddot x_b } \\[5pt] \end{array}} \right \} \end{equation}

Substituting Equation (21) into Equation (14), the difference $(y_G + (L_{30} + L_{50}))$ leads to fourth-order terms in $k_y$ , which are neglected in $K_e$ . A proportional damping assumption is made [24] and the damping matrix is assumed proportional to $K_e$ . It is worth noticing that the matrices are constant and the equations of motion describing the forced vibrations due to robot arm motion are a system of linear coupled differential equations.

The natural frequencies and the modes of vibration of this model are calculated solving the eigenvalue problem. The first mode is a pendulum mode of the whole system, whereas the second and third modes are the hook modes dominated by the transverse vibrations of the hook; they are similar to the modes found in [36] where a 5 DOF model of the free vibrations of a cable robot was developed.

As reported in Table II, very different values of the human arm stiffness $k_h$ and damping factor $c_h$ are present in the literature. The influence of arm stiffness on the natural frequencies of the modes of vibration is depicted in Fig. 6, considering the parameters of the planar model reported in Table III. The natural frequency of the pendular mode increases from 0.72 Hz (when arm stiffness is 0) to 1.48 Hz when the arm stiffness reaches the maximum value of Table II. Conversely, the variations in the natural frequencies of the 2 $^{nd}$ and 3 $^{rd}$ modes are negligible.

Table III. Parameters of the mathematical model for the transverse forced vibration due to robot arm motion.

Figure 6. Influence of the human arm stiffness $k_x$ on the natural frequency of first transverse mode (pendulum mode) (a) and of the second and third modes (hook modes) (b).

3.2. Response to robot arm motion

The simulated input motion of the robot $x_b$ is a sinusoidal signal:

(27) \begin{equation} x_b(t) = x_0 sin(\omega t) \end{equation}

with amplitude $x_0$ set to 0.150 m and angular frequency $\omega$ set to 0.6283 rad/s, corresponding to a frequency of 0.1 Hz. The reason for choosing such input motion is twofold:

• It is a realistic input;

• Since the system is linear, harmonic analysis can be used to study the effect of more complex inputs (e.g., periodic inputs)

Figure 7 deals with the forced vibration due to the robot motion considering four different models of arm stiffness and damping. Figure 7a and b represents the amplitude of the forcing torque $T$ for each DOF. Figure 7c and d represents the time histories of $\theta _3$ and $\alpha _3$ , respectively. $\alpha _4$ is not represented since it behaves as $\alpha _3$ (since the only difference in the equations of motion disappears due to the assumption $m_3=m_4$ ).

Figure 7. Transverse forced vibration due to robot motion for four values of stiffness and damping of patient’s arm: (a) amplitude of forcing $T_{\theta _3}$ for $\theta _3$ ; (b) amplitude of forcing $T_{\alpha _3}$ for $\alpha _3$ ; (c) time history of $\theta _3$ ; (d) time history of $\alpha _3$ .

The first stiffness and damping values are the ones from Doria et al. [11] and from Dyck et al. [12], since they present maximum stiffness (Doria et al.) and minimum stiffness (Dyck et al.), respectively. The third set of values includes the average values of $k_h$ and $c_h$ ( $k_{avg}$ = 78.66 N/m and $c_{avg}$ = 7.40 Ns/m). In the fourth model, the stiffness and damping values of the arm and forearm are set to zero, representing the absence of any reactions of the human arm.

The forcing torques are harmonic functions. The forcing torque $T_{\theta _3}$ on $\theta _3$ largely depends on arm stiffness and damping. In fact, for $k_h$ and $c_h$ = 0 the value is very small (in Fig. 7a it is amplified by a factor of 20). Moreover, since in this case $T_{\theta _3}$ depends only on $\ddot{x}_b$ , it is in opposition with respect to the torques referring to the other cases. The amplitude of the forcing torque $T_{\alpha _3}$ acting on $\alpha _3$ is very small and constant, since Equation (26) states that it depends only on $L_3$ , $m_3$ and on the acceleration of the robot $\ddot{x}_b$ .

The time domain response does not highlight resonance phenomena. $\theta _3$ chiefly oscillates at the forcing frequency (0.1 Hz). When arm damping is small, small amplitude transient vibrations at the frequency of the pendular mode appear as well. The response of $\alpha _3$ depicted in Fig. 7d is characterized by the presence of vibrations at the forcing frequency. There are small high-frequency vibrations at the natural frequency of the hook modes. When stiffness is large and damping is small, there are some vibrations at the frequency of the pendular mode as well. The amplitude of the main harmonic at the forcing frequency decreases if arm stiffness decreases and becomes very small when $k_h$ and $c_h$ are null (the corresponding curve in Fig. 7d is multiplied by a factor of 20). It is worth noticing that this phenomenon happens because, even if torque $T_{\alpha _3}$ does not depend on the stiffness and damping of the human arm, the equations of motion (22) are coupled by the mass and stiffness matrices.

These results are confirmed when the responses are analyzed in the frequency domain (Fig. 8). The main peak of the spectrum appears at the forcing frequency, and the peaks at the natural frequencies of the pendular mode and of the hook modes are very small.

Figure 8. Transverse forced vibration due to robot motion: (a) FFT of $\theta _3$ for the four cases considered and (b) of $\alpha _3$ for the four cases considered.

4. Forced vibrations due to cable motion

4.1. Model description

The two-cable model can be used to study the transverse forced vibrations due to cable motions having amplitude $\delta$ . In this case, the source of excitation is the periodic variation in the length of the upper cables, whereas the robot is considered steady ( $x_b = 0$ ), as shown in Fig. 9. In typical rehabilitation exercises, the cables move with the same motion law to keep the orthosis horizontal [28]; therefore, in the model $L_4(t) = L_3(t)$ . Kinematic equations (1-9) and (19-21) dealing with the position of the masses and the rotations of the cables still hold true if $L_{30}$ is replaced by $L_3(t)$ and if $x_b$ is set to zero.

Figure 9. Model of transverse forced vibration due to cables’ elongation. The input is represented in red, and the DOFs $\theta _3$ , $\alpha _3$ , and $\alpha _4$ are represented in blue.

The equations of the elastic potential energy (14) and gravity potential energy (15) still hold true, but it is assumed that the human spring is unloaded when the orthosis is at the lowest height.

The equation of the kinetic energy (10) still holds true, but new expressions of velocity are needed to take into account the elongation velocity of the cables ( $\dot L_4(t) = \dot L_3(t)$ ). The hooks have both elongation velocity ( $\dot L$ ) and tangential velocity,

(28) \begin{equation} \begin{cases} \dot x_G = \dot L_3 \theta _3 + \left ({L_{3} + L_{50} } \right )\dot \theta _3 - \dfrac{1}{2} b \theta _2 \dot \theta _2 \\[9pt] \dot y_G = \theta _3 \left ( L_3 + L_{50} \right ) \dot \theta _3 - \dfrac{1}{2} \dot L_3 + \dfrac{1}{4} \left ( -2\dot L_3 + 2 \theta _3^2 \dot L_3 \right ) + \dfrac{2\alpha _4L_3(L_3+L_{50})\dot \alpha _4 + \alpha _4^2L_3\dot L_3 + \alpha _4^2\left ( L_3 + L_{50} \right ) \dot L_3}{4L_{50}} \\[17pt] \quad \quad \quad +\dfrac{2\alpha _3L_3(L_3+L_{50})\dot \alpha _3 + \alpha _3^2L_3\dot L_3 + \alpha _3^2\left ( L_3 + L_{50} \right ) \dot L_3}{4L_{50}} \\[5pt] \end{cases} \end{equation}

(29) \begin{equation} \left \{{\begin{array}{*{20}c}{\dot x_3 = L_{3} \left ( \dot \alpha _3 + \dot \theta _3 \right ) cos\left ( \alpha _3 + \theta _3 \right ) + \dot L_3 sin\left ( \alpha _3 + \theta _3 \right )} \\[5pt]{\dot y_3 = L_{3} \left ( \dot \alpha _3 + \dot \theta _3 \right ) sin\left ( \alpha _3 + \theta _3 \right ) - \dot L_3 cos\left ( \alpha _3 + \theta _3 \right )} \end{array}} \right. \end{equation}

(30) \begin{equation} \left \{{\begin{array}{*{20}c}{\dot x_4 = L_{3} \left ( \dot \alpha _4 + \dot \theta _3 \right ) cos\left ( \alpha _4 + \theta _3 \right ) - \dot L_3 sin\left ( \alpha _4 + \theta _3 \right )} \\[5pt]{\dot y_4 = L_{3} \left ( \dot \alpha _4 + \dot \theta _3 \right ) sin\left ( \alpha _4 + \theta _3 \right ) + \dot L_3 cos\left ( \alpha _4 + \theta _3 \right )} \end{array}} \right. \end{equation}

The equations of motions of forced damped vibrations in matrix form are derived using the Lagrange’s approach and neglecting third (and higher) order terms. They have the form of Equation (22). Again the stiffness matrix $K_S$ is split into two matrices $K_e$ and $K_g$ including the elastic and gravitational terms, respectively. The mass, damping, and stiffness matrices and the forcing torque of the model have the following forms:

(31) \begin{equation} \begin{array}{*{20}c}{{\bf M}_{\bf s} = \left [{\begin{array}{*{20}c}{L_{3} ^2 \left ({m_3 + m_4 } \right ) + M\left ({L_{3} + L_{50} } \right )^2 } & \quad {L_{3} ^2 m_3 } & \quad {L_{3} ^2 m_4 } \\[9pt]{L_{3} ^2 m_3 } & \quad {L_{3} ^2 m_3 } & \quad 0 \\[9pt]{L_{3} ^2 m_4 } & \quad 0 & \quad {L_{3} ^2 m_4 } \\[9pt] \end{array}} \right ]} \end{array} \end{equation}

(32) \begin{equation} \begin{array}{*{20}c}{{\bf C}_{\bf s} = \left [{\begin{array}{*{20}c}{c_x\left ( L_{30} + L_{50} \right )^2 - c_y\left (L_3 + L_{50} \right )\left ( L_3 - L_{30} - \delta \right )}& \quad\!\!\!\! 0 & \quad\!\!\!\! 0 \\[9pt] 0 & \quad\!\!\!\! \dfrac{-c_y\left (L_3+L_{50}\right ) \left ( L_3 - L_{50} - \delta \right )L_3}{{2L_{50} }} & \quad\!\!\!\! 0 \\[9pt] 0 & \quad\!\!\!\! 0 & \quad\!\!\!\! \dfrac{-c_y\left (L_3+L_{50}\right ) \left ( L_3 - L_{50} - \delta \right )L_3}{{2L_{50} }} \\[9pt] \end{array}} \right ]} \\[9pt] \end{array} \end{equation}

(33) \begin{equation}{\bf K}_{\bf g} = \left [{\begin{array}{*{20}c}{g\left [{L_{3} \left ({m_3 + m_4 } \right ) + M\left ({L_{3} + L_{50} } \right )} \right ]}& \quad\!\!\!\! {m_3 gL_{3} } & \quad\!\!\!\! {m_4 gL_{3} } \\[9pt]{m_3 gL_{3} } & \quad\!\!\!\! {\dfrac{{\left ({2L_{50} m_3 + M\left ({L_{3} + L_{50} } \right )} \right )L_{3} g}}{{2L_{50} }}} & \quad\!\!\!\! 0 \\[9pt]{m_4 gL_{3} } & \quad\!\!\!\! 0 & \quad\!\!\!\! {\dfrac{{\left ({2L_{50} m_4 + M\left ({L_{3} + L_{50} } \right )} \right )L_{3} g}}{{2L_{50} }}} \\[9pt] \end{array}} \right ] \end{equation}

(34) \begin{equation}{\bf K}_{\bf e} = \left [{\begin{array}{*{20}c}{k_x\left ({L_{3} + L_{50} } \right )^2 -k_y\left (L_3 + L_{50} \right )\left ( L_3 - L_{30} - \delta \right )} & \quad\!\!\!\! 0 & \quad\!\!\!\! 0 \\[9pt] 0 & \quad\!\!\!\! {\dfrac{-k_y\left (L_3+L_{50}\right ) \left ( L_3 - L_{30} - \delta \right )L_3}{{2L_{50} }}} & \quad\!\!\!\! 0 \\[9pt] 0 & \quad\!\!\!\! 0 & \quad\!\!\!\! {\dfrac{-k_y\left (L_3+L_{50}\right ) \left ( L_3 - L_{30} - \delta \right )L_3}{{2L_{50} }}} \\[9pt] \end{array}} \right ] \end{equation}

(35) \begin{equation}{\bf T} = - \left \{{\begin{array}{*{20}c}{2{L_{3} \dot L_3 \left ( m_3 \left ( \dot \alpha _3 + \dot \theta _3 \right )+ m_4 \left ( \dot \alpha _4 + \dot \theta _3 \right )\right ) + 2M\left ( L_{3} + L_{50} \right ) \dot L_3 \dot \theta _3 }} \\[12pt]{2m_3L_3\dot L_3\left ( \dot \alpha _3 + \dot \theta _3 \right ) - \dfrac{1}{2}\dfrac{M\left (L_3 + L_{50}\right )\alpha _3 L_3 \ddot L_3}{L_{50}}}\\[12pt]{2m_4L_3\dot L_3 \left ( \dot \alpha _4 + \dot \theta _3 \right ) - \dfrac{1}{2}\dfrac{M\left (L_3 + L_{50}\right )\alpha _4L_3\ddot L_3}{L_{50}}} \\[5pt] \end{array}} \right \} \end{equation}

Due to the variation in cable lengths, the matrices of the systems depend on the configuration. It is worth noticing that the forcing torque on $\theta _3$ includes non-linear Coriolis terms, whereas the forcing torques on $\alpha _3$ and $\alpha _4$ include both Coriolis terms and terms related to cable acceleration. All the forcing terms depend either on $\dot \theta _3$ , $\dot \alpha _3$ , $\dot \alpha _4$ or on $\alpha _3$ , $\alpha _4$ . Hence, a non-null initial condition is needed to initiate the forced vibrations.

The natural frequencies and the modes of vibrations of this model are calculated by solving the eigenvalue problem for assigned configurations (cable lengths). Figure 10 shows the effect of cable length ( $L_4 = L_3$ ) on the natural frequencies of the transverse modes of vibration. It is interesting to note that the natural frequencies of the second and third transverse mode $f_{n2}$ and $f_{n3}$ assume the same value regardless of $k_h$ for a certain cable length, which corresponds to the lowest position of the orthosis (where the cable length is maximum). Indeed, in that point $L_3 = L_{30} + \delta$ , meaning that the terms in $k_y$ in the elastic part $K_e$ of the stiffness matrix are null.

Figure 10. Influence of the human arm stiffness $k_h$ on the natural frequency of first (a), second (b), and third (c) transverse mode.

On the other hand, the natural frequencies of the first transverse mode $f_{n1}$ for different values of $k_h$ do not reach a common value, but change during time keeping different values. This is reasonable as the first term of the elastic part $K_e$ of the stiffness matrix presents a component ( $k_x\left (L_3 + L_{50}\right )^2$ ) which does not become zero during the cycle.

Figure 11 represents influence of $L_3$ on the natural frequencies of the transverse modes. All the natural frequencies decrease as cable length increases. The natural frequencies of the second and third transverse modes (hook modes) assume the same value for the maximum elongation regardless of the human arm stiffness $k_h$ . On the other hand, the curves representing the value of the natural frequency of the first transverse mode for different values of $k_h$ do not intersect. It is worth noticing that the values of the natural frequencies for null values of the human arm stiffness and for $L_3 = 0.380$ m are equal to the ones obtained in [36].

Figure 11. Influence of the cable length $L_3$ on the natural frequency of first transverse mode (a), of the second one (b), and of the third one (c) for different human arm stiffness $k_h$ .

4.2. Response to cable motion

The simulated input motion of the cables is sinusoidal:

(36) \begin{equation} L_3(t) = L_{30}+\delta \sin (\omega _Lt) \end{equation}

amplitude $\delta$ is set to 0.150 m and angular frequency $\omega _L$ to 0.6283 rad/s, corresponding to 0.1 Hz. Initial cable length $L_{30}$ is set to 0.345 m.

Figure 12 deals with the forced vibration due to cable motion considering different models of arm stiffness and damping. Differently from Section 3.2, only two stiffness and damping models are considered, that is, the model provided by Doria et al. [10], which features large stiffness and minimum damping, and the model with null stiffness and damping. The results obtained considering the other models are deemed uninteresting since they show only an initial transient, due to the large damping values.

Figure 12. Transverse forced vibration due to cable motion considering two values of stiffness and damping of patient’s arm: (a) amplitude of forcing $T_{\theta _3}$ for $\theta _3$ ; (b) amplitude of forcing $T_{\alpha _3}$ for $\alpha _3$ ; (c) FFT of $F_{\theta _3}$ ; (d) FFT of $F_{\alpha _3}$ .

Moreover, only the forcing torques and the responses of $\theta _3$ and $\alpha _3$ are discussed, since $\alpha _4$ has almost the same behavior of $\alpha _3$ due to the symmetric configuration of the system. Figure 12a and b depict the forcing torques pf $\theta _3$ and $\alpha _3$ in the time domain, whereas Fig. 12c and d depict the same torques in the frequency domain. Differently from the previous scenario (Fig. 7), the forcing torques have complex waveforms with the presence of a beating phenomenon. This result can be explained by looking at the typical forcing term due to Coriolis acceleration, for example, the term

(37) \begin{equation} T_{\theta _3}^* = 2M(L_3(t) + L_{50})\dot L_3 \dot \theta _3 \end{equation}

in Equation (35).

If Equation (36) is inserted in (37), the following result holds:

(38) \begin{equation} T_{\theta _3}^* = 2M(L_3 + L_{50}) \delta \omega _L cos(\omega _L t) \dot \theta _3 + 2M \delta ^2 sin(\omega _L t)\omega _L cos(\omega _Lt)\dot \theta _3 \end{equation}

Since harmonic oscillations of $\theta _3$ at the natural frequency ( $\omega _n = 2\pi f_{n1}$ ) are expected

(39) \begin{equation} \theta _3 = \theta _{30}sin\omega _{n1}t \end{equation}

Equation (38) becomes:

(40) \begin{equation} T_{\theta _3}^* = 2M (L_{30}+L_{50})\delta \theta _{30}\omega _L\omega _{n1}cos(\omega _Lt)cos(\omega _{n1}t)+M\delta ^2\theta _{30}\omega _L\omega _{n1}sin(2\omega _Lt)cos(\omega _{n1}t) \end{equation}

where the second term has been modified using the double-angle formula. Then, the Werner’s formulae are used to transform the products of trigonometric functions into sums:

(41) \begin{equation} \begin{split} T_{\theta _3}^* &= 2M \left [ L_{30} + L_{50} \right ]\delta \theta _{30}\omega _L\omega _{n1} \frac{1}{2}\left [ cos((\omega _{n1}+\omega _L)t) + cos((\omega _{n1} - \omega _L)t) \right ] + \\ & + M\delta ^2\theta _{30}\omega _L\omega _{n1}\frac{1}{2}\left [ sin((\omega _{n1} + 2\omega _L)t) - sin((\omega _{n1} + \omega _L)t)\right ] \end{split} \end{equation}

Since in the present case $\omega _L \lt \lt \omega _{n1}$ , Equation (40) states that the terms $cos(\omega _Lt)$ and $sin(2\omega _Lt)$ modulate the input torque generating the beating phenomenon. Moreover, Equation (41) states that the modulated input torque consists of two harmonics, the former having frequency $\omega _{n1}-\omega _L$ , the latter having frequency $\omega _{n1}+\omega _L$ (if the modulation frequency is $2\omega _L$ the two frequencies are $\omega _{n1}-2\omega _L$ and $\omega _{n1}+2\omega _L$ ). This result is in agreement with the spectra depicted in Fig. 12c and d that show the presence of clusters of harmonics around the natural frequencies of the system. It is worth noticing that the behavior of the actual system is made more complex by some factors:

• the presence of damping;

• the coupling between the various DOFs;

• the variation of natural frequencies with cable length, which is a non-linear effect.

For the above-mentioned reasons, there are clusters of harmonics instead of two harmonics around the natural frequencies.

Figure 12 shows that the main effect of arm stiffness is the shift of the cluster of harmonics corresponding to $f_{n1}$ toward higher frequencies. Figure 12 highlights that the introduction of arm damping has a large effect on the torques on $\theta _3$ and $\alpha _3$ causing in the frequency domain a large reduction in the torque components having frequencies around $f_{n2}$ and $f_{n3}$ .

Figure 13. Transverse forced vibration due to cable motion for two values of stiffness and damping of patient’s arm: (a) time history of $\theta _3$ ; (b) time history of $\alpha _3$ ; (c) FFT of $\theta _3$ ; (d) FFT of $\alpha _3$ .

Figure 14. Comparison between the polynomial motion law and the sinusoidal one adopted in the previous sections.

Figure 15. Comparison of the responses due to robot arm motion for different input motions.

The responses of the system in the time and in the frequency domains are depicted in Fig. 13. The amplitude of the pendulum motion of the whole system (rotation $\theta _3$ ) is much smaller than the rotation of the upper cable caused by the lateral motion of the hook (rotation $\alpha _3$ ). If there is no arm damping, the vibrations are periodic and are characterized by beating phenomenon caused by the particular waveforms of the forcing torques. In the presence of arm damping, the oscillations quickly extinguish. In particular, the spectra depicted in Fig. 13c and d show that the high-frequency vibrations of the upper cables are more affected by arm damping than the low-frequency pendulum vibrations.

5. Forced vibrations caused by a polynomial motion law

From a practical point of view, it is interesting to analyze the response of the system to a generic input typical of an actual rehabilitation scenario [27, 28]. The input consists of a polynomial motion law having the same amplitude and fundamental frequency of the sinusoidal input. It is calculated interpolating with cubic splines $n$ via points defined along a cycle of oscillation. This kind of interpolation yields very regular trajectories with null speeds at each via point. In Fig. 14, the polynomial motion law and the sinusoidal one are compared. Figure 15 deals with the excitation caused by robot motion which is described by a linear model. The polynomial motion law has only the beneficial effect of reducing the vibrations at the natural frequencies of the transverse modes.

The excitation by cable motion, which is described by a non-linear model, could be more affected by the shape of the motion law. Figure 16 shows the forcing torques $T_{\theta _3}$ and $T_{\alpha _3}$ and the responses of $\theta _3$ and $\alpha _3$ in the time domain. The overall behavior of the system is similar to the one obtained with sinusoidal input in terms of vibration amplitudes and frequencies. But the different motion law causes modifications in the waveform of the torques that, in turn, cause modifications in the waveforms of the responses.

Figure 16. Generic cable motion. Transverse forced vibration due to cable motion for two values of stiffness and damping of patient’s arm: (a) time history of $T_{\theta _3}$ ; (b) time history of $T_{\alpha _3}$ ; (c) time history of $\theta _3$ ; (d) time history of $\alpha _3$ .

6. Discussion and conclusions

A planar two-cable model for the study of transverse vibrations of the CDRR has been developed. The matrices of the model can be quickly modified to simulate both the forced vibrations due to the lateral motion of the arm and the forced vibrations due to the variations in cable length (cable motion). First harmonic inputs were considered since they are realistic and useful to highlight the main features of the system’s response.

The system excited by arm motion has a linear behavior and shows the presence of lasting vibrations at the forcing frequency. The most important vibration is the pendulum motion of the whole cable system; there are also small transverse vibrations of the hooks.

The system excited by cable motion exhibits a non-linear behavior since the natural frequencies depend on cable lengths and the forcing torques include terms depending on cable rotations and angular velocities (Coriolis terms). With this excitation, hook modes are more excited than the four-bar mode. If realistic values of arm damping are considered, transverse vibrations quickly extinguish and have little effect on the rehabilitation exercise. The effect of the waveform of the input motion was considered as well. Results show that, if the amplitude of the input motion does not change, a sinusoidal and a polynomial motion law generate similar responses in terms of amplitudes and frequencies.

Future development will be the experimental validation of the forced model, assuming different motion laws and durations of the rehabilitation exercises. Moreover, a 3D model of the forced vibrations of the cable system will be developed, exploiting the symmetry properties of the system.