2. Mathematical model of the under-actuated robot

The robots studied in this paper are designed according to the theory of the differentially flat system [Reference Franch, Agrawal and Sangwan23]. A robot must satisfy certain criteria to be differentially flat. Specifically, in the case of an $n$ -DOF robot, the following conditions must be met:

• the center of mass of the final link ( $n$ ) must be positioned along the $n$ -th joint axis;

• the center of mass of both links $n$ and $n-1$ must be located along the $(n-1)$ -th joint axis.

These conditions are iteratively applied until the center of mass of the final $j$ links (i.e., $n, n-1, \ldots, n-j+1$ ) is positioned on the $(n-j+1)$ -th joint axis, as described in ref. [Reference Franch, Agrawal and Sangwan23]. As a result, the last $j$ links are called “fully balanced.” The actuators for the system are located on the $(n-j+1)$ -th joint, while the $(j-1)$ preceding joints are equipped with torsional springs. In this paper, only one non-actuated joint is considered, hence $j=2$ . Therefore, it is necessary to balance the robot up to the $n-j+1 = n-1$ joint. Fig. 1 shows the scheme of an $n$ -DOF robot with the last passive joint. $q_1$ represents the rotation of link 1 with respect to the base frame $x-y$ , and $q_i$ is the relative rotation of link $i$ with respect to $i-1$ link, with $i = 2,3,\dots, n$ . Only the last joint is not actuated, and the last link $n$ is connected with link $n-1$ by a torsional spring with stiffness $k_n$ , and by a viscous damper with damping coefficient $c_n$ .

Figure 1. Scheme of the mechanical system with $n-$ DOF.

Coming to mass distribution, the term $m_k$ is the mass of the $k$ -th link, $I_{Gk}$ is the barycentric moment of inertia of the $k$ -th element, $a_{Gk}$ is the distance of the center of mass of $k$ -th link from the $k$ -th joint and $a_k$ is the total length of the $k$ -th links (with $k = 1,2,\dots, n$ ). Regarding the last $n-1$ and $n$ balanced links, $m_{c(n-1)}$ and $m_{c(n)}$ are the counterbalancing masses of $n-1$ and $n$ links, respectively. Similarly, $a_{C(n-1)}$ and $a_{C(n)}$ are the distances of the counterbalancing masses from the $n-1$ and $n$ joint axis, respectively.

The dynamic model of the system can be derived by means of the Lagrangian approach, the result is a system of equations that can be written in matrix form:

(1) \begin{equation} \boldsymbol{M}_n(\boldsymbol{q}) \begin{bmatrix} \ddot{q}_1 \\[5pt] \ddot{q}_2 \\[5pt] \vdots \\[5pt] \ddot{q}_{n-2}\\[5pt] \ddot{q}_{n-1}\\[5pt] \ddot{q}_{n}\\[5pt] \end{bmatrix} + \boldsymbol{C}_n \begin{bmatrix} \dot{q}_1 \\[5pt] \dot{q}_2 \\[5pt] \vdots \\[5pt] \dot{q}_{n-2}\\[5pt] \dot{q}_{n-1}\\[5pt] \dot{q}_{n}\\[5pt] \end{bmatrix} + \boldsymbol{K}_n \begin{bmatrix} q_1 \\[5pt] q_2\\[5pt] \vdots \\[5pt] q_{n-2}\\[5pt] q_{n-1}\\[5pt] q_n \end{bmatrix} + \begin{bmatrix} b_1(\boldsymbol{q},\dot{\boldsymbol{q}}) \\[5pt] b_2(\boldsymbol{q},\dot{\boldsymbol{q}}) \\[5pt] \vdots \\[5pt] b_{n-2}(\boldsymbol{q},\dot{\boldsymbol{q}}) \\[5pt] 0\\[5pt] 0 \end{bmatrix}+ \begin{bmatrix} g_1(\boldsymbol{q}) \\[5pt] g_2(\boldsymbol{q}) \\[5pt] \vdots \\[5pt] g_{n-2}(\boldsymbol{q}) \\[5pt] 0\\[5pt] 0 \end{bmatrix} = \begin{bmatrix} \tau _1 \\[5pt] \tau _2 \\[5pt] \vdots \\[5pt] \tau _{n-2}\\[5pt] \tau _{n-1}\\[5pt] 0 \end{bmatrix} \end{equation}

in which the vector $ \boldsymbol{b}(\boldsymbol{q},\dot{\boldsymbol{q}}) = [b_1(\boldsymbol{q},\dot{\boldsymbol{q}}),b_2(\boldsymbol{q},\dot{\boldsymbol{q}}) \dots b_{n-2}(\boldsymbol{q},\dot{\boldsymbol{q}}),0,0]^T$ contains the Coriolis-centrifugal terms, while $ \boldsymbol{g}(\boldsymbol{q}) = [g_1(\boldsymbol{q}),g_2(\boldsymbol{q}) \dots g_{n-2}(\boldsymbol{q}),0,0]^T$ is the vector of gravity terms. It is worth noting that the last two elements of the vectors are null since the last $[n-1,n]$ links are fully balanced. $\tau _i$ is the motor torque applied on $i-$ th motor (with $i=1,2,\dots, n-1$ ). All the actuated joints are considered to be infinitely stiff and frictionless. Mass matrix $\boldsymbol{M}_n$ , damping matrix $\boldsymbol{C}_n$ , and stiffness matrix $\boldsymbol{K}_n$ are defined as follows:

(2) \begin{equation} \boldsymbol{M}_n (\boldsymbol{q}) = \left [ \begin{array}{c@{\quad}c@{\quad}c} & I_{n-1}^* & I_{n}^* \\[5pt] & I_{n-1}^* & I_{n}^* \\[5pt] \boldsymbol{I}^*_{(n-2)\times (n-2)}(\boldsymbol{q}) & \vdots & \vdots \\[5pt] & I_{n-1}^* & I_{n}^* \\[5pt] I_{n-1}^* \;\;\; \cdots \;\;\; I_{n-1}^* & I_{n-1}^* & I_{n}^*\\[5pt] I_{n}^* \;\;\; \cdots \;\;\; I_{n}^* & I_{n}^* & I_{n}^*\\[5pt] \end{array} \right ], \boldsymbol{C}_n = \left [ \begin{array}{c@{\quad}c} & 0 \\[5pt] \textbf{0}_{(n-1)\times (n-1)} & \vdots \\[5pt] & 0 \\[5pt] 0 \;\;\; \cdots \;\;\; 0 & c_n \end{array} \right ], \boldsymbol{K}_n = \left [ \begin{array}{c@{\quad}c} & 0 \\[5pt] \textbf{0}_{(n-1)\times (n-1)} & \vdots \\[5pt] & 0 \\[5pt] 0 \;\;\; \cdots \;\;\; 0 & k_n \end{array} \right ] \end{equation}

In Eq. (2), $\textbf{0}_{(n-1)\times (n-1)}$ is the zero matrix with dimension $n-1$ . Since the last $[n-1,n]$ links are fully balanced, the mass matrix of the system includes a $(n-2)\times (n-2)$ configuration-dependent sub-matrix, two constant rows and two constant columns. The constant columns are the transpose of the rows.

From Eq. (1), differential flatness can be used to define a set of flat variables [Reference Agrawal and Sangwan1]. For a $n$ -DOF system with one non-actuated joint, the vector of flat variables $\boldsymbol{Y} = [y_1,\ldots,y_{n-1}]$ can be defined as:

(3) \begin{equation} \begin{gathered} y_1 = \sum _{i=1}^n q_i\\[5pt] y_{i} = q_{i-1}\quad \text{with} \quad i = 2,3,\dots,n-1 \end{gathered} \end{equation}

for example, $n=4$ , $\boldsymbol{Y} = [q_1+q_2+q_3+q_4,q_1,q_2]$ . It is worth noting that flat variable $y_1$ has the physical meaning of orientation of the last link with respect to the fixed reference frame.

Flat variables can be used to control point-to-point motions of under-actuated systems and depend on the constraints on the joint values (position and derivatives) and on the total motion time $t_f$ . Although $\boldsymbol{Y}$ can be defined by means of any function, in most cases, polynomials are used [Reference Biagiotti and Melchiorri38]. In order to achieve the full controllability of the system, it is necessary to express all torques as a function of flat variables $\boldsymbol{Y}$ .

If $c_n$ is equal to zero, starting from the last $[n-1,n]$ rows of Eq. (1), it is possible to derive the angular position $q_n$ of the non-actuated joint as a function of the first flat variable:

(4) \begin{equation} q_n = -\frac{I_n^*\left (\ddot{q}_1+\ddot{q}_2+\dots +\ddot{q}_{n-1}+\ddot{q}_n\right )}{k_n} = -\frac{I_n^*}{k_n}\ddot{y}_1 \end{equation}

The motor torque applied to the $n-1$ joint can be calculated from the $n-1$ row of Eq. (1) and using Eq. (4). It results:

(5) \begin{equation} \tau _{n-1} = I_{n-1}^* \ddot{y}_1 +\frac{I_n^*(I_{n-1}^*-I_n^*)}{k_n}y^{(4)}_1 \end{equation}

Since $I_n^*$ and $I_{n-1}^*$ are constant, $\tau _{n-1}$ is a linear function of the derivatives of $y_1$ . The other motor torques can be calculated starting from $\tau _{(n-1)}$ by means of the Newton-Euler formula [Reference Sciavicco and Siciliano39]:Footnote 1

(6) \begin{equation}{\tau }_i = \sum _{k=1}^{n-2} \left (I_{ik}^*(\boldsymbol{q})-I_{(i+1)k}^*(\boldsymbol{q})\right ) \ddot{y}_{k+1} +(b_i(\boldsymbol{q},\dot{\boldsymbol{q}})-b_{i+1}(\boldsymbol{q},\dot{\boldsymbol{q}})) +(g_i(\boldsymbol{q})-g_{i+1}(\boldsymbol{q}))+{\tau }_{i+1} \end{equation}

with $i=1,2,\dots n-2$ . The flat variable ${y}_{k+1}$ has the subscript $k+1$ because the first flat variable is defined as the sum of all angular positions of the robot. Owing to the presence of the Coriolis-centrifugal terms, the gravity terms, and the configuration-dependent inertial terms, these motor torques are non-linear functions of flat variables.

Only under the hypothesis of a frictionless non-actuated joint, Equation (4) can be obtained. If $c_n \neq 0$ , the analytical expression of $q_n$ as a function of only $y_1$ is difficult to obtain. In the next section, an analytical method in the Laplace domain to express the joint variable $q_n$ as a function of a flat variable in the presence of viscous friction is presented.

3. Laplace domain analysis of the last two motion equations

The last two rows of Eq. (1) are linear since the last $[n-1,n]$ links satisfy the differential flatness rules and are fully balanced. Hence, the Laplace transform can be rigorously calculated.

3.1. Motor torque calculation in the presence of viscous damping

The flat variable $y_1 = \sum _{i=1}^n q_i$ with $i=1,2,\dots,n$ , and its derivatives can be written in the Laplace domain as follows:

(7) \begin{equation} \begin{aligned} \mathcal{L}(y_1 (t))& = Y_1(s) = \sum _{i=1}^n Q_i(s) \\[5pt] \mathcal{L}(\dot{y}_1(t))& = s Y_1(s) =s \sum _{i=1}^n Q_i(s)\\[5pt] \mathcal{L}(\ddot{y}_1(t))& = s^2 Y_1(s) = s^2 \sum _{i=1}^n Q_i(s) \end{aligned} \end{equation}

in which $\mathcal{L}$ means Laplace transform.

The last two equations of motion considering the damping coefficient $c_n\neq 0$ can be rewritten as follows:

(8) \begin{equation} \begin{cases} I_{n-1}^*s^2(Y_1(s)-Q_n(s))+ I_n^* s^2Q_n(s) = T_{n-1}(s)\\[5pt] I_{n}^*s^2Y_1(s)+(c_n s+k_n)Q_n(s) = 0\\[5pt] \end{cases} \end{equation}

in which $T_{n-1}(s)$ is the motor torque of the last non-actuated joint. The joint variable $Q_n(s)$ can be obtained from the last row of Eq. (8):

(9) \begin{equation} Q_n(s) =-\frac{ I_{n}^*}{c_n s+k_n}s^2Y_1(s) \end{equation}

The torque that the last actuated joint has to generate can be obtained from the first row of Eq. (8) introducing Eq. (9):

(10) \begin{equation} \tilde{T}_{n-1}(s) = I_{n-1}^* s^2Y_1(s) +\frac{I_n^*(I_{n-1}^*-I_n^*) }{c_n s+k_n} s^4 Y_1(s) \end{equation}

in which $\tilde{T}_{(n-1)}(s)$ is the motor torque that takes into account the damping coefficient and $\tilde{k}_n = k_n+sc_n$ is the complex stiffness of the non-actuated joint.

Equation (10) can be rewritten collecting the stiffness $k_n$ from the second term, it results:

(11) \begin{equation} \tilde{T}_{n-1}(s) = I_{n-1}^* s^2Y_1(s) +\frac{I_n^*(I_{n-1}^*-I_n^*) }{k_n} \left (1+\frac{c_n}{k_n}s\right )^{-1} s^4 Y_1(s) \end{equation}

Usually, the damping coefficient $c_n$ is much smaller than torsion spring stiffness, especially if low-friction bearings are used. For this reason, the term $(1+c_n/k_n\,s)^{-1}$ can be approximated using its Taylor expansion arrested at the first order. Hence, Eq. (11) can be written as follows:

(12) \begin{equation} \tilde{T}_{n-1}(s) = I_{n-1}^* s^2 Y_1(s) +\frac{I_n^*(I_{n-1}^*-I_n^*) }{k_n} \left (1-\frac{c_n}{k_n}s\right ) s^4 Y_1(s) \end{equation}

Hence, the expression of the motor torque $\tilde{T}_{n-1}$ in the Laplace domain takes the final form:

(13) \begin{equation} \tilde{T}_{n-1} (s) = I_{n-1}^* s^2 Y_1(s) +\frac{I_n^*(I_{n-1}^*-I_n^*) }{k_n} s^4 Y_1(s) - \frac{I_n^*(I_{n-1}^*-I_n^*) c_n }{k_n^2} s^5 Y_1(s) \end{equation}

Finally, it is possible to express motor torque $\tilde{T}_{n-1}(s)$ in the time domain by means of the inverse Laplace transform (the initial condition are zero):

(14) \begin{equation} \tilde{\tau }_{n-1} (t) = I_{n-1}^* \ddot{y}_1 + \frac{ I_n^*(I_{n-1}^*-I_n^*) }{k_n} y_1^{(4)} - \frac{c_n I_n^*(I_{n-1}^*-I_n^*)}{k_n^2} y_1^{(5)} \end{equation}

It is worth noting that the introduction of viscous friction leads to an increase in the order of the derivative of the flat variable that appears in the expression of the motor torque. Hence, the degree of the polynomial representing $y_1$ must increase to satisfy constraints on the fifth derivative. Equation (6) is still valid if the damped motor torque $\tilde{\tau }_{n-1}$ is used instead of $\tau _{n-1}$ :

(15) \begin{equation} \tilde{\tau }_i(t) = \sum _{k=1}^{n-2} \left (I_{ik}^*-I_{(i+1)k}^*\right ) \ddot{y}_{k+1}+(b_i(\boldsymbol{q},\dot{\boldsymbol{q}})-b_{i+1}(\boldsymbol{q},\dot{\boldsymbol{q}}))+(g_i(\boldsymbol{q})-g_{i+1}(\boldsymbol{q}))+ \tilde{\tau }_{i+1} \end{equation}

with $i=1,2,\dots n-2$ .

3.2. Natural frequencies

Free vibrations analysis makes it possible to identify the resonances of the system and is preliminary to forced vibrations analysis.

In the actual operation of the robot, the sum of the joint variables is an output that depends on the motor torques which are functions of the commanded flat variables. Hence, in the analysis of natural frequencies and forced response, the sum of joint variables in the Laplace domain is named $\Theta _n(s)$ .

The last equations of motion considering the damping coefficient ( $c_n\neq 0$ ) can be rewritten as follows in the Laplace domain:

(16) \begin{equation} \begin{cases} s^2 I_{n-1}^*(\Theta _n(s)-Q_n(s))+s^2 I_n^* Q_n(s) = T_{n-1}(s)\\[5pt] s^2 I_{n}^*\Theta _n(s)+(c_n s+k_n)Q_n(s) = 0\\[5pt] \end{cases} \end{equation}

in which $T_{n-1}(s)$ is the motor torque for the last non-actuated joint.

If the motor torque $T_{n-1}(s)$ and $c_n$ are set to zero, the equations of free undamped vibrations are obtained. The natural frequencies of the undamped system can be determined by finding the roots of the determinant of the matrix that multiplies the vector of unknowns:

(17) \begin{equation} \begin{bmatrix} I_{n-1}^*s^2 & (I_{n}^*-I_{n-1}^*)s^2 \\[5pt] I_{n}^*s^2 & k_n \\[5pt] \end{bmatrix} \begin{bmatrix} \Theta _1(s)\\[5pt] Q_{n}(s) \end{bmatrix} = \begin{bmatrix} 0\\[5pt] 0 \end{bmatrix} \end{equation}

The determinant of the matrix at the left-hand side of Eq. (17) is:

(18) \begin{equation} s^2(I_{n-1}^* k_n - I_{n}^*(I_{n}^*-I_{n-1}^*)s^2) = 0 \end{equation}

Equation (18) shows that one solution is null since it corresponds to a rigid motion of the system. The other non-null solution is the natural frequency of the system:

(19) \begin{equation} \omega _n = \sqrt{\frac{k_n I_{n-1}^*}{I_{n}^*(I_{n-1}^*-I_{n}^*)}} \end{equation}

3.3. Forced vibrations

Coming back to Eq. (16), the transfer function of the system between the joint variable $Q_n(s)$ and the generic input $T_{n-1}(s)$ can be calculated as follows:

(20) \begin{equation} Q_n(s) = \frac{I_n^*}{ -I_{n}^*(I_{n}^*-I_{n-1}^*)s^2+(c_n s+k_n)I_{n-1}^*} T_{n-1}(s) \end{equation}

Since Eq. (20) shows that the transfer function $Q_n(s)/T_{n-1}(s)$ has a pair of complex conjugate poles, in general, the response is characterized by damped oscillations.

If the particular torque input, calculated with the flat variable (Eq. (10)), is considered:

(21) \begin{equation} \tilde{T}_{n-1}(s) = \frac{ I_n^*(I_{n}^*-I_{n-1}^*)s^2-(c_n s+k_n)I_{n-1}^*}{c_ns+k_n} s^2 Y_1(s) \end{equation}

when Eq. (21) is introduced in Eq. (20), a zero-pole cancelation takes place:

(22) \begin{equation} Q_n(s) = -\frac{I_n^*}{c_n s+k_n} s^2 Y_1(s) \end{equation}

which corresponds to the disappearance of the oscillations.

Conversely, if the torque input is calculated by setting the damping coefficient $c_n$ equal to $0$ , Eq. (10) takes the following form:

(23) \begin{equation} T_{n-1}(s) = \frac{ I_n^*(I_{n}^*-I_{n-1}^*)s^2- k_n I_{n-1}^*}{k_n} s^2 Y_1(s) \end{equation}

When this torque is introduced in Eq. (20), the zero-pole cancelation does not take place, and owing to the complex conjugate poles oscillations may occur:

(24) \begin{equation} Q_n(s) = \frac{I_n^*}{ -I_{n}^*(I_{n}^*-I_{n-1}^*)s^2+(c_n s+k_n)I_{n-1}^*} \cdot \frac{ I_n^*(I_{n}^*-I_{n-1}^*)s^2- k_n I_{n-1}^*}{k_n} s^2 Y_1(s) \end{equation}

4. Numerical and experimental results for a $2$ -DOF robot

The robot designed for the experimental validation of the mathematical model has $2$ -DOF ( $n=2$ ), with one motor and one non-actuated joint. It moves on the horizontal plane, and link 2 is fully balanced. The first link does not need to be balanced, since its moment of inertia about joint $1$ is always constant [Reference Agrawal and Sangwan1].

For such a robot, Eq. (1) takes the following form:

(25) \begin{equation} \boldsymbol{M}_2 \begin{bmatrix} \ddot{q}_1 \\[5pt] \ddot{q}_2 \end{bmatrix} + \boldsymbol{C}_2 \begin{bmatrix} \dot{q}_1 \\[5pt] \dot{q}_2 \end{bmatrix} + \boldsymbol{K}_2 \begin{bmatrix} q_1 \\[5pt] q_2 \end{bmatrix} = \begin{bmatrix} \tau _1 \\[5pt] 0 \end{bmatrix} \end{equation}

where $\tau _1$ is the motor torque and the matrices are defined as follows:

(26) \begin{equation} \boldsymbol{M}_2 = \begin{bmatrix} I_1^* & I_2^* \\[5pt] I_2^* & I_2^* \\[5pt] \end{bmatrix}, \boldsymbol{C}_2 = \begin{bmatrix} 0 & 0 \\[5pt] 0 & c_2 \\[5pt] \end{bmatrix}, \boldsymbol{K}_2 = \begin{bmatrix} 0 & 0 \\[5pt] 0 & k_2 \\[5pt] \end{bmatrix} \end{equation}

in which the moments of inertia are $I_1^* = I_{G1} + m_1 a_{G1}^2 + (m_2+m_{c2}) a_{1}^2 + I_{G2}+m_2a_{G2}^2+m_{c2}a_{c2}^2$ , and $I_2^* = I_{G2}+m_2a_{G2}^2+m_{c2}a_{c2}^2$ , and $c_2$ is the damping coefficient of the second joint. It is worth noting that, in this case, the mass matrix is constant and the Coriolis-centrifugal and gravity terms are null because the second link is balanced and the system works in the horizontal plane.

Only one flat variable is required: $y_1=q_1+q_2$ . If $c_2 = 0$ , the motor torque of the first joint is calculated by means of Eq. (5):

(27) \begin{equation} \tau _1 = I_1^* \ddot{y}_1 +\frac{I_2^*(I_1^*-I_2^*)}{k_2}y^{(4)}_1 \end{equation}

whereas the motor torque for the damped systems is calculated by means of Eq. (14):

(28) \begin{equation} \tilde{\tau }_{1} = I_{1}^* \ddot{y}_1 + \frac{I_2^*(I_{1}^*-I_2^*) }{k_2} y_1^{(4)} - \frac{c_2 I_2^*(I_{1}^*-I_2^*) }{k_2^2} y_1^{(5)} \end{equation}

The flat variable is chosen to obtain a specific position, velocity, acceleration, jerk, and snap of $q_1$ and $q_2$ at the beginning and at the end of the motion ( $t_f$ ). For the undamped system, a $9$ -th degree polynomial is sufficient, and the corresponding constraints on the flat variable $y_1$ are as follows:

(29) \begin{equation} \begin{cases} y_1(t_i) = q_{1i}+q_{2i}\\[5pt] y_1(t_f) = q_{1f}+q_{2f} \\[5pt] \dot{y}_1(t_i) = \ddot{y}_1(t_i) = y_1^{(3)}(t_i) = y_1^{(4)}(t_i) = 0\\[5pt] \dot{y}_1(t_f) = \ddot{y}_1(t_f) = y_1^{(3)}(t_f) = y_1^{(4)}(t_f) = 0 \end{cases} \end{equation}

where $q_{1i}$ , $q_{1f}$ , $q_{2i}$ , and $q_{2f}$ are the initial and final angular positions of joint 1 and joint 2, respectively. For the damped system, at least an $11$ -th degree polynomial must be used to obtain a continuous fifth derivative $y_1^{(5)}$ . Hence, $y_1$ is defined by an $11$ -th degree polynomial, and the boundary conditions are the same as Eq. (29) with the addition of the constraints on the crackle:

(30) \begin{equation} \begin{cases} y^{(5)}_1(t_i) = 0\\[5pt] y^{(5)}_1(t_f) = 0 \end{cases} \end{equation}

The constants of the polynomial are obtained by solving the linear system:

(31) \begin{equation} \boldsymbol{A}\boldsymbol{x} = \boldsymbol{b} \end{equation}

in which $\boldsymbol{A}$ is the matrix associated with the boundary conditions, $\boldsymbol{x}$ is the vector of the polynomial coefficients and $\boldsymbol{b}$ is the vector that contains the boundary conditions [Reference Jazar40].

The prototype was built using additive manufacturing and moves in the horizontal plane as depicted in Fig. 2. The first joint is actuated by means of a direct drive Portescap 35NT2R82 426SP brushed DC motor. The non-actuated joint is equipped with a small torsional spring with a stiffness of 0.0026 $Nm/rad$ identified via modal analysis directly on the prototype [Reference Tonan, Bottin, Doria and Rosati41]. The inertial properties of links 1 and 2 are summarized in Table I. It should be noted that the $m_i$ term represents the mass of the $i$ -th link without considering any balancing mass $m_{ci}$ . The details of stiffness, damping coefficient, and the commanded point-to-point motion are listed in Table II. It is worth noting that, for a point-to-point motion, the initial and final angles of the under-actuated joint (i.e., $q_{2i}$ and $q_{2f}$ ) must be null [Reference Bottin and Rosati26]. The damping coefficient of the under-actuated joint $c_2$ is identified by means of free vibration tests on the prototype using the logarithmic decrement method [Reference Kelly42]. The experimental natural frequencies of the prototype robot are $0$ and $1.4\,\,Hz$ . They are in good agreement with the analytical values calculated by means of the formulas of Section 3.2.

Figure 2. The $2$ -DOF robot used in the experimental tests. The motor of the device is located under the black plane to avoid interference with the vision system.

The response of the prototype robot to the torques calculated by means of the standard and the improved diffeomorphism (with viscous friction) is studied both experimentally and numerically. For the experimental validation, the measurement of joint variables is needed. The angle of the actuated joint ( $q_1$ ) is obtained by means of an incremental optical encoder LIKA Electronic I58-L-5000ZCU16R installed on the motor of joint 1. A second encoder, for the measurement of the passive joint angle ( $q_2$ ), would represent an added mass. Hence, the second link is equipped with two markers, and its trajectory is tracked by the industrial camera Dalsa Genie Nano GM30-M2050, having a resolution of $2064\times 1544$ pixels and using a sensor Sony IMX252.

For the simulations, the model of Eq. (25) with the parameters of Tables I and II is implemented in Simulink.

Three scenarios are studied:

• Control based on the simple diffeomorphism ( $c_2=0$ ) with $y_1$ represented by a $9$ -th degree polynomial.

• Control based on the simple diffeomorphism ( $c_2=0$ ) with $y_1$ represented by an $11$ -th degree polynomial.

• Control based on the new diffeomorphism ( $c_2\neq 0$ ) with $y_1$ represented by an $11$ -th degree polynomial.

Figs. 3 and 4 depict the simulated and measured motion when the joint 1 torque is calculated according to the simple diffeomorphism ( $c_2=0$ ) with a $9$ -th degree polynomial, and with an $11$ -th degree polynomial, respectively. In both cases, there is a good agreement between numerical and experimental results. The end-effector of the robot describes a trajectory with a cusp. It is worth noting that the area swept by the robot is smaller than the area swept by a $1$ -DOF robot made by a unique link with length $a_1+a_2$ . The increase in the order of the polynomial modifies the shape of the trajectory near the cusp. In both cases, the end-effector has an overshoot and before stopping oscillates around the final position. It is worth noting that the robot is not aligned to the $x$ -axis at the end of the planned motion ( $t=0.6\,\,s$ ), since it is still oscillating.

Table I. Inertial parameters of the simulated and experimental robot.

Table II. Stiffness and damping coefficient of the non-actuated joint used in the simulated and experimental robot and the motion parameters.

Figure 3. $2$ -DOF robot trajectory of the end-effector in $x-y$ plane obtained with $y_1$ described by a $9$ -th degree polynomial, without the viscous friction on the non-actuated joint ( $c_2=0$ ). The black trajectory is obtained by a $1$ -DOF robot with only one link of length $a_1+a_2$ .

Figure 4. $2$ -DOF robot trajectory of the end-effector in $x-y$ plane obtained with $y_1$ described by an $11$ -th degree polynomial, without the viscous friction on the non-actuated joint ( $c_2=0$ ). The black trajectory is obtained by a $1$ -DOF robot with only one link of length $a_1+a_2$ .

Fig. 5 shows the effect of the torque calculated taking into account joint viscous friction ( $c_2\neq 0$ ). In this case, the overshoot and the oscillations are strongly reduced and at the end of the planned motion, the two links are aligned. The agreement between experimental and numerical results is good.

Figure 5. $2$ -DOF robot trajectory of the end-effector in $x-y$ plane obtained with $y_1$ described by an $11$ -th degree polynomial, with the viscous friction on the non-actuated joint ( $c_2\neq 0$ ). The black trajectory is obtained by a $1$ -DOF robot with only one link of length $a_1+a_2$ .

It is worth noticing that the shape of the trajectory changes with the polynomial degree and with the inclusion of viscous friction. End-effector trajectory is also influenced by the non-actuated joint stiffness and the motion time, as shown in ref. [Reference Tonan, Doria, Bottin and Rosati9].

To highlight the effect of the diffeomorphism on the oscillation of the robot, the simulated and measured joint variables are depicted in Fig. 6. There is a good agreement between the two sets of results. The comparison between $q_1$ and $q_2$ depicted in Fig. 6 shows that the oscillations at the end of the planned motion chiefly involve the second joint which is non-actuated. The right side of Fig. 6 shows that if the simple diffeomorphism is adopted, there are large oscillations of $q_2$ whose frequency is equal to $\omega _n$ . The increase in the order of the polynomial has a small effect on the amplitude of such oscillations. Conversely, the adoption of the new diffeomorphism leads to a nearly complete disappearance of oscillations of $q_2$ . The suppression of oscillations at the end of the trajectory is due to the zero-pole cancelation described in Section 3.3.

Figure 6. $2$ -DOF robot numerical and experimental joint rotations (left: $q_1$ ; right: $q_2$ ) in the case of three different open-loop controls: frictionless control with $9$ -th degree polynomial law for $y_1$ (top), frictionless control with $11$ -th degree polynomial (middle), friction compensation control with $11$ -th degree polynomial (bottom).

5. Numerical results for a $4$ -DOF robot

This section aims to show the application of the general theory to a more complex under-actuated robot with $4$ -DOF having non-linear terms and operating in the horizontal plane. Trajectories are simulated to highlight the effect of torques calculated with and without viscous friction. For this robot, Eq. (1) becomes:

(32) \begin{equation} \boldsymbol{M}_4 (\boldsymbol{q}) \begin{bmatrix} \ddot{q}_1 \\[5pt] \ddot{q}_2 \\[5pt] \ddot{q}_3 \\[5pt] \ddot{q}_4 \end{bmatrix} + \boldsymbol{C}_4 \begin{bmatrix} \dot{q}_1 \\[5pt] \dot{q}_2 \\[5pt] \dot{q}_3 \\[5pt] \dot{q}_4 \end{bmatrix} + \boldsymbol{K}_4 \begin{bmatrix} q_1 \\[5pt] q_2 \\[5pt] q_3 \\[5pt] q_4 \end{bmatrix} + \begin{bmatrix} b_1(\boldsymbol{q},\dot{\boldsymbol{q}}) \\[5pt] b_2(\boldsymbol{q},\dot{\boldsymbol{q}}) \\[5pt] 0 \\[5pt] 0 \end{bmatrix}= \boldsymbol{\tau } \end{equation}

in which $\boldsymbol{M}_4 (\boldsymbol{q})$ , $\boldsymbol{C}_4$ , and $\boldsymbol{K}_4$ are the mass, damping, and stiffness matrices, respectively. These matrices are defined as follows:

(33) \begin{equation} \boldsymbol{M}_4 (\boldsymbol{q}) = \begin{bmatrix} I_{11}^*(\boldsymbol{q}) & I_{12}^*(\boldsymbol{q}) & I_{3}^* & I_{4}^* \\[5pt] I_{21}^*(\boldsymbol{q}) & I_{22}^*(\boldsymbol{q}) & I_{3}^* & I_{4}^* \\[5pt] I_{3}^* & I_{3}^* & I_{3}^* & I_{4}^* \\[5pt] I_{4}^* & I_{4}^* & I_{4}^* & I_{4}^* \end{bmatrix} \boldsymbol{C}_4 = \left [ \begin{array}{cc} & 0 \\[5pt] \textbf{0}_{3\times 3} & \vdots \\[5pt] & 0 \\[5pt] 0 \;\;\; \cdots \;\;\; 0 & c_4 \end{array} \right ], \quad \boldsymbol{K}_4 = \left [ \begin{array}{c@{\quad}c} & 0 \\[5pt] \textbf{0}_{3\times 3} & \vdots \\[5pt] & 0 \\[5pt] 0 \;\; \cdots \;\; 0 & k_4 \end{array} \right ] \end{equation}

The elements of the mass matrix $\boldsymbol{M}_4 (\boldsymbol{q})$ and of the vector of Coriolis-centrifugal terms are rather complex and are reported in the Appendix.

The vector $\boldsymbol{\tau }$ contains the motor torques of the first, second, and third actuated joints. According to Eq. (3), the three flat variables are $\boldsymbol{Y}=[y_1,y_2,y_3] = [q_1+q_2+q_3+q_4, q_1, q_2]$ . The motor torques without considering the viscous friction can be calculated from Equations (5) and (6) imposing $n=4$ .

(34) \begin{equation} \begin{gathered}{\tau }_3 = I_{3}^* \ddot{y}_1 + \frac{ I_4^*(I_{3}^*-I_4^*) }{k_4} y_1^{(4)}\\[5pt] {\tau }_2 = \left (I_{21}^*(y_3)-I_{3}^*\right ) \ddot{y}_{2}+\left (I_{22}^*(y_3)-I_{3}^*\right ) \ddot{y}_{3}+b_2(y_3,\dot{y}_2)+{\tau }_{3} \\[5pt] {\tau }_1 = \left (I_{11}^*(y_3)-I_{21}^*(y_3)\right ) \ddot{y}_{2}+\left (I_{12}^*(y_3)-I_{22}^*(y_3)\right ) \ddot{y}_{3}+\left (b_1(y_3,\dot{y}_2,\dot{y}_3)-b_2(y_3,\dot{y}_2)\right )+{\tau }_{2}\\[5pt] \end{gathered} \end{equation}

The motor torques with viscous friction can be calculated from Equations (14) and (15) imposing $n=4$ , the expressions are as follows:

(35) \begin{equation} \begin{gathered} \tilde{\tau }_3 = I_{3}^* \ddot{y}_1 + \frac{ I_4^*(I_{3}^*-I_4^*) }{k_4} y_1^{(4)} - \frac{c_4 I_4^*(I_{3}^*-I_4^*)}{k_4^2} y_1^{(5)}\\[5pt] \tilde{\tau }_2 = \left (I_{21}^*(y_3)-I_{3}^*\right ) \ddot{y}_{2}+\left (I_{22}^*(y_3)-I_{3}^*\right ) \ddot{y}_{3}+b_2(y_3,\dot{y}_2)+\tilde{\tau }_{3} \\[5pt] \tilde{\tau }_1 = \left (I_{11}^*(y_3)-I_{21}^*(y_3)\right ) \ddot{y}_{2}+\left (I_{12}^*(y_3)-I_{22}^*(y_3)\right ) \ddot{y}_{3}+\left (b_1(y_3,\dot{y}_2,\dot{y}_3)-b_2(y_3,\dot{y}_2)\right )+\tilde{\tau }_{2}\\[5pt] \end{gathered} \end{equation}

The flat variables $\boldsymbol{y}$ are the same for both cases and are calculated by means of polynomial functions, in particular, an $11$ -th degree polynomial for $y_1$ and a $5$ -th degree polynomial for $y_2$ and $y_3$ . The polynomial coefficients are calculated by means of Eq. (31), and the initial and final conditions are defined as follows:

(36) \begin{equation} \begin{cases} y_1(t_i) = \sum _{k=1}^4 (q_k)_i \\[5pt] y_1(t_f) = \sum _{k=1}^4 (q_k)_f \\[5pt] y_2(t_i) = q_{1i}\\[5pt] y_2(t_f) = q_{1f}\\[5pt] y_3(t_i) = q_{2i}\\[5pt] y_3(t_f) = q_{2f}\\[5pt] \dot{y}_1(t_i) = \ddot{y}_1(t_i) = y_1^{(3)}(t_i) = y_1^{(4)}(t_i) = y_1^{(5)}(t_i) = \dot{y}_2(t_i) = \ddot{y}_2(t_i) = \dot{y}_3(t_i) = \ddot{y}_3(t_i) = 0\\[5pt] \dot{y}_1(t_f) = \ddot{y}_1(t_f) = y_1^{(3)}(t_f) = y_1^{(4)}(t_f) = y_1^{(5)}(t_f) = \dot{y}_2(t_f) = \ddot{y}_2(t_f) = \dot{y}_3(t_f) = \ddot{y}_3(t_f) = 0 \end{cases} \end{equation}

Table III. Parameters of the simulated $4$ -DOF robot.

Figure 7. Simulated joint variables of the $4$ -DOF robot considering the torque calculated without and with the viscous friction. An $11$ -th degree polynomial has been used for $y_1$ and a $5$ -th degree polynomial for $y_2$ and $y_3$ in both cases.

The simulations are performed considering a motion time equal $t_f=0.5 \,s$ . The start configuration is $\boldsymbol{q}_i=[\!-\!20, -20, -20, 0]\,deg$ , while the end configuration is $\boldsymbol{q}_f=[185,35,20,0]\,deg$ . The parameters of $4$ -DOF robot are summarized in Table III.

Figure 8. Simulated trajectories of the end-effector of the $4$ -DOF robot in $x-y$ plane considering the torque calculated without and with the viscous friction. An $11$ -th degree polynomial has been used for $y_1$ and a $5$ -th degree polynomial for $y_2$ and $y_3$ in both cases.

Fig. 7 depicts the joint variables of the robot. The behavior of the $4$ -DOF robot is very similar to the one of the $2$ -DOF robot in which motor torque has been calculated by means of $11$ -th degree polynomials (Fig. 6). The introduction of viscous friction has a direct effect on torque $\tilde{\tau }_3$ and an indirect effect on the other motor torques through the Newton-Euler recursive formula (Eq. (15)). As a result, $q_1$ and $q_2$ are poorly affected by viscous friction, $q_3$ is slightly affected by viscous friction, and $q_4$ is strongly affected by viscous friction, especially at the end of the motion, since $q_4$ shows oscillations about the final value if viscous friction is not considered in the calculation of torques.

Fig. 8 shows the simulated trajectories of the end-effector considering the torque calculated without and with viscous friction. The difference between the trajectories is relevant near the cusp and at the end of the motion. Moreover, the zoom of the end of the trajectories (depicted in Fig. 8) clearly shows that the end-effector presents oscillations at the end of the planned motion if viscous friction is not taken into account in torque calculation. Such a behavior is not suitable for many applications, in which the motion is considered complete when the amplitude of residual end-effector vibrations is reduced within a small range. Conversely, if torques are calculated taking into account the viscous friction, the system does not show residual vibrations.

These results show that large improvements in point-to-point motion can be achieved with a proper open-loop system without any feedback from the non-actuated joint.