## Abstract

An enhanced inverse dynamics approach is here presented for feedforward control of underactuated multibody systems, such as mechanisms or robots where the number of independent actuators is smaller than the number of degrees of freedom. The method exploits the concept of partitioning the independent coordinates into actuated and unactuated ones (through a QR-decomposition) and of linearly combined output, to obtain the internal dynamics of the nonminimum-phase system and then to stabilize it through proper output redefinition. Then, the exact algebraic model of the actuated sub-system is inverted, leading to the desired control forces with just minor approximations and no need for pre-actuation. The effectiveness of the proposed approach is assessed by three numerical test cases, by comparing it with some meaningful benchmarks taken from the literature. Finally, experimental verification through an underactuated robotic arm with two degrees of freedom is performed.

## 1 Introduction

### 1.1 Motivations and State of the Art.

Precise trajectory tracking is a relevant research topic in the field of mechanisms, multibody systems, and robotics. To achieve it, controllers exploiting both feedforward (open-loop) and feedback (closed-loop) actions are usually adopted: the first contribution is aimed at ensuring faster responses with lower tracking errors and low feedback gains (and hence higher robustness margins); the second one is aimed at compensating for model-plant mismatches, at changing the dynamic properties (e.g., the eigenstructure [1–3]) and at ensuring disturbance rejection. As for feedforward, its implementation requires the computation of the inverse dynamic model of the system under investigation. In the case of fully actuated mechanisms and robots, this computation is usually straightforward and just relies on simple algebraic calculations. On the contrary, when the number of independent actuators is smaller than the number of degrees of freedom (DOF), i.e., underactuated multibody systems are considered, the computation of inverse dynamics is challenging and often relies on the solution of both algebraic and differential equations. Underactuation usually arises because of the presence of passive joints [4–6] or flexible links [7–9]. The explicit, algebraic solution of the inverse dynamics problem is possible if the system is differentially flat [10]. In differentially flat systems, inverse dynamics is solved through an algebraic model, and the control forces are algebraically specified by the required output trajectory. Therefore, no pre- and post-actuations are required. Flatness is a system property that is not straightforward to assess; i.e., it requires a lot of algebraic manipulations and finding a flat output definition is not trivial [11], since no systematic approach exists. Additionally, exploiting flatness in underactuated multibody systems imposes very smooth reference trajectories, usually with four continuous derivatives [12]. For example, up to ninth-degree polynomials are usually adopted in the literature. In contrast, high-degree polynomials are not widely adopted in motion planning of mechatronic systems because of large requirements of speed, acceleration (both peak and average) and hence of torque, power, and energy that make these motion profiles often unfeasible when considering limits of the components (in particular, motor and gearbox [13–15]).

The difficulties in inverting the dynamic model of underactuated multibody systems are exacerbated for some definitions of desired output, such as in the case of noncolocation of actuators and controlled outputs [16], which leads to unstable internal dynamics. In the case of linear systems, for example, this is related to the presence of right half-plane zeros: by inverting the model to compute the suitable control forces to track the desired output trajectory, these zeros would become right half-plane poles, causing the divergence of the control actions. Generally speaking, with reference to both linear and nonlinear systems, the solution of the inverse dynamic problem of these so-called nonminimum-phase systems leads to unbounded solutions, unless noncausal inversion is adopted [17], or approximated causal solutions are sought [16].

To overcome these difficulties, and given the relevance of the topic in the field of motion control, several papers have been proposed in recent years. An inversion-based approach to exact nonlinear output tracking control is proposed in Ref. [17] where a feedback term is introduced to stabilize the system along the desired trajectory; the problem is solved by applying the Byrnes–Isidori regulator to a specific trajectory, introducing boundedness and integrability requirements. However, the solution is noncausal, i.e., the control forces begin before the starting of the trajectory, and therefore, this technique results to be not implementable through an online procedure. In Ref. [18], an extension of stable inversion to nonlinear time-varying systems is carried on, considering both minimum-phase and nonminimum-phase systems. This technique consists of linearizing the system dynamics and subsequently partitioning it into time-varying stable and unstable sub-systems, which are used to design time-varying filters that ensure a bounded inverse input-state trajectory. However, this approach is local to the time-varying path, and furthermore, it requires that internal dynamics varies slowly. An interesting analysis of the concept of approximated model inversion for nonminimum-phase systems is proposed in Ref. [19], with reference to linear systems. Three stable approximate model inversion techniques are examined: the nonminimum-phase zeros ignore (NPZ-Ignore), the zero-phase-error tracking controller (ZPETC), and the zero-magnitude-error tracking controller (ZMETC). Since the proper choice of one of these techniques highly depends on the system under investigation, this work also discusses how the position of the zeros can indicate the more effective technique to be used to maximize the performances. When dealing with nonlinear systems, an effective and simple approach is the output redefinition, i.e., the approximation of the desired output with a formulation that ensures a simpler problem solution and stable internal dynamics. This idea has been first introduced in the milestone paper [16] and then also exploited with a different formulation in Refs. [20] and [21]. In both the papers, an algebraic-differential scheme is adopted, where the differential part can be easily solved through a standard numerical scheme for integration of ODEs (ordinary differential equations), representing the internal dynamics. A remarkably different approach is proposed in Ref. [22], by handling dynamic models formulated either with ODEs or with DAEs (differential algebraic equations), and without requiring the explicit derivation of the internal dynamics model. In such a paper, inverse dynamics is formulated as an optimal control problem, and a noncausal solution is obtained.

### 1.2 Contributions of This Paper.

In this paper, a feedforward computation algorithm for nonminimum-phase underactuated multibody systems, modeled through ODEs, is proposed by exploiting and extending the idea of linearly combined output adopted in Refs. [16,20,21]. Compared to the usual formulation of such an approximation, a more general formulation of the linearly combined output is handled in this paper. The goal of such an output representation is to easily obtain the internal dynamics for stabilizing it and then to obtain a stable model inversion, without pre-actuation. The method includes a systematic approach to stabilize the unstable internal dynamics, by means of performing through proper output redefinition. The resulting differential scheme is solved through the forward dynamics of the stabilized internal dynamics, by exploiting any scheme for ODE numerical integration. Then, the approximated formulation of the controlled output is replaced by its actual nonlinear formulation (at the position, speed, and acceleration levels) to compute the commanded value of the actuated variables and to perform exact algebraic inversion of the dynamic model governing the sub-system made by the actuated coordinated. Since no output approximation is done in solving such an algebraic model, the proposed method provides an improvement compared to the existing methods: the approximation assumed through the linearly combined and redefined output is just used in the internal dynamics and hence causes small errors in a large variety of underactuated multibody and robotic systems.

Finally, to tackle different topologies of actuations, an approach for partitioning the model coordinates into actuated and unactuated ones is proposed, by exploiting the *QR-decomposition* of the input matrix, which allows for handling a large variety of underactuated systems.

The effectiveness of the method is proved through three test cases, aimed at highlighting different features of the proposed method and also by providing a fair comparison with several benchmarks.

Additionally, the experimental verification through a 2-DOF underactuated robotic arm is proposed, thus providing a further improvement of the literature which usually just provides numerical applications.

As a result, a comprehensive, and computationally simple, method for precise and accurate inverse dynamics is formulated, which can be applied to several underactuated systems, by improving the existing state of the art.

## 2 Method Description

### 2.1 Definitions.

*n*DOFs, formulated through a set of

*n*ODEs

$z\u2208Rn$ is the vector of the independent generalized coordinates, $Mz(z)\u2208Rn\xd7n$ is the mass matrix, $cz(z\u02d9,z)\u2208Rn$ contains the Coriolis terms, $Kz(z)\u2208Rn\xd7n$ is the stiffness matrix, and $fz(z\u02d9,z,t)\u2208Rn$ are the external forces projected on the directions of $z$. *t* denotes the time; along the paper, the dependence on *t* is usually not written, unless the cases where it makes the comprehension clearer. The system is controlled through *m* < *n* independent control forces collected in vector $u\u2208Rm$, and $Bz\u2208Rn\xd7m$ is the force distribution matrix, which is assumed to be full column-rank; such matrix can be dense since this paper will propose a technique to handle arbitrary topologies of $Bz$.

Since $rank(Bz)=m<n$, the system is said to be underactuated. In underactuated systems, the dimension of the configuration space is greater than the one of the control input space.

The goal of inverse dynamics is to compute the time-history of the control input $u(t)$ to make $y(t)$ track the desired output trajectory $ydes(t)\u2208Rm$ as closely as possible. It is assumed that $ydes(t)$ is at least twice differentiable, as usually done in motion planning of flexible systems ([23,24]).

It is also assumed that the system's initial conditions $z(0)$ and $z\u02d9(0)$ are known, as well as the external forces $fz$.

In fully actuated systems, i.e., when *n* = *m*, model inversion is straightforward and just relies on simple algebraic calculations by taking advantage of the inversion of matrix $Bz$. In contrast, in the case of underactuated systems, $Bz$ is rectangular and cannot be inverted; therefore, model inversion should rely on more complicated solutions, as the one discussed in the following of this paper.

### 2.2 Coordinate Partitioning.

*m*not-null entries, the model in Eq. (1) does not allow for immediate partition of the coordinates into actuated and unactuated ones (e.g., coordinate permutation is not an effective solution). A coordinate transformation should be therefore performed by exploiting the

*QR-decomposition*of $Bz$. It yields an orthonormal matrix $Q\u2208Rn\xd7n$ ($QTQ=QQT=I$, where $I$ denotes the identity matrix of proper dimension) and an upper-triangular matrix $B\u2208Rn\xd7m$. $B$ is upper-triangular; hence, it can be written through its sub-matrix $BA\u2208Rm\xd7m$ which collects all its non-null entries and, moreover, is nonsingular, i.e., $rank(BA)=rank(Bz)=m$

**0**denotes a null matrix of proper dimension. A new vector of the generalized coordinates $q\u2208Rn$ is introduced, such that

Although the entries of $q$ have sometimes no straightforward physical interpretation, this transformation partitions the vector of the generalized coordinates into *m* actuated coordinates, $qA\u2208Rm$, and *n* − *m* unactuated coordinates, $qU\u2208Rn\u2212m$.

In the following, it is assumed that **Q** is constant; otherwise, the definition of $q$ should change during the motion or in time. In contrast, $BA$ can depend on **q**: this is the occurrence solved in this paper, which is a common situation in several multibody and robotic systems (as shown through the three examples). Nonetheless, this technique provides a more general formulation than the papers in the literature, such as Refs. [16,20,21], which just handles models that are already partitioned.

_{,}and $f$ have the following dimensions: $MAA,KAA\u2208Rm\xd7m$, $MAU,KAU\u2208Rm\xd7(n\u2212m)$, $MAUT,KAUT\u2208R(n\u2212m)\xd7m$, $MUU,KUU\u2208R(n\u2212m)\xd7(n\u2212m)$, $cA,fA\u2208Rm$, and $cU,fU\u2208R(n\u2212m)$. The output equation in Eq. (2) is modified as well, by expressing

**y**as a function of the new coordinate vector

**q**through function $h:Rn\u21a6Rm$

*m*equations of the model in Eq. (6) represent the dynamics of the actuated sub-system, which is directly affected by the control forces $u$ through the input matrix $BA$ and by the external forces $fA$. Conversely, the remaining

*n*−

*m*equations represent the dynamics of the unactuated coordinates that are forced by the external disturbance forces $fU$ and are indirectly controlled by the evolution of the actuated coordinates

Equation (8) is a set of second-order nonholonomic constraints [25] and hence cannot be used to reduce the system model dimension by removing some generalized coordinates. In order to be feasible, any required trajectory should satisfy such constraints at any time. For this reason, the number of independent outputs has been assumed to be equal to the number of independent control inputs.

### 2.3 Approximation of the Internal Dynamics Through Linearly Combined Output.

This use of linear approximations of the output map accomplishes several tasks: on the one hand, it easily allows formulating the inverse dynamics; on the other one, it can be exploited to stabilize the internal dynamics whenever it is unstable. The tracking errors that could arise because of such an approximation are treated, in practice, as those due to the unavoidable mismatches between the model and the dynamics of the actual plant, and hence compensated by the necessary feedback control schemes. Indeed, as mentioned in the introduction, feedforward control should be supported by feedback control if precise motion is required.

On the other hand, Eq. (9) might degrade performances in some particular systems where such an approximation is rough; nonetheless, the use of local approximations, such as piecewise-linear models, overcomes this issue.

The term in Eq. (11) is not constant since $MUU=MUU(ydes,qU)$, $cU=cU(y\u02d9des,q\u02d9U,ydes,qU)$, $fU=fU(y\u02d9des,q\u02d9U,ydes,qU,t)$, $MAUT=MAUT(ydes,qU)$, and $KAUT=KAUT(ydes,qU)$ (these dependencies have been omitted in Eq. (11) for brevity of representation).

The stability of the internal dynamics in Eq. (11) strongly depends on the chosen output $y$. Usually, the presence of noncolocated pairs of control forces and desired output leads to zeros of the input–output normal form having a positive real part.

The ODEs in Eq. (13) can be solved by forward integration to compute $q\xa8U(t)$, $q\u02d9U(t)$, and $qU(t)$, i.e., the time-history of unactuated coordinates in the execution of the desired motion profile $y\xa8des(t)$, in the presence of the external forces in $fU(t)$. The selection of the discretization time-step and of the numerical scheme adopted for integration relies on the usual rules adopted in the numerical simulations of multibody systems [26].

### 2.4 Solution of the Algebraic Equation for the Actuated Coordinates.

*m*algebraic equations representing the inverse dynamics problem of the actuated sub-system. Rather than using again the approximations of the desired output provided by Eqs. (9), (10), and (12), this second step can be performed by inverting the actual map $ydes=h(qA,qU)$, together with its first- and second-time derivatives $y\u02d9des=H(q)q\u02d9$ and $y\xa8des=H(qA,qU)q\xa8+$$H\u02d9(qA,qU,q\u02d9A,q\u02d9U)q\u02d9$, to compute the reference trajectory for the actuated coordinates $qAdes(t)$, $q\u02d9Ades(t)$, $q\xa8Ades(t)$, where $H(qA,qU)\u2208Rm\xd7n$ is the Jacobian matrix relating the desired output with the actuated and unactuated coordinates, while $H\u02d9(qA,qU,q\u02d9A,q\u02d9U)\u2208Rm\xd7n$ indicates its time derivative. This original approach proposed in this paper partially compensates for the approximation introduced by the linear approximation, as well as by its stabilizing form in Eq. (12). By defining as $g:Rn\u21a6Rm$ such an inverse map, the motion of the actuated coordinates at the position level is inferred as follows:

**G**is the Jacobian matrix of

**g**, i.e., $q\u02d9Ades=$$\u2202g(ydes,qU)\u2202ydesy\u02d9des+\u2202g(ydes,qU)\u2202qUq\u02d9U$.

*m*equations in Eq. (6), which represent the forced dynamics of the actuated sub-system, the required control forces can be inferred from the following algebraic computation, which provides exact inversion of the actuated part of the dynamic model

Besides computing the feedforward control force, the proposed method can be exploited to generate consistent references to the actuated variables $(qAdes(t),q\u02d9Ades(t))$ that can be fed to the feedback control loop (or loops in case of decoupled feedback controllers) to compute the closed-loop corrections required to compensate for the unavoidable model mismatches.

## 3 Test Case 1: Flexible Beam With Lumped Masses

### 3.1 Description and System Model.

The first test case is a linear system composed of a flexible beam, modeled through the finite-element approach, plus some lumped spring-damper-mass systems, that recall the model of linear vibratory feeders often adopted in manufacturing plants to convey small components [27]. Despite the apparent simplicity of the system, due to its linearity, the proposed test case has several challenging issues to be properly tackled.

*EJ*= 1.93 × 10

^{5}N · m

^{2}, and linear mass density

*ρA*= 22.87 kg · m

^{−1}). The beam is connected to the frame through two springs (

*k*

_{l}=

*k*

_{r}= 1.5 × 10

^{5}N · m

^{−1}). Three actuators exert three pairs of independent forces, as shown in Fig. 2, and are modeled through three lumped masses (

*m*

_{ai}= 23 kg, with

*i*= 1, …, 3) connected to the beam by means of three linear springs (with stiffness

*k*

_{ai}= 4.6

*e*5 N · m

^{−1}, with

*i*= 1, …, 3). The resulting model has 13 DOFs, which include five absolute transversal translations (

*η*

_{j}, with

*j*= 1, …, 5), five absolute rotations (

*φ*

_{j}, with

*j*= 1, …, 5) of the beam, and three absolute translations of the actuators $(\eta a1,\eta a2,\eta a3)$

*QR-decomposition*of $Bz=QB$, leading to

Matrix $QT\u2208R13\xd713$ in Eq. (22), which maps **z** into $q=[qAqU]T$, with $qA\u2208R3$ and $qU\u2208R10$, reveals that the actuated coordinates $qA$ represent the relative translations of the nodes where each actuator is placed. In contrast, many coordinates in $qU$ have not a clear physical meaning.

Since $ydes$ is a subset of $z$, by considering the linear state transformation due to the *QR-decomposition*, $z=Qq$, it can be easily inferred that $\Gamma A\u2208R3\xd73$ and $\Gamma U\u2208R3\xd710$ are proper sub-matrices of $Q$. Hence, the use of linearly combined output is an exact representation of the test case under consideration.

Inverse dynamics is exploited to compute the actuator forces that make $y$ track a periodic, non-sinusoidal, asymmetric reference waveform, as sometimes adopted in feeders [29]. Each period of the wave consists of a rising phase lasting 25 ms, followed by a drop phase of 5 ms, and the absolute value of the peak displacement is 5 mm. The period of oscillation is 35 Hz, as often adopted in feeding small parts ([28,30]). The position references of the three controlled output are set with a $40deg$ phase shift.

### 3.2 Stabilization of the Internal Dynamics.

The eigenvalue analysis of the internal dynamics reveals that it is unstable, as it can be inferred from Fig. 3, because of the presence of poles with positive real parts.

The literature sometimes suggests a simple approach to stabilize the nonminimum-phase system by removing from the linear model the nonminimum-phase zeros, and hence the unstable poles of the internal dynamics, leading to the so-called “nonminimum-phase zeros ignore” (NPZ-ignore) approach [19]. Although this option could seem effective at a first glance because the unstable poles are quite far from the reference harmonic components, its application leads to large errors. The results of the application of the NPZ-ignore approach are therefore omitted to keep the paper concise.

Therefore, the stabilization approach by exploiting approximated values of $\Gamma A$ and $\Gamma U$ in Eq. (12), hereafter denoted $\Gamma ~A$ and $\Gamma ~U$, is adopted. Such matrices should be wisely chosen to ensure asymptotic stability with the negligible spillover on the other poles, while allowing for a precise approximation of the tracked output. To simplify the analysis, $\Gamma A$ has been kept constant while the approximation of $\Gamma U$ is parametrized through just a scalar *α*, as $\Gamma ~U=\alpha \Gamma U$. Figure 4 shows that *α* = 1 is the boundary condition between stable and unstable internal dynamics.

### 3.3 Numerical Results.

To prevent instability arising from numerical integration error, it has been assumed *α* = 0.99. The ODEs in Eq. (25) have been integrated through the fourth-order Runge–Kutta scheme with a sampling time of 10 *μ*s, due to the presence of some high-frequency poles. The time-history of the control forces $u$ is reported in Fig. 5, while the tracking responses are reported in Fig. 6 together with the corresponding tracking errors. The maximum tracking error is 5.2 × 10^{−2} mm, that is the $1.0%$ of the amplitude of the periodic motion; the average RMS (root mean square) error is instead 2.1 × 10^{−2} mm. Such errors are due to the output redefinition to stabilize the internal dynamics; additionally, the presence of underdamped dominant dynamics in the system model exacerbates such errors, by creating some residual, although small, oscillations after motion completion.

It should be noted that no pre-actuation is required. In contrast, post-actuation is required after the end of the required motion, with a small force to be exerted by the actuators (and a very small position error that quickly settles) since the system is not differentially flat.

## 4 Test Case 2: Overhead Cartesian Crane

### 4.1 Description and System Model.

*θ*

_{X}is the oscillation angle projected on the

*XZ*plane, while

*θ*

_{Y}is the angle measured from the

*XZ*plane. By assuming that the cable is taut with tensile stress and that hoisting is not allowed (i.e., the cable length

*h*is constant), the following four nonlinear equations describe the system's dynamic behavior [31]

The following definitions have been introduced in Eq. (27): *u*_{X}, *u*_{Y} are the forces that are applied to the platform along the *X* and *Y* directions; *m* is the lumped mass of the suspended load; $cX,cY,c\theta $ are the viscous friction coefficients; *M*_{X}, *M*_{Y} are the translational masses of the platform along the *X* and *Y* directions respectively; *g* is the gravity acceleration. The values of all the parameters of the system model are reported in Table 1.

Parameter | Value |
---|---|

M_{X}, M_{Y} | 30 kg |

m | 0.7 kg |

h | 1 m |

c_{X}, c_{Y} | 0.5 N · s · m^{−1} |

$c\theta $ | 0.25 N · m · s · rad^{−1} |

Parameter | Value |
---|---|

M_{X}, M_{Y} | 30 kg |

m | 0.7 kg |

h | 1 m |

c_{X}, c_{Y} | 0.5 N · s · m^{−1} |

$c\theta $ | 0.25 N · m · s · rad^{−1} |

In this example, both the dynamic model and the output map are nonlinear. The approximation of $y$ as a linear combination of $q=[qATqUT]T$ is obtained by linearizing Eq. (28) about the vertical equilibrium configuration (i.e., $\theta X=\theta Y=0deg$), leading to $\Gamma A=I\u2208R2\xd72$ and $\Gamma U=hI\u2208R2\xd72$.

### 4.2 Stabilization of the Internal Dynamics.

**0**denotes a null matrix).

To study the stability of the internal dynamics, the equations of motion of the unactuated sub-system $(MAUTq\xa8A+MUUq\xa8U+cU=fU)$ are linearized around the stable downward equilibrium point defined by $\theta X=\theta Y=0deg$, and secondly, $q\xa8A$ is replaced thanks to the linearly combined output technique described in Sec. 2.3, with $y\xa8des=q\xa8A+\Gamma Uq\xa8U$.

If the exact $\Gamma U$ is used, and therefore *α* = 1, $MID$ becomes a null matrix and therefore the internal dynamics features poles with infinite positive values. To assess the effect of α on the stability, Eq. (30) can be treated as two independent scalar equations by considering that $MID$, $CID$, and $KID$ are diagonal matrices, and the two characteristic polynomials are separately discussed. Since all the entries of $CID$ and $KID$ are positive, Descartes’ Rule of Signs reveals that the internal dynamics is stable if and only if *α* < 1. As in the previous test case, the parameter *α* should be, therefore, chosen as close as possible to the unitary value, without exceeding it, ensuring a good reconstruction of the desired output and making negligible the pole spillover, while accounting for numerical issues of the ODE solver that can make the system unstable even with *α* < 1. In the following, it has been chosen *α* = 0.99 that leads to the following poles of the internal dynamics: *λ*_{1} = *λ*_{2} = −17.86 + 25.73*i* and *λ*_{3} = *λ*_{4} = −17.86 − 25.73*i*.

### 4.3 Numerical Results.

The test consists of a planar circle with a diameter of 0.5 m in the *XY* plane. The entire circular path is set to be completed in 10 s through a fifth-degree polynomial timing law; to assess the behavior of the load at steady-state conditions, two idle intervals of 4 s are introduced both at the beginning and at the end of the trajectory. The ODEs in Eq. (25) have been integrated through the fourth-order Runge–Kutta scheme with a sampling time of 10 *μ*s.

The displacements of the load and the trolley are shown in Fig. 8, which also provides details on the tracking errors: the reference trajectory is tracked with high precision, with a maximum transient error of 0.11 mm, and with a negligible residual oscillation after motion completion.

The path-tracking response is shown in Fig. 9. Figure 9(a) compares the motion of the trolley with the one of the load, by clearly showing that the proposed technique properly shapes the trolley path (by means of the control forces applied to it), by considering that centrifugal effects will make the load path differ from the one of the load. Figure 9(b) details the contour error of the load, with respect to the reference path: the maximum value is 0.12 mm, the RMS error during the motion is 0.04 mm, while the residual error after motion completion is negligible. Finally, Fig. 10 displays the computed forces applied to perform the desired task.

A final assessment of the capability of the proposed method is obtained by comparing it with the method presented in Ref. [20], by adopting the same stabilized linearly combined output. The results are summarized in Fig. 11, which shows the tracking and the contour errors sported by the two different ways of calculating $uFFW$, by adopting the same stabilization scheme for the internal dynamics. It is evident the benefit of the approach proposed in this paper, which remarkably reduces the errors: using Ref. [20] leads to a maximum tracking error equal to 0.35 mm, and a maximum contour error equal to 0.35 mm as well. Indeed, due to the centrifugal forces, the mismatch between the nonlinear model and the linearized one becomes relevant.

## 5 Test Case 3: Underactuated Robotic Arm With Flexible Passive Joint

### 5.1 Description and System Model.

The third test case is a planar flexible robotic arm lying in the vertical plane, as shown in Fig. 12, made by two rigid links connected through a passive revolute joint (B) with a linear torsional spring. Another joint (A) is directly actuated by a DC (direct current) motor. The actuated link is of length *l*_{1}, mass *m*_{1}, and moment of inertia *J*_{1} (defined with respect to point A). The moment of inertia of the motor, *J*_{m}, and of the rigid coupling, *J*_{d}, should be accounted for as well. The other link is of length *l*_{2}, mass *m*_{2}, and a moment of inertia *J*_{2} (defined with respect to point B). A nodal mass *m*_{B} is placed at joint B to include the mass contribution of the rotational encoder. Two linear springs couple the two links, resulting in an equivalent rotational spring with rotational stiffness *k*_{s}. All the values of the mentioned parameters are reported in Table 2.

Parameter | Value | Parameter | Value |
---|---|---|---|

l_{1} | 0.17 m | l_{2} | 0.155 m |

m_{1} | 0.05 kg | m_{2} | 0.021 kg |

J_{1} | 4.82 × 10^{−4} kg · m^{2} | J_{2} | 1.68 × 10^{−4} kg · m^{2} |

J_{m} | 2.7 × 10^{−5} kg · m^{2} | J_{d} | 2.3 × 10^{−5} kg · m^{2} |

m_{B} | 0.178 kg | k_{s} | 0.177 N · m/rad |

Parameter | Value | Parameter | Value |
---|---|---|---|

l_{1} | 0.17 m | l_{2} | 0.155 m |

m_{1} | 0.05 kg | m_{2} | 0.021 kg |

J_{1} | 4.82 × 10^{−4} kg · m^{2} | J_{2} | 1.68 × 10^{−4} kg · m^{2} |

J_{m} | 2.7 × 10^{−5} kg · m^{2} | J_{d} | 2.3 × 10^{−5} kg · m^{2} |

m_{B} | 0.178 kg | k_{s} | 0.177 N · m/rad |

The computed torque is commanded by translating it as the reference current through the nominal torque constant of the motor. Such a reference current is fed to the drive current loop through Simulink Real Time scheme that communicates the current to be delivered at each instant to the motor through an external input–output board. It should be noted that the current loop is the sole closed-loop controller adopted in these experiments: neither the position nor the speed feedback loops are used, since trajectory tracking just relies on an open-loop scheme to assess the effectiveness of the proposed approach.

*u*) with respect to the vertical configuration, and the relative position of the passive link with respect to the actuated one, denoted

*θ*

_{1}and

*θ*

_{2}, respectively, as shown in Fig. 12. The following nonlinear model of the system is formulated:

### 5.2 Stabilization of the Internal Dynamics.

*l*

_{1}+

*l*

_{2}, linearized around the stable downward equilibrium point defined by

*θ*

_{1}=

*θ*

_{2}= 0

*M*

_{ID},

*C*

_{ID},

*K*

_{ID},

*B*

_{ID}∈ ℝ assume the following meaning:

*λ*

_{1,2}∈ ℂ of the homogeneous characteristic equation associated with Eq. (38) are calculated

To provide a similar treatment as the previous test cases, defining $\Gamma ~U=\alpha \Gamma lim$ (with *α* ∈ ℝ ) and $\Gamma lim=2l23l1+2l2$ recasting the inequality in (41) as *α* < 1. The same result can be obtained by considering the stability features of the associated zero-dynamics [16], which is unstable in nonminimum-phase systems.

As previously seen, the parameter *α* should be chosen as close as possible to the unitary value, without exceeding it, ensuring a good reconstruction of the desired output and, at the same time, making the spillover effect of the poles negligible.

### 5.3 Numerical Application.

The test consists of a pick-and-place repetitive cycle, from the pick-position to the place one, and then back again to the pick-location. Indeed, pick-and-place tasks have great importance for robots employed in manufacturing plants [32,33], and hence, the end-effector motion should be properly controlled. Two different amplitudes of the motion have been considered, corresponding to two different opening angles of the equivalent rigid single-link mechanism equal to $30deg$ and $50deg$. The motion time of each displacement is 1.2 s, and fifth-degree polynomial timing laws have been assumed. Additionally, rest intervals, lasting 1 s for the $30deg$ case and 2 s for the $50deg$ one, are required between each motion.The method application has led to the torques shown in Fig. 14, which are continuous and do not require pre-actuation. The ODEs in Eq. (25) have been integrated through the fourth-order Runge–Kutta scheme with a sampling time of 1 ms.

Figure 15 shows the result obtained, by applying the computed torque to the simulator of the nonlinear mechanism: precise tracking is obtained for the three references, and just some small errors are present since references and actual displacement are almost overlapped.

A clearer view is provided by the analysis of the tracking error for both the tests and through benchmarking: the outcomes of the proposed method have been compared with those sported by the methods proposed in Refs. [16] and [34], which are assumed as the benchmarks for this test case. In the application of Ref. [16], the same output redefinition is adopted to stabilize the internal dynamics. The results are shown in Fig. 16, where the tracking errors are detailed to corroborate the effectiveness of the method proposed in this paper and the advancement it provides compared with the state-of-the-art methods. In the case of an opening angle equal to $30deg$, the proposed method leads to a maximum error equal to 1.3 mm and to an RMS error equal to 0.45 mm ; the benchmark I, i.e., the method proposed in Ref. [16], leads to a maximum error equal to 2.1 mm ($+62%$) and to an RMS error equal to 0.9 mm ($+100%$). Even larger errors are sported by benchmark II, i.e., the method proposed in Ref. [34]: the maximum error is equal to 17.6 mm ($+1254%$) and the RMS error is 7.1 mm ($+1478%$), with visible residual oscillations.

In the presence of the $50deg$ opening angle, the proposed method still outperforms the benchmarks since it leads to a maximum error equal to 2.5 mm and to an RMS error equal to 0.8 mm ; benchmark I leads to a maximum error equal to 3 mm $(+20%)$ and to an RMS error equal to 1.2 mm $(+50%)$; again, remarkably larger errors are sported by benchmark II: the maximum error is equal to 54.2 mm $(+2068%)$, and the RMS error is 28.1 mm $(+3412%)$.

Finally, the comparison of the computed torques is provided in Fig. 17: very similar values are provided since all these methods are aimed at best approximating the exact torques; in particular, benchmark I, that is the most accurate one between the two benchmarks, is almost overlapped with the proposed method.

### 5.4 Robustness Considerations Against Model Uncertainties.

Since dynamic models are always affected by uncertainties with respect to the real setups, the robustness of the results achieved by the proposed method is assessed in the following by applying some perturbations to the dynamic model emulating the real system, to which feedforward torques computed through the nominal model are applied. The uncertain parameters of the dynamic models are the motor moment of inertia, the masses and lengths of the links (and hence the related moment of inertia), the stiffness of the rotational spring, and the entries of the damping matrix. First, the decoupled effect of the uncertainty of just one parameter at a time, by keeping the other ones to their nominal values, has been evaluated to assess which is the most sensitive parameter in terms of trajectory tracking error. The results of this analysis are reported in Fig. 18, for both the opening angles considered (30 deg and 50 deg), where ±10% variations of the nominal parameters are assumed. It can be noticed that the lengths of the two links are the most influential parameter; this is due mainly to the error in the kinematics (Eq. (35)), which would also severely affect the inverse dynamic in a rigid link system. Among the two opening angles, smaller errors are obtained in the case of 30 deg, because of the smaller accelerations.

A second robustness analysis is done with respect to the uncertainty of the motor torque loop dynamics. In this paper, as well as in all the papers on inverse dynamics of multibody and robotic systems quoted in the references, model inversion is done by assuming the ideal response of the actuator current loop, i.e., infinite bandwidth of such a loop. This assumption is often reasonable since modern servomotors can reach up to 5000 − Hz bandwidth [35]. On the other hand, “cheap” or old motors might have a slower dynamic response. To evaluate this issue, robustness analysis is also performed by simulating the system response in the presence of a current loop modeled as a first-order system with a bandwidth ranging from 500 Hz to 5000 Hz. The results are summarized in Fig. 19 that shows the RMS error as a function of the bandwidth (by assuming no uncertainty on the other parameters). It is evident that the tracking error approaches the ideal one as the bandwidth increases, and for the case under investigation, bandwidth greater than 1500 Hz just causes a negligible increase. In contrast, the harmonic distortion and the lag in reproducing the commanded torque profile remarkably amplify the tracking error in the case of 500 − Hz bandwidth.

To further assess the robustness, a Monte Carlo analysis (with 2000 random tests) has been performed by considering the simultaneous perturbations of all the above-mentioned parameters. As for the current loop bandwidth, a uniform distribution in the interval (500 Hz, 5000 Hz ) has been assumed. As for the other parameters, Gaussian distributions have been assumed, with a standard deviation equal to the 5% of the nominal value. The average RMS value of the tracking error over a cycle is equal to 1.0 mm when the test case with the opening angle equal to $30deg$ is considered, while the maximum RMS error achieved in such an analysis is 3.0 mm. In the case of a 50 deg opening angle, the average RMS value is 1.8 mm while the worst-case simulated leads to 5.6 mm. The Monte Carlo analysis is summarized through the trajectory tracking responses reported in Fig. 20, to highlight the bands of uncertainty for both the opening angles evaluated, when the proposed method is adopted.

The same perturbations have been applied to Benchmark I, by obtaining worse results because of the higher approximation it exploits in model inversion: the average and maximum RMS errors are, respectively, 1.6 mm (+60% compared with the proposed method) and 3.8 mm (+27%) with the 30 deg opening angle, and 2.5 mm (+38%) and 5.9 mm (+6%).

### 5.5 Experimental Application.

The motor torques shown in Fig. 14, and computed through the method proposed in this paper, by exploiting the nominal model parameters, have been applied to the experimental setup. As far as friction is concerned, the viscous term is considered in Eq. (32), and hence, it affects also the internal dynamics. Additionally, a Coulomb term related to the friction at the active joint (and hence including the contribution of the motor rotor) and equal to ±0.01 Nm has been included in the algebraic Eq. (17) by adding the term $0.01sign(\theta \u02d91)$ Nm. Such a value has been estimated through experimental identification. A more complicated friction model can be employed in the algebraic part, by taking advantage of the wide literature on friction modeling and compensation; this goes, however, beyond the scope of the paper. No feedback terms are adopted, and therefore, the showed results just rely on open-loop control through model inversion, thus providing a fair and severe evaluation of its effectiveness.

The achieved displacements are reported in Fig. 21, for each opening angle, and have been measured through a high-resolution Vicon motion capture exploiting four cameras Vicon Vantage V5. A precise tracking of the references is obtained, although the tracking error is greater than the one expected by the numerical analysis: the RMS values of the error are equal to 3.5 mm with the 30 deg opening angle, and 4.3 mm with the 50 deg one. This aspect is mainly due to the unavoidable presence of model mismatches between the nominal one, used to compute the torque, and the real setup. In particular, among the parameters investigated through the Monte Carlo analysis, it should be noted that the current loop bandwidth of the actuator is about 500 Hz, thus approaching the worst-case scenario for such a parameter. A second relevant uncertainty source is friction, which has been remarkably simplified in the model and whose compensation has been performed in a very simple algebraic way. Finally, the deformation of the wire of the encoder placed at the passive joint significantly perturbs the system dynamics both in terms of mass distribution, by creating an equivalent stiffness (that slightly changes as the wire moves, since the diameter is 5 mm) and by acting as an external disturbance force acting on the system itself. Nonetheless, by considering that no feedback is adopted, tracking is accurate if compared with the resolution sported by the two encoders mounted on the mechanism, which is equal to 4.1 mm.

## 6 Conclusions

This paper proposes enhanced inverse dynamics techniques for non-differentially flat, nonminimum-phase underactuated multibody systems, to obtain a feedforward control term without pre-actuation. The method is suitable for mechanisms and robots with passive joints or flexible links and allows boosting the design of closed-loop control schemes with high robustness, thanks to the possibility of lowering the feedback gains.

By taking advantage of the idea already proposed in the literature, of properly approximating the desired output as a linearly combined output, and by performing the output redefinition to stabilize the internal dynamics, stable inversion of such a dynamic is obtained. Then, by exploiting the exact model of the actuated sub-system, the computation of the feedforward forces is obtained with just a minor approximation. This possibility is allowed by a general approach that exploits the *QR-decomposition* of the input matrix to partition the independent coordinates into actuated and unactuated ones, thus simplifying the definition of the internal dynamics and hence its stabilization.

The concurrent use of these tools results in a comprehensive, and computationally simple, method for precise and accurate inverse dynamics, which can be applied to several underactuated systems where linear (or piecewise-linear) output approximations are reasonable, by improving the existing state of the art.

The effectiveness of the proposed method has been assessed numerically with three different test cases, highlighting different features of the method also by means of comparison with several state-of-the art techniques. In all the studied cases, the proposed method has shown higher accuracy in tracking the desired trajectory with tracking errors that are smaller than all the benchmarks. Furthermore, the proposed technique has been experimentally applied on a 2-DOF underactuated robotic arm, made by two rigid links and a passive joint. The experimental results highlight the enhanced performances achieved through the proposed method even in the presence of plant uncertainties and model mismatches. As a result, the mean experimental tracking errors lay inside the uncertainty band defined by the resolution of the encoders mounted on the mechanism.

The absence of pre-actuation and the simple proposed formulation, which exploits some algebraic calculations and the numerical integration of the equation of motions of the stabilized internal dynamics without the need of optimization procedures, make the method suitable for applications where no pre-calculation is allowed and, therefore, where a real-time calculation of feedforward is required.

## Funding Data

This research was funded by the Italian Ministry of University and Research through the research grant “SISTEMA—Dipartimenti di Eccellenza” and by the PRIN 2020 “Extending Robotic Manipulation Capabilities by Cooperative Mobile and Flexible Multi-Robot Systems (Co-MiR)” (Prot. 2020CMEFPK).

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The data sets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.