Robust transient oscillation reduction for rest-to-rest motion of underactuated multibody systems

Conventional model-based design methods are often limited in their effectiveness by model-plant discrepancy. A solution to this problem is proposed in this work to enhance the robustness of motion planning solution for systems affected by parametric uncertainty. The method exploits a variational formulation in the form of a two-point boundary value problem (TPBVP) in which the robustness is achieved as a constraint enforced at the two boundaries. The formulation, which is specifically targeted at underactuated systems, aims at reducing both the transient and residual vibrations, as well as at mitigating the actuation effort. The development of the method is supported by its application to two numerical test-cases in the form of a double pendulum on a cart and a translating flexible beam.


Introduction
The execution of a precise and repeatable motion of underactuated systems while confining both transient and residual oscillations represents a classic challenge of mechatronics. This problem can be dealt with either as a control design problem or as a trajectory design one. Formulating it as the latter is rather convenient in all situations in which the effectiveness of the control system is limited by practical issues such as reduced actuator bandwidth or by the impossibility of altering pre-programmed commercial control devices, just to cite two common occurrences.
The literature on motion planning is extremely large and deeply rooted into engineering practice and research, and hence a detailed analysis of such a large corpus of research is well outside the scope of this work. However, it is still useful to mark a distinction between model-free and model-based approaches.
Model-free approaches rely on the definition of trajectories/motion profiles independently of the details of the dynamics of the system to be applied to. Commonly, such meth-P. Boscariol paolo.boscariol@unipd.it 1 Dipartimento di Tecnica e Gestione dei sistemi industriali, Università degli Studi di Padova, Vicenza, Italy ods make use of such methods as interpolation [14,18,32], and vibration reduction is mainly achieved by enhancing the smoothness of the motion profile [12,17]. These methods have gained wide popularity mainly due to their simple implementation and for being suitable to a wide array of applications. Currently they represent the de-facto standard choice for industrial robots [6,21].
Model-based methods, as the name suggests, make an explicit use of the dynamic model of the system for which the trajectory is planned. In many cases such methods rely on translating a motion design problem into an optimal control one, enabling the use of numerous and well-established methods belonging to the mature literature on optimal control [2,31].
Model-based methods can be divided into direct and indirect trajectory optimization methods [30]. Whenever a method uses a parametrized description of the kinematic quantities of the system to transform an optimal control problem into a parameter optimization problem, that method is commonly referred to as a direct one. Hence a finite-dimensional optimization problem is solved in lieu of a more challenging one of infinite size, making use of one of the countless solvers available for the purpose. The space of possible solution is, within this method, reduced with the implicit acceptance of a trade-off between ease of solution and optimality.
Indirect trajectory optimization methods, on the other hand, are based on the solution of a 'pure' optimal control problem, usually resorting to the formulation and solution of a twopoint boundary value problem (TPBVP) using calculus of variations. A rather heterogeneous collection of examples can cite applications of this method to multi-agent systems [13], spacecrafts [34], electric vehicles [23], and industrial robots [15]. Indirect methods are often praised for their accuracy but are often plagued by a small radius of convergence [5] and the applicability to 'small size' problems. The latter is generally referred to as the 'curse of dimensionality' [3]. Nonetheless, they still find some use whenever such limitations can be circumvented. Another commonly recognized drawback of indirect methods, which is shared with all optimal control methods, is the generally limited robustness to model-plant mismatches. The literature on motion planning has been, over the years, less sensitive to this issue in comparison with the control literature as robust control techniques are, by orders of magnitude, more frequently investigated than robust motion design techniques.
One area in which the need for robustness has met the actual development of robust solutions is input shaping. The techniques based on input shaping exploit some clever filtering techniques to alter a predefined motion or force profile to obtain a specific target, usually the limitation of residual oscillations. The literature on the topic has quickly recognized the robustness issue, and several robust shapers have been introduced over the years [22,25,29,33]. The popularity of such methods is due to their several advantages such as the applicability to a large range of applications and their straightforward implementation, as well as their applicability in conjunction with standard motion planning techniques.
Shapers, both in the form of nominal and robust ones, are used in this work as a benchmark against which the proposed robust planning method is tested. The proposed method, which is an indirect one, focuses on the definition of motion profiles for underactuated systems. In particular, it is aimed at mitigating not only residual oscillations, but also transient ones with the addition of robustness constraints, which are included in the design problem through sensitivity functions.
The application to two numerical test-cases, namely a double pendulum on a cart and a translating beam, will showcase the noticeable improvement over classic shaping techniques.

System model formulation
The aim of this section is to establish a system dynamics formulation to be used for formulating and solving the robust motion design problem. Let us take into account an underactuated mechanical system with n degrees of freedom (DOFs), which is represented by the vector of n independent coordinates q. The dynamics of this system can be represented as a set of n differential equations in which M ∈ R n×n is the mass matrix, K ∈ R n is the vector of position-dependent forces, and G ∈ R n collects gyroscopic and centrifugal forces, as well as damping forces. The force distribution matrix B(q) ∈ R n×m weighs the effects of the external forces collected in F ∈ R m . According to the hypothesis of dealing with an underactuated system, m must be smaller than n. Accordingly, q can be split into a vector q m of actuated DOFs and a vector q u of the n − m underactuated DOFs, as in Equation (2) can be split into two contributions, one of which expresses the dynamics of the n − m underactuated coordinates as follows: Whenever B is a sparse matrix with more than n not-null entries, a formulation consistent with the one in Eq. (2) can be defined by a proper QR decomposition, as described in the work [4].
The formulation in Eq. (3) does not include the vector of actuation forces F, but it does include the vector of accelerations of the actuated ones, i.e.,q a . This testifies that the motion of unactuated coordinates can be determined directly by 'shaping' the motion of the actuated ones. By assuming this point of view, what is usually treated as a control problem can be viewed as a motion design problem. Hence, the focus can be set on just Eq. (3), and the effect of the actuation forces can be neglected. Therefore the motion of actuated DOFs alone can be planned to achieve the desired properties for the unactuated ones.
Some care must, however, be taken as the dynamics of the actuators and of the control loop system that ensures the tracking on the actuated DOFs q a cannot be always neglected. In this case the planned and the executed motion of the actuated variable are different from each other. To cope with this possibility, it can be assumed that the value of the actuated variable q a is related to its reference value q ref a by means of a generic function h, as in Once a specific function h is defined, it can be included in the formulation of the system dynamics: in which, again, the explicit dependence on the actuation force vector F is avoided. This approach brings several practical advantages, the main one being the relaxation of the need for a detailed and accurate modeling of the system dynamics, which is often dependent on complex and nonlinear effects such as those due to friction forces. A formulation which does not explicitly depend on exogenous forces is also inherently more robust as each additional parameter to be determined for modeling purposes brings an additional source of uncertainty in the determination of its value. Moreover, a strategy that revolves around the definition of a suitable reference profile is generally of more straightforward practical implementation on the vast majority of commercial/industrial devices, in which the control loop cannot be altered or fine-tuned to achieve the desired performance improvement.
As stated before, the relationship h in Eq. (4) is arbitrary. Here it is chosen to balance the complexity of representation with its physical meaningfulness by assuming that the control loop can achieve a tracking without steady-state errors and with a finite bandwidth, as in a system with a single dominant pole. Therefore the actuated degree of freedomq a (the cart acceleration in this case) is related to its reference signalq which can also be written in the time domain as follows: Using this relationship in Eq. (5) provides the description of the dynamics of the system, which includes the effect of the limited bandwidth of the controller and of the actuators using the time constant τ as follows: This model has, as the only exogenous input, the acceleration reference signalq ref a , which will be designed according to the methodology that will be outlined in the next section. This methodology requires a system representation in the form of a system of differential equations indicated as f(x, t, v) (either a linear or a nonlinear one) with the state vector x and the input vector v in the forṁ f(x, t, v) being the first-order form of Eq. (8) obtained from it by adding a suitable number of integrators:

Variational formulation of the robust motion design problem
Having defined a description of the dynamics of the system according to either a linear or nonlinear state space form, we can use a variation formulation of the motion design problem. In particular, the motion design problem can be formulated and solved as a two-point boundary value problem, in which the initial and final instants of the motion are captured by the 'left' and 'right' boundary conditions. Defining either partially or completely the two boundary conditions does not, however, allow the definition of a specific trajectory as still infinite choices are available. Among this infinite set of solutions forq ref a (t), the one to be found will be the one that minimizes the cost function The scalar cost function g can be chosen arbitrarily, but its design should reflect the desired performance associated with the trajectory under development. The minimization of the cost function in Eq. (11) must comply also with the dynamics of the system, otherwise the planned motion would be infeasible. Collecting all constraints together results in the following optimization problem: subject to Most commonly this kind of problem is solved using calculus of variations and Pontryagin's minimum principle (PMP) [24,26]. The problem in Eqs. (12)- (14), once solved, produces a trajectory which is optimal, feasible, and compliant with the prescribed boundary values. However, this is true also in practical terms only if the actual plant is perfectly described by the dynamic model used in Eq. (14). For this reason, the problem in Eqs. (12)- (14) can be referred to as a nominal trajectory design problem. Any mismatch between the two systems-the actual one and the one used for planning-is neglected in the formulation outlined so far.
A robust counterpart of the aforementioned problem can be formulated using the framework defined in the previous work of one of the authors [8]. The methodology, which is briefly recalled here, is well suited whenever one (or more) of the parameters included in the plant dynamics formulation is affected by unmeasured or unmodeled variations, or simply it cannot be estimated with the needed level of accuracy. If the uncertain parameter is referred to as μ, then its influence on the system dynamics can be represented by including it within the system first-order differential equation description as follows: If f is continuous in (x, t, μ) and is continuously differentiable w.r.t. x and μ for any value of (x, t, u, μ) in the time frame [t 0 , t], then the system response can be evaluated as Since our aim is to quantify the effects of the change of μ on the system response, we can compute the partial derivative of x w.r.t. μ as in S(t) is referred to as the vector of sensitivity functions of the system as it gauges the effect of any deviation of the uncertain parameter μ on the state vector x(t). S(t) has many entries as x(t) has times the number of uncertain parameters. Sensitivity functions can hence be treated as the elements of a 'supplementary' state vector, as their dynamics is described by a set of differential equations exactly as the original state vector does. For the same reason, some metric of the sensitivity functions and some constraints on them can be included in the trajectory design problem structured as in Eq. (12)- (14) after having recognized that reducing the magnitude of sensitivity functions has the effect of reducing the impact of the related uncertainties.
Although not shown here (the reader might refer to [19] for the details), it can be shown that sensitivity functions do capture the first-order effects of the alteration of the parameter μ to the system trajectory, as in μ 0 being the nominal value of μ is recalled here to highlight that reducing the magnitude of S(t) has the effects of reducing the impact of the mismatch between μ and μ 0 on the time evolution of the plant.
Here, sensitivities are altered by acting on the motion design, but a recent work has shown that a similar target can also be tackled by introducing some physical alterations to the plant by an inverse structural modification approach [9].
Since sensitivity functions are meant to be weighted and constrained within the robust trajectory design problem, an augmented system dynamics is to be defined by stacking x(t) and S(t) as the 'robust' state vector whose dynamics is condensed as follows: The resulting robust trajectory design problem, as suggested in [8], can therefore be cast as (20) subject to The function g r in Eq. (20) has the same role as g in Eq. (12), i.e., to define the cost function to be minimized.
When compared with its nominal counterpart, the problem above includes an additional constraint (24) and two additional boundary conditions (22), which force sensitivities to be null at the two most critical points of the trajectory, i.e., the initial and final one. The latter, in particular, is rather important as the optimally designed trajectory will bring the robustified system at t = t f to a configuration with zero sensitivity and hence with minimal final error x(t f , μ) − x(t f − μ 0 ), as can be shown using Eq. (17) when written for t = t f and for S(t f ) = 0. This feature translates, in a vibrating system, into a trajectory, in which residual oscillations due to parametric uncertainties are minimized, as will be shown in the next sections through two numerical test-cases.
There are several ways of solving the optimization problem in Eqs. (20)-(23); here, the classical method based on the enforcement of the PMP conditions is used. This requires the formulation of the Hamiltonian of the system: The scalar Hamiltonian function includes the vector of the Lagrangian multipliers λ that constrains the solution to the augmented system dynamics defined in Eq. (19). The three conditions are, according to the PMP: As an analytic solution for this kind of problems is limited to a few simple problems, it is far more common to seek for a numerical solution using a collocation method [7,11] or a shooting method [1,10] just to cite a couple of possible options. The explicit integration of the system dynamics of Eq. (23)- (24) is not strictly needed as the TPBVP solver takes care of enforcing the three PMP conditions of Eq. (26) with or without performing it as, for example, in spectral methods. Here, the standard Matlab's solver bvp5c, which is based on collocation [26], is used for computing all the trajectories defined and tested in the next sections. The solver parametrizes the system state and the control action as piecewise sequences of fifth-order segments with C 1 continuity by means of a collocation algorithm [20]. The time needed to compute the trajectories that will be presented and analyzed in Sects. 4 and 5 ranges between 1.27 s and 2.32 s: such results have been obtained by using a laptop running Windows 10 with an Intel i5-8265U CPU with 8 GB of RAM.
As a final remark, it should be stated that the proposed method focuses on the mitigation of the effects of parametric uncertainties, so it does not apply, at least directly, to problems with other sources (or representations) of uncertainties such as neglected dynamics or additive noise.

System description and dynamic model
The first test-case consists of the planar double pendulum shown in Fig. 1. Two suspended masses, m 1 and m 2 , are connected to a translating cart through two cables whose lengths are According to the formulation of Eq. (8), it is assumed that the control input, to be designed through the proposed method, is the cart acceleration reference signal (i.e., v = q ref a =ÿ ref cart ). The matrices of the dynamics models are the following ones: with the state vector x: The nominal values used in all tests are: L 1 = 0.47 m, L 2 = 0.391 m, m 1 = 0.192 kg, and m 2 = 0.203 kg, while the time constant τ is set to 0.002 s. The uncertain parameter is assumed to be L 1 , which is chosen as its perturbation does significantly affect both natural frequencies of the system ω 1 and ω 2 . Therefore the sensitivity functions, formulated according to Eq. (24), are partial derivatives made with reference to L 1 : The augmented state vector is therefore Its dynamics is captured by the 'robust' first-order system that captures also the dynamics of the sensitivity functions:

Problem definition
The performance evaluation of the proposed method is conducted through numerical simulation and the comparison with two benchmark input shaping methods: the negative zero vibration (NZV) shaper [28] and the zero vibration derivative (ZVD) shaper [27]. The goal of the motion planning is to design a motion profile for rest-to-rest motion with negligible residual and reduced transient oscillations, as well as with enhanced robustness. The requirement of zero residual vibrations for the nominal case is obtained by imposing, in the optimization problem, null values for ϑ 1 , ϑ 2 and for their time derivatives at final time t f . As for the transient behavior, it is required to reduce both the 'average' and the peak values of the relative displacement between the cart and the end-point mass while ensuring reasonably small cart accelerations. To satisfy these requirements, the cost function g to be used in Eqs. (12) and (20) has been defined as follows: Such a definition of g trades off between the reduction of the acceleration reference signal v(t), i.e., the control effort, and the load oscillation, which is weighed by the positive scalar weight β and by the scaling factor γ . The robust counterpart of the motion planning algorithm is set up with the same cost function, i.e., g r = g. The same formulation is adopted in defining the cost function for the robust design problem, i.e., g r = g. The use of the exponential emphasizes the peaks of the oscillation in the cost function, thus keeping them small. In practice, it approximates a kind of minmax optimization problem within the frame of the classical optimal control while allowing for a straightforward numerical solution. The possibility of defining different objective functions to accomplish various secondary goals, besides the usual requirement of zero residual vibration, is an important advantage of the proposed approach over input shaping. Indeed, input shaping motion planning techniques work by convolving an arbitrary reference signal with a sequence of impulses whose number, amplitudes, and times of application depend on the natural frequency and damping of the vibrational modes to control, and on the desired robustness. Other control objectives, such as reduction of the transient oscillations and of the control effort, should be obtained by a clever selection of the original reference signal to be convolved with the shaper.
The problem of reducing transient oscillations and control effort by means of input shaping is exacerbated in three cases. First, if the desired motion time is comparable with the semiperiod of the damped oscillation in the case of nonrobust shapers or with the period in the case of ZVD. Again, if increased desired robustness is needed, and in the presence of more than one vibrational modes to be controlled. Indeed, the convolution of the unshaped reference with the shaper impulses introduces delays that are higher as the damped frequency decreases. The number of such impulses is then higher for robust shapers and for multi-mode shapers. The latter, indeed, are obtained by cascading two or more shapers, each one tuned for a vibrational mode. Hence, the only way to retain the original motion time is to pre-compensate the overall delay by reducing the time duration of the reference signal prior to shaping. This feature can be a limiting factor when the overall time delay is comparable with the required motion time since the unshaped motion profiles have very high accelerations that cause high transient oscillations.
As far as the choice of shaper is concerned, the NZV shaper is here chosen for being the nominal shaper that introduces the minimum amount of overall time delay, while the ZVD shaper is chosen as the standard robust one. As far as the choice of the unshaped motion profile is concerned, both shapers are used together with either a harmonic motion profile or a fifth-degree polynomial profile (henceforth denoted poly5 for brevity). The first one is chosen for being one of the standard profiles with lesser peak and RMS (root mean square) acceleration, a feature that reduces transient oscillations. The poly5 profile is chosen since it provides a smooth motion thanks to its continuous acceleration. The same design specifications, i.e., an overall displacement of the cart equal to 0.3 m in 3 s, are used in designing all the motion profiles.

Method application with unperturbed plant
The first comparison is set by focusing on the cart acceleration and jerk through their peak and RMS values. These data are summarized in Table 1 that compares seven motion profiles: four of them refer to the combination of the two shapers with the two mentioned unshaped motion profiles; the other three refer to three different tunings of the robust TPBVP approach proposed in this paper. In particular, β = 0, β = 20, β = 50, and γ = 500 have been set to reflect with null, 'small', and 'large' weighting of the peak load sway according to the formulation of the cost function in Eq. (35). The scaling factor γ is then set to 500 in all the numerical results reported in this work as well.
Jerk values are not defined when the harmonic profile is involved since such a motion primitive has a discontinuous acceleration. The analysis of the peak acceleration values shows that using the NZV shaper with the poly5 reduces the acceleration; despite the smaller accelerations of the unshaped harmonic motion primitive, the 'overlapping' acceleration profiles due to the convolution with the shaper make the peak acceleration higher in the case of harmonic profile. In contrast, the ZVD shaper results in the highest peak and RMS accelerations as a result of the need to compensate for the delays introduced by the cascaded shapers ω 1 = 3.64 rad/s and ω 2 = 9.02 rad/s, the total delays introduced by the shapers are equal to 0.808 s and 1.789 s, respectively. In the latter case the original motion profile has to be 'compressed' to a duration of a mere 1.211 s.
The peak and RMS cart accelerations required by the robust TPBVP profiles significantly change as β changes, since such a tuning parameter trades off between actuation effort (i.e., cart accelerations) and peak load sway reduction. The data in Table 1, however, show that for any of the robust TPBVP profiles the peak value of the cart acceleration is still lower than the value found for any of the robust shapers, at least as long as β does not exceed 50. The value of β defines also the measures of the cart jerk: comparing the RMS jerk, which is a common measure of the actual smoothness, shows that all three robust TPBVP profiles produce a smoother motion profile than the one generated by the ZVD shapers.
The commanded cart speed and acceleration profiles are shown in Figs. 2 and 3, respectively. Each figure shows in the upper plot the results related to input shaping and in the lower graphs the results obtained by the robust TPBVP profiles. The acceleration plot  clearly shows the discontinuous acceleration sported by the use of the harmonic with both the shapers and the jerky motion required by the poly5 together with the ZVD. The results of the application of the seven motion profiles designed are shown in Fig. 4, which shows the time evolution of the oscillation of mass m 2 around its equilibrium position (computed as L 1 sin(ϑ 1 ) + L 2 sin(ϑ 1 + ϑ 2 )) in the case of an unperturbed plant. As far as transient oscillation is concerned, the minimal peak value is found for the combination of the harmonic motion profile and the NZV shaper. A just slightly larger peak oscillation is found for the robust TPBVP profiles when β are 20 and 50, i.e., when the 'minmax' cost function (35) is adopted. In this case, however, the small peak transient load sway is obtained together with the smallest RMS transient oscillations and the increased robustness, which cannot be achieved by the NZV shaper.
As for the load response after the motion completion, the residual oscillations are almost zero for all the robust TPBVP profiles thanks to their accuracy that is often recognized as one of their main advantages. In contrast, larger oscillations are found in all the tests exploiting the shapers. This result is mainly caused by the rounding-off of the time of application of the impulses used in the shapers that is approximated with a sample time equal to 1 ms. The actual values of the load oscillations, written in terms of peak and RMS transient and residual load sway, are reported for a more detailed evaluation of the performances in Table 2.

Method application with perturbed L 1
The data in Table 2 also include the peak residual oscillation when the uncertain parameter L 1 is increased by 20%. This result provides a first indication on the robust performance with respect to the variable adopted for defining the sensitivities. While all the profiles produce rather small residual oscillation under nominal conditions, only robust design approaches (both TPBVP and ZVD) ensure small residual oscillations, whose peaks are about one third or one fourth of those obtained with the nominal methods.
The robustness performance of the investigated motion profiles can be analyzed in more detail by observing the data plotted in Fig. 5. The graph on the left shows the peak residual oscillation vs. the change of L 1 , whose nominal value is equal to 0.47 m. This graph shows the sensitivity curves for all the motion profiles and highlights that the choice of the specific cost function of Eq. (35) introduces also a trade-off between robustness enhancement and transient oscillation reduction.
The graph on the right-hand side of Fig. 5 shows instead the sensitivity curves computed for the peak value of transient oscillation, i.e., the load oscillation that happens during motion, which is an often overlooked measure in the literature on input shaping. The light blue line, which refers to the robust TPVBP profile with β = 50, is the closest one to the 'best' one, which is the one related to the harmonic/NZV combination, at least in the proximity of nominal conditions. The robust performance of the robust TPBVP profiles degrades sensibly whenever L 1 is lower than 0.4 m: this feature highlights the complexity of the relationship found between several features and targets of the motion profile design procedure, which must balance several objectives such as limited transient and residual oscillation and the moderation of the 'control' effort. The availability of the tunable weighting factor β, however, aids the achievement of the desired trade-off.

Method application with random perturbations
Finally, all the motion profiles have been tested against several sources of uncertainties to provide some information on the capability of each method of reacting to unmodeled perturbations. The outcomes of this Monte Carlo analysis are summarized in Fig. 6: the results on the bar chart on the left report the peak value of transient oscillations in groups of eight tests: each test, from left to right, refers to the perturbation, by −10% and then by +10%, to the parameters m 1 , m 2 , l 1 , and finally l 2 from their nominal values. As far as the peak transient oscillation is concerned, the best solutions are still the ones based either on the NVD shaper and the harmonic motion profiles or on the TBPVP solution with weighted transient oscillations (for β = 50). It should be noted that almost identical results have been obtained for all the perturbations.
The chart on the right-hand side of Fig. 6, in contrast, does show that the amplitude of residual oscillation can be significantly affected in a different way by different perturbations.
In particular, such a graph shows that the nonrobust motion profiles obtained through the NZV shaper are remarkably sensitive to changes of the rope lengths L 1 and L 2 , while they are less affected by alterations of similar magnitude of the two masses. As far as the robust motion profiles are concerned, the two robust TPBVP profiles behave slightly worse than the ones based on the ZVD: this occurrence might suggest that the ZVD shows the best 'overall' robustness.
All the results presented so far nonetheless showcase the effectiveness and the flexibility of the proposed method, which can combine the robustness properties of a ZVD shaper with the reduced transient oscillation amplitude achieved by a nonrobust shaper.

System description and dynamic model
The second test-case consists of designing the optimal rest-to-rest motion for a planar translating slender beam actuated by a cart to obtain null residual oscillations and transient oscillations of minimal amplitude. The dynamic model of the translating beam is developed according to the FEM-ERLS model, firstly developed and tested in the work [16]. This model uses a finite element discretization to describe the elasticity of the beam, while the concept of equivalent rigid-link system (ERLS) is used to represent the gross motion of the cart and to introduce a moving reference for defining the small elastic displacements of the beam. Such a formulation has been adopted since it directly leads to ordinary differential equations in a straightforward way. To reduce model dimension, a single Euler-Bernoulli beam element, with an additional tip mass, is adopted in this work and null longitudinal displacements are reasonably assumed. As a result, the inertial and elastic properties of the beam are represented, respectively, by a 4 × 4 mass matrix M and a 4 × 4 stiffness matrix K. The elastic displacement of the beam as a result of its elasticity is represented therefore as a set of four nodal elastic displacements: a lateral and a rotational finite displacement at the base and at the tip of the beam.
The position of the end-effector of the cart is represented by y cart (t). Four nodal elastic coordinates are initially adopted to represent the beam with free-free boundary conditions and are expressed in the local reference frame of the ERLS: The coordinates u y i are the nodal elastic deflections of the ith node and in the y direction, while the u z i are the elastic rotation in the xy plane with respect to the z axis, as shown in Fig. 7. According to the model presented in [16], the dynamics of the system can be first written as a system of five second-order differential equations: in which appears also the mass of the cart M c , the force exerted on it by the actuator F y , and the Jacobian matrix J , which is referred to as the matrix of nodal sensitivity coefficients according to the nomenclature of [16]. Such a Jacobian matrix is used to relate the speeds of the ERLS to the speeds of the generalized coordinates of the ERLS. In this case the only generalized coordinate needed to describe the motion of the ERLS is the position of the cart y cart (t). The system in Eq. (37) must then be changed to reflect the clamping conditions of the beam: in this case both nodal displacements of the beam at the first node, i.e., u y 1 and u z 1 must be set to zero and removed from Eq. (37). In all cases in which J T MJ M c , the beam dynamics does not affect significantly the cart dynamics, as in the case under consideration here, hence the generalized coordinate y cart (t) can be considered, again, as the control input to the unactuated subsystem made by the translating beam. Under this assumption, Eq. (37) can be replaced by the system of differential equations Recalling the formulation already developed in Sect. 2, it is straightforward to recognize thatq u =ü and thatq a =q a =ÿ cart are the vectors of unactuated and actuated accelerations, respectively. If, again, it is assumed that the dynamics of the closed-loop control that ensures the cart motion tracking can be approximated by the first-order model already assumed in Eq. (6), the equivalent of Eq. (8) can be written for the translating beam as follows: As for the previous test-case, the external control input is the acceleration reference signal The model in Eq. (9) requires a first-order formulation, which is obtained by adding integrators to it according to the following state vector: The model used for planning the motion of the beam is, as already stated, based on a single finite element. This choice is motivated by two main issues: the first one is related to the limitations of the numerical method used to solve the variational problems; the other one is related to the frequency range imposed by the modeled actuator bandwidth. It must be pointed out that the chosen state vector x(t) has seven entries, and each additional finite element increase its size by 4. The same increase is also reflected in the number of sensitivity functions, meaning that x r is increased by eight elements for each additional finite element. Since each element of x r requires a Lagrangian multiplier to enforce the dynamic constraints (see Eq. (26)), the size of the optimal control problem escalates quickly, with the number of finite elements being of size 22, 38, and 54 for 1, 2, and 3 finite elements, respectively.
The increase in the number of finite elements, moreover, does not necessarily lead to a better performance of the planning method. The main effect of a finer discretization of the finite element modeling is the inclusion of a larger number of vibrational modes, but such additional modes are actually relevant only if they lie within the bandwidth of the actuator, which in the case under consideration is set by the time constant τ . Modes that exceed the actuator bandwidth are minimally excited by the cart motion since frequency components beyond the bandwidth cannot be executed accurately by the actuator; so, including them in the model used for planning does not alter significantly the outcomes in real systems.

Problem definition
The planning procedure applied to the translating beam is, again, targeted at defining a restto-rest motion with null residual oscillations and reduced transient oscillations. Moreover, a robust performance is sought: the uncertain parameter is assumed to be the weight of the concentrated mass located at the free end of the beam m P , whose value is embedded in the mass matrix M used in Eqs. (38) and (9). m P being the uncertain parameter, the sensitivity functions to be included in the extended dynamics formulations are: The state vector is defined as in Eq. (43) since, again, the motion of the cart is assumed to be unaffected by the alteration of the physical properties of the beam. The augmented state vector used to set up and solve the associated TPBVP will therefore be The dynamic model of the translating beam augmented with sensitivities, i.e., to be used for the robust motion planning procedure, is written as in Eq. (34), with the only difference of having partial derivatives performed with reference to m P rather than to L 1 . As in the other test-case, the approximate 'minmax' cost function of choice weighs the actual lateral displacement of the beam tip u y 2 using this formulation: which trades off, through the weight β, between the actuation effort and the peak transient oscillations. The value γ = 500 has been used for all the computation of all the results provided in this work. As in the other test-case, here g r = g, so both the nominal and the robust motion planning problems use the same cost function.
The performance evaluation will involve the same benchmark planning procedures already defined for the previous test-case, i.e., using the same kind of motion profiles and the same types of shapers. The design specifications call for a tighter planning as the motion of the beam is set to cover a translation of the cart over a distance equal to 0.4 m in just 1 s. The beam is made of aluminum, has a rectangular cross section equal to 2 × 23 mm, is 0.915 m in length, and the nominal tip mass is m P = 10 g.
As already mentioned, the model used for planning is based on a single finite element discretization, and its natural frequencies are f n,1 = 1.53 Hz and f n,2 = 14.34 Hz. The model used for evaluating the performance of the planned motion is, however, based on two finite elements. As a result, the two models differ mainly by the fact that the second one has two additional vibrational modes, located at f 3 = 33.3 Hz and at f 4 = 87.8 Hz: both of them exceed the actuator's bandwidth set by the value τ = 0.02 s. The first two natural frequencies are, instead, located at 1.57 Hz and 14.42 Hz The model used for testing the motion profile, therefore, reflects a model-plant mismatch through both parametric uncertainties and unmodeled dynamics.

Results
The first comparison of the outcome of the benchmark and proposed planning methods is defined by the analysis of the data collected in Table 3. The data include some simple met- rics on the cart motion profile: the first column refers to the peak acceleration and shows that the application of the ZVD shaper results in two largest values among the seven shown, meaning that it implies a rather large actuation effort. The robust TPBVP profiles provide, according to the same metric, a better performance, but the absolute best behavior is found when planning with the NZV shaper. Differently from what has been observed in the other test-case, and in particular shown in Table 1, the specific choice of the motion profile just minimally alters the peak accelerations. A similar trend is observed also for the RMS acceleration (see column 3): again the variational formulation offers a good middle ground between the two solutions based on shapers. The cost to be paid for the reduction of acceleration is the noticeable increase in the peak and RMS jerk, meaning that the variational formulation defines a profile that is less smooth than the one generated by the ZVD.
The actual planned motion profiles are represented in Figs. 8 and 9. The first figure shows that the variational formulation defines a motion profile with a rather low peak speed, while  the second one highlights the pronounced accelerations observed, for the same trajectories, during the initial and final phases of the motion task. The application of the planned motion profiles to the translating beam, performed in a numerical environment, highlights the time evolution of the tip sways reported in Fig. 10.
As observed for the double pendulum test-case, the variational formulation is rather effective in confining the peak sway: the proposed robust formulation can reduce them below 5 cm, thus achieving values that are comparable with the ones based on the NZV shaper and well below the ones associated with the ZVD shaper. A more precise evaluation of this issue is provided in Table 4. Figure 10 reveals the presence of some visible residual oscillations in the case of input shaping, as corroborated by Table 4, whose frequency is both the one of the two modeled vibrational modes (in particular the first one is the most evident in Fig. 10) and of the two unmodeled ones as well. The excitation of the unmodeled dynamics is more evident in the case of input shaping as the actuator bandwidth limitation is not taken into account, and thus the planned motions can include high frequency content. This is exacerbated when the harmonic motion profile is used because of the presence of acceleration discontinuities. On the other hand, the robust TPBVP profiles, as shown in the left-hand side graph of Fig. 10 minimally excite the two unmodeled natural frequencies: in this case their presence is minimally relevant to the outcome of the simulations. The robustness sported by the methods based on the NZV shaper is so limited that a residual oscillation is observed even when the endpoint mass m P is left unaltered: as visible in Fig. 11 and in Table 4, such an amplitude is around 1.31 mm. This occurrence is due to the error introduced by the rounding-off of the time instants of application of the impulses of the shapers and by the model-plant mismatches. Figure 11 also shows the sensitivities of the residual oscillations to the change of tip mass m P : such sensitivities are rather uniform for all robust profiles as in the worst case the amplitude due to a +50% alteration of m p is kept below 2.4 mm, as listed in detail in Table 4. As far as residual oscillations are concerned for null parametric perturbation, their amplitude is reduced by one order of magnitude when switching from the nonrobust NZV shapers to the ZVD shapers and then again to the robust TPBVP profiles. The latter, in practical, has the highest insensitivity to both unmodeled dynamics and parametric uncertainties.
Finally, sensitivity curves have been evaluated and displayed in the right-hand side plot of Fig. 11 to show how the amplitude of transient oscillation changes with the change of endpoint mass m P within a ±50% range. The figure shows that the most effective way of reducing transient oscillations is either using a robust TPBVP profiles with a large value of β or using the NZV shaper and a harmonic motion profile. The second option, however, is not suitable in all cases in which enhanced robustness is required by severe uncertainties. On the other hand, the ZVD shaper provides a robust action although it does not have the capability of confining peak transient oscillation as sported by the robust TPBVP profiles. Therefore, the approach proposed in this paper seems capable of combining the best features of the two other solutions, incorporating both the robustness of the ZVD shaper and the reduced transient oscillations of the NZV shaper.

Method application with random perturbations
Finally, a general evaluation of the robustness to several sources of unmodeled uncertainties is provided by means of Fig. 12 that summarizes the results of the Monte Carlo analysis. In particular, all results are provided for changed to the beam length L by ±5% from its nominal values, while endpoint mass m P and the Young modulus E are altered by ±10% from their nominal values. As highlighted in the previous test-case, peak transient oscillations are minimally affected by any of the perturbations. On the other hand, residual oscillations are remarkably affected by the type of perturbation, as highlighted by the chart on the right-hand side of Fig. 12. Changing the mass of the end-point m P affects by a very small amount the amplitude of residual oscillations for any of the four robust motion profiles-in such cases the ZVD and the proposed robust TPBVP profiles behave almost identically. They sport also similar performances as the other parameters are changed, i.e., when either the beam length L or the Young modulus E are altered, with a general improvement brought by the robust TPBVP profiles. The robustness sported by the solutions based on the NZV is very limited-this applies to any of the alteration included in Fig. 12, as observed in Fig. 11 for the perturbation of just mass m P .

Conclusion
In this work a method for the design of motion profiles for underactuated multibody systems with enhanced parametric robustness has been presented. The method is based on the definition of a robust optimization problem constrained to the plant dynamics, which is set up with the aid of parametric sensitivity functions extracted from the dynamic model of the plant and by exploiting Pontryagin's minimum principle. The performances are assessed by comparison with the application of NZV and ZVD shapers, showing that, for the case under consideration, the proposed method has similar robustness properties of the widely adopted ZVD shaper while requiring lower actuator effort and bandwidth. Additionally, a suitable definition of the cost function adopted for the synthesis of the motion profile allows a significant reduction of the peak transient oscillation at the cost of a minor increase of the control effort without compromising the robustness properties.
Funding Note Open access funding provided by Università degli Studi di Padova within the CRUI-CARE Agreement.

Competing Interests
The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.