Simulation and Experimental Validation of Novel Trajectory Planning Strategy to Reduce Vibrations and Improve Productivity of Robotic Manipulator

: This paper aims at investigating vibrational behaviors of the industrial manipulator Racer 7-1.4, designed and manufactured by COMAU S.p.A., with the target of new trajectory planning strategies to improve productivity rate without any loss of positioning accuracy. Starting from the analysis of a 9DoF multi-body system with lumped parameter, the ﬁrst natural frequency of the robot was calculated in seven reference positions. Then, static and dynamic simulations were run by applying saturated ramp input and large motions to analyze the vibrational behavior of the manipulator. This research underlines that the optimal way to design the robot move is to set its duration at twice a period of free oscillation according to the ﬁrst vibrational mode. Due to strong analogy of dynamic response of both 1DoF and 9DoF robot models, the closed-form solution of the 1DoF undamped system—featured by natural frequency equal to the ﬁrst frequency of the 9DoF system—may be successfully adopted by the real-time trajectory planning process to predict residual vibration at move end-condition. This strategy was conﬁrmed by experimental tests, allowing either residual vibration decrease and execution time reduction as well.


Introduction
In the last years, the use of robots increased in several fields. Besides the industrial applications, robots are being used in agriculture, high risk military operations, underwater exploration, surgery and transports. In particular, robots carry out many activities in industrial field: pick and place, assembly, welding, transport, painting and much more.
The trajectory planning of a robotic arm is extremely important in production systems because it allows to carry out all the operations it is programmed to do in less time as possible, increasing the productivity and reducing the production costs [1][2][3]. Moreover, the trajectory planning can be optimized to minimize the distance traveled [4], the energy used [5], avoid obstacles [6,7].

Trajectory Generation
The main purpose of the trajectory planning task is the generation of reference inputs for the control system of the manipulator, in order to perform the motion. The inputs of any trajectory planning algorithm are the geometric path and the kinematic and dynamic constraints; the output is

Manipulators Design and Control
Apart from the specific strategy, the motion time-history generated by the trajectory planner must fulfill the constraints set a priori on the maximum values of the joint generalized torques, with no mechanical resonance excitation as fundamental target. This can be achieved by forcing the trajectory planner to generate smooth trajectories, i.e., trajectories with good continuity features: in particular, it would be desirable to design trajectories with continuous accelerations, so that the absolute value of the jerk (i.e., time-derivative of acceleration) keeps bounded. Jerk limitation is quite crucial since high jerk values may result in quick wear of the robot mechanical parts and heavily excite its resonance frequencies. Vibrations induced by non-smooth trajectories can damage robot actuators and introduce large errors while a robot is performing tasks such as trajectory tracking. Moreover, low-jerk trajectories can be executed more rapidly and accurately [8].
Most existing robotic manipulators are designed and built in a manner to maximize stiffness in an attempt to minimize the vibration of the end-effector to achieve a good position accuracy. This high stiffness is achieved by using heavy material and a bulky design. Hence, the existing heavy rigid manipulators are shown to be inefficient in terms of power consumption or speed with respect to the operating payload. Also, the operation of high precision robots is severely limited by their dynamic deflection, with strong residual oscillations even after the move end. The settling time required for this residual vibration delays subsequent operations, thus conflicting with the demand of high productivity. These conflicting requirements between high speed and high accuracy are the core of challenging research problems about robot assemblies. Also, many industrial manipulators face the problem of arm vibrations in high speed motion. In order to improve industrial productivity, it is imperative to achieve light-weight design and/or speed increase of robot operation. For these purposes, it is very desirable to build flexible robotic manipulators. Compared to the conventional heavy and bulky robots, flexible link manipulators have the potential advantage of lower cost, larger work volume, higher operational speed, greater payload-to-manipulator-weight ratio, smaller actuators, lower energy consumption, better maneuverability, better transportability and safer operation due to reduced inertia. But the greatest disadvantage of these manipulators is the vibration problem due to low stiffness. Flexible manipulators can find many applications but since the main problem is to control their vibrations, many researchers have tried to solve this problem by improving the dynamic models and incorporating different control strategies [16].
Not only the design, but also the control strategies aim to improve the accuracy and the speed of robot operations. Several techniques have been developed about trajectory tracking of robot manipulators: terminal sliding mode (TMS) control strategy for global approximate fixed-time tracking of robot manipulators [17], geometric homogeneity technique (GHT) for global finite-time tracking [18], TSM control for finite-time tracking [19]. Basically, GHT control requires the knowledge of robot dynamics, contrary to TMS control.

Aim of the Research
This research aims to identify a new strategy to minimize vibrations in an industrial robot during the execution of an assigned spatial trajectory, preserving or reducing the execution time, acting on the trajectory generation. The original contribution of this work is given by the approach adopted in the new strategy identification method: the vibrational behavior of a complex non-linear system was approximated with a hybrid solution deriving from a simplified 1DoF elastic system.
In this paper we will refer to two types of trajectory, called "non-short move" and "short move". In the first case, with "non-short move" we indicate the most general trajectory shape generated by the manipulator trajectory planner according to the following kinematic variables: where a', a" and v' are, respectively, the peak values of acceleration, deceleration and velocity defining the generated trajectory, while acc max , dec max and v max are the maximum assignable values of acceleration, deceleration and velocity determined by technological constraints of the system (motors and gearboxes). Such a trajectory is generated starting from a seven-part jerk curve, assumed with parabolic profile in part 1, 3, 5 and 7, and zero value in part 2, 4 and 6. The displacement profile is obtained by subsequent integrations of the jerk curve. Consequently, the displacement curve shows a profile depicted by fifth order polynomial. In the second case, with "short move" we introduce trajectory featuring the following kinematic variables: a' < acc max , a" > dec max , causing the lack of part 2 and 6 in the jerk curve. In particular, say, type 1 short move when: Otherwise, the type 2 short move is obtained with the lack of part 4, when: v' < v max . Figure 1 shows type 1 short moves, with t 2 and t 6 equal to 0. Starting from the jerk profile ( Figure 1a) the displacement profile ( Figure 1d) is obtained by further integrations.
This paper presents the extension of previous research, Ref. [20], with generalization of trajectory planning strategy as main objective. Such a result can be adopted both for short and non-short moves, through appropriate assumptions herein discussed.
In this analysis, the industrial manipulator Racer 7-1.4 ( Figure 2) produced by COMAU SpA was assumed as reference industrial manipulator; the mathematical modelling involves Racer 7-1.4 design data for numerical analysis and ensuing validation. This robot architecture is a 6 controlled axis industrial manipulator with 7.0 kg (69 N) maximum wrist payload, whereas maximum allowed Electronics 2020, 9, 581 4 of 17 horizontal displacement is 1436 mm (Figure 3). It is composed by 7 links forming an open kinematic chain; each link is identified by progressive number starting from Link 0 (the base) until Link 6 (the end effector). The links are connected to each other by means of joints actuated by electrical motors coupled to gearboxes. This manipulator has a wide range of applications: handling and manipulation of objects, welding, measurement and assembly purpose [21].
Electronics 2020, 9, x FOR PEER REVIEW 4 of 17 The investigation was carried out according to the four following phases: 1. multibody modeling of the industrial manipulator Racer 7-1.4 in Simscape Multibody environment; 2. identification of a quantitative index for the online calculation of the end effector oscillation, with the aim of optimizing the input signals generation process in case of "short move"; 3. experimental validation of analytic and numerical results obtained for short moves simulations; 4. extension of point 2 to "non-short move" case.   The investigation was carried out according to the four following phases: 1. multibody modeling of the industrial manipulator Racer 7-1.4 in Simscape Multibody environment; 2. identification of a quantitative index for the online calculation of the end effector oscillation, with the aim of optimizing the input signals generation process in case of "short move"; 3. experimental validation of analytic and numerical results obtained for short moves simulations; 4. extension of point 2 to "non-short move" case.

Definition and Validation of the Elastic Model in Simscape Multibody Environment
As first step, an elastic lumped parameter model of the manipulator was realized, in order to run numeric simulations of the vibrational behavior. Then, the first natural oscillation frequency was calculated in seven known configurations.
The seven natural frequencies were compared to those evaluated by the trajectory generating software now adopted by COMAU on the manipulator, in order to validate the virtual model.
The robot model was developed in Matlab Simulink Simscape Multibody environment, starting from the CAD 3D model of the manipulator, imported in Simulink environment associated to 2017b Matlab release. The graphic reproduction tool during virtual simulation, though redundant by the elaborative point of view, allows to display the spatial evolution.
The model is made assuming rigid links connected to each other by spherical joints that simulate the lumped elastic behavior. In particular, each spherical joint presents three lumped elastic stiffnesses (related to the three orthogonal directions) obtained by FEM analysis and no dumping, while the rigid links present values of mass and inertia obtained by CAD analysis on the real robot.
To calculate the first natural frequency of the virtual manipulator in the seven reference positions, the "Spectrum Analyzer" tool of Simulink was adopted. It accepts the displacement signals measured at the end effector as input and the output is the frequency spectrum of the three input signals related to the three orthogonal axes.
An exciting force was simulated with application point at the end effector. The force signal is a ramp with maximum value of 1000 N and saturation time of 1 s. This kind of excitation was independently applied along three orthogonal directions.
The results are summarized in Table 1 Table 2 Table 3, where they are compared with the values calculated with a model implemented in another software environment previously adopted by COMAU. It is possible to find 3.17% maximum difference about the first vibration mode frequency with position 7, demonstrating the validity of the proposed model. The investigation was carried out according to the four following phases: 1.
identification of a quantitative index for the online calculation of the end effector oscillation, with the aim of optimizing the input signals generation process in case of "short move"; 3.
experimental validation of analytic and numerical results obtained for short moves simulations; 4. extension of point 2 to "non-short move" case.

Definition and Validation of the Elastic Model in Simscape Multibody Environment
As first step, an elastic lumped parameter model of the manipulator was realized, in order to run numeric simulations of the vibrational behavior. Then, the first natural oscillation frequency was calculated in seven known configurations.
The seven natural frequencies were compared to those evaluated by the trajectory generating software now adopted by COMAU on the manipulator, in order to validate the virtual model.
The robot model was developed in Matlab Simulink Simscape Multibody environment, starting from the CAD 3D model of the manipulator, imported in Simulink environment associated to 2017b Matlab release. The graphic reproduction tool during virtual simulation, though redundant by the elaborative point of view, allows to display the spatial evolution.
The model is made assuming rigid links connected to each other by spherical joints that simulate the lumped elastic behavior. In particular, each spherical joint presents three lumped elastic stiffnesses (related to the three orthogonal directions) obtained by FEM analysis and no dumping, while the rigid links present values of mass and inertia obtained by CAD analysis on the real robot.
To calculate the first natural frequency of the virtual manipulator in the seven reference positions, the "Spectrum Analyzer" tool of Simulink was adopted. It accepts the displacement signals measured at the end effector as input and the output is the frequency spectrum of the three input signals related to the three orthogonal axes.
An exciting force was simulated with application point at the end effector. The force signal is a ramp with maximum value of 1000 N and saturation time of 1 s. This kind of excitation was independently applied along three orthogonal directions. The results are summarized in Table 1 Table 2 Table 3, where they are compared with the values calculated with a model implemented in another software environment previously adopted by COMAU. It is possible to find 3.17% maximum difference about the first vibration mode frequency with position 7, demonstrating the validity of the proposed model. Table 1. Comparison of the first natural frequency in the seven reference positions obtained by applying a force with saturated ramp profile at the end effector along x-axes (no damping).

Position
First  Table 3. Comparison of the first natural frequency in the seven reference positions obtained by applying a force with saturated ramp profile at the end effector along z-axes (no damping).

Simulations Analysis and Results
Once defined and validated the 3D elastic manipulator model, it was employed to run static simulations (to study the vibrational behavior in case of excitation that cause local oscillations around the equilibrium position) and dynamic simulations with "large" motions in different operative conditions, in order to identify a smart way to predict the oscillation at move planning stage.

Analytic Quantification of the End Effector Oscillation Entity for the Optimization of "Short Move" Calculation
In a first stage, the analysis was carried out on simple 1-degree-of-freedom (or 1DoF) and 2-degrees-of-freedom (or 2DoF) models, to arrive, later, to a complete 9DoF (related to Joint 1, 2 and 3, respectively the spherical joints that connect links 0-1, 1-2 and 2-3) model to evaluate any possible analogy with the simpler systems. As will be shown, the 9DoF system presents a strong similarity to the 1DoF system. For this reason, the analytic solution of the ODE describing 1DoF system motion was later adopted to calculate the synthetic index that quantifies the end effector oscillation.

Modification of the Simscape Multibody Model for Motion Actuation
In order to simulate "large" motions in Simscape Multibody software environment and, mainly, short and non-short moves, a model upgrade was implemented. A new element called "bushing" was incorporated in the model as part of Joint 1, 2 and 3. The bushing is made of a negligible mass disc in revolute pair with the rigid link; the addition of the negligible mass disc doesn't affect the inertial properties of the model. Joints 1, 2 and 3 now present three elements in series connection each: a revolute joint, a massless disc and a spherical joint. In this configuration, the bushings simulate the electric motors, allowing the free rotation with respect to the controlled axis, while the spherical joints keep containing the elastic properties by way of the finite stiffnesses. The joints are still assumed with no damping and the gravity effect was disabled to automatically compensate the oscillation component due to the weight force of the several links acting on the manipulator itself.

End Effector Vibration Entity Quantification During Simulation
During simulations, it was necessary to introduce a comparative index allowing to quantify and compare the entity of the end effector oscillations according to the excitation. This index, also called "cost function", was defined as follow: where x(t), y(t) and z(t) represent the measured position of the end effector along the three directions and x ob (t), y ob (t) and z ob (t) the target position the end effector should nominally reach. The cost function accumulates the position error along the three directions during the whole simulation.
In the same way, three cost functions were defined, each one related to a different direction: Because of the lack of damping, the vibration persists during time and, considering an execution time reasonably longer than oscillation period, it is possible to neglect the contributions of oscillation period fractions at the integration range edges.

Displacement Response Spectrum of 1DoF and 2DoF Models
Next to the usual representation in time domain, it is very effective to display in graphic form the locus of maximum amplitudes of oscillation with respect to different input signals characteristic timings and natural oscillation period ratios, through the Displacement Response Spectrum (DRS). In this study, the cost function int was represented as function of the characteristic time/natural oscillation period ratio. The characteristic time of an arbitrary motion or force input is usually related to the input signal duration or its period. For a saturated ramp input the characteristic time corresponds to the time range needed to reach the saturation value. For a semi-sinusoidal input, the characteristic time coincides with the completion of the semi-sinewave.
If we define an adimensional time index τ as: τ = t 1 /T n where t 1 is the characteristic time and T n is the natural oscillation period, we know from the theory that the DRS of a 1DoF system excited by a saturated ramp input (both in force or position) with characteristic time t 1 shows a typical trend with minima for integer values of the adimensional time [22]. As expected, the Racer manipulator multibody model has shown the same dynamic behavior after making it similar to a simpler 1DoF elastic system, with all the degrees of freedom constrained except for vertical rotation at Joint 1. It is worth nothing that at this stage the bushings are clamped to carry out a static analysis. First, the (only) natural frequency of such a system was calculated in position 1, or "calibration" position: F n = 36.121 Hz. Then, the virtual manipulator in the same position was excited applying a saturated ramp input force at the end effector along local z direction. A total of 12 simulations were run with different saturation time and same maximum force value after saturation, equal to 5000 N. During simulation, the comparative indexes int, int x , int y and int z were calculated in the integration range (10t 1 ; T sim ) where t 1 is the saturation time and T sim is the period of the simulation, fixed at 30 s. The cost function int was represented in graphic form as function of the adimensional time τ, where T n is obtained as 1/F n (Figure 4). The vibrations are null when the adimensional time τ assumes an integer value, i.e., the saturation time t 1 is an integer multiple of the natural oscillation period T n . According to the theory, if we excite in position a dynamic 1DoF system we obtain, as well, null vibrations for integer τ.  Similarly, the multibody model was configured to have two degrees of freedom allowing the rotation around the local z-axes at Joint 1 and 2 (i.e., the controlled axes of the manipulator at Joint 1 and 2). In this configuration, the first natural frequency is Fn1 = 17.898 Hz. Twelve simulations were carried out in the same operative way as shown for the 1DoF model, with a variable saturation time t1 and calculating the cost functions. In Figure 5, the cost function int as function of τ is shown. The adimensional time is referred to the natural oscillation period related to the first mode of vibration. For this reason, the local minima of the cost function int for a 2DoF system occur for τ = k/2 with k ∈ N. When k is an even number, t1 is an integer multiple of the natural oscillation period related to the first mode of vibration and, even if the vibrations are not null as they are for the 1DoF system, we obtain local minima of the cost function. The reason why we have additional local minima with respect to the 1DoF system (when k is an odd number) is that the vibration level is affected by the second mode of vibration, with a natural frequency equal (in this case) to twice the first one.  Similarly, the multibody model was configured to have two degrees of freedom allowing the rotation around the local z-axes at Joint 1 and 2 (i.e., the controlled axes of the manipulator at Joint 1 and 2). In this configuration, the first natural frequency is F n1 = 17.898 Hz. Twelve simulations were carried out in the same operative way as shown for the 1DoF model, with a variable saturation time t 1 and calculating the cost functions. In Figure 5, the cost function int as function of τ is shown. The adimensional time is referred to the natural oscillation period related to the first mode of vibration. For this reason, the local minima of the cost function int for a 2DoF system occur for τ = k/2 with k ∈ N. When k is an even number, t 1 is an integer multiple of the natural oscillation period related to the first mode of vibration and, even if the vibrations are not null as they are for the 1DoF system, we obtain local minima of the cost function. The reason why we have additional local minima with respect to the 1DoF system (when k is an odd number) is that the vibration level is affected by the second mode of vibration, with a natural frequency equal (in this case) to twice the first one.

Displacement Response Spectrum of 9DoF Model
Afterwards, the analysis moved to a Racer manipulator multibody model with nine degrees of freedom and the bushings still clamped: an iterative simulation process was performed in the seven reference positions with variable saturation time and maximum value of the force applied at the end effector, in the hypothesis of no damping and no gravity. The adimensional time τ was calculated as τ = t 1 /T n1 , where T n1 is the natural oscillation period related to the first mode of vibration, equal to T n1 = 1/F n1 . It is worth highlighting that the first natural frequency changes in each position of the manipulator.
The investigation was led assuming 30 values of τ equally spaced in the range [0+; 3] and 5 values of the maximum force F max equally spaced in the range [1000; 5000] N. In total, 150 simulations for each of the seven positions were run. In every simulation, the execution time was fixed at T sim = 30 s and the cost functions were calculated in the integration range [10t 1 ; T sim ]. The adimensional index int was displayed as function of τ and F max on tridimensional graphs. In Figure 6, the results related to the calibration position are shown. It appears evident a net similarity with the DRS of the 1DoF system: for integer values of τ we obtain the minimum values of the cost function int. Moreover, it is clear the proportionality to the applied force intensity F max .
Similarly, the multibody model was configured to have two degrees of freedom allowing the rotation around the local z-axes at Joint 1 and 2 (i.e., the controlled axes of the manipulator at Joint 1 and 2). In this configuration, the first natural frequency is Fn1 = 17.898 Hz. Twelve simulations were carried out in the same operative way as shown for the 1DoF model, with a variable saturation time t1 and calculating the cost functions. In Figure 5, the cost function int as function of τ is shown. The adimensional time is referred to the natural oscillation period related to the first mode of vibration. For this reason, the local minima of the cost function int for a 2DoF system occur for τ = k/2 with k ∈ N. When k is an even number, t1 is an integer multiple of the natural oscillation period related to the first mode of vibration and, even if the vibrations are not null as they are for the 1DoF system, we obtain local minima of the cost function. The reason why we have additional local minima with respect to the 1DoF system (when k is an odd number) is that the vibration level is affected by the second mode of vibration, with a natural frequency equal (in this case) to twice the first one.   Afterwards, the analysis moved to a Racer manipulator multibody model with nine degrees of freedom and the bushings still clamped: an iterative simulation process was performed in the seven reference positions with variable saturation time and maximum value of the force applied at the end effector, in the hypothesis of no damping and no gravity. The adimensional time τ was calculated as τ = t1/Tn1, where Tn1 is the natural oscillation period related to the first mode of vibration, equal to Tn1 = 1/Fn1. It is worth highlighting that the first natural frequency changes in each position of the manipulator.
The investigation was led assuming 30 values of τ equally spaced in the range [0+; 3] and 5 values of the maximum force Fmax equally spaced in the range [1000; 5000] N. In total, 150 simulations for each of the seven positions were run. In every simulation, the execution time was fixed at Tsim = 30 s and the cost functions were calculated in the integration range [10t1; Tsim]. The adimensional index int was displayed as function of τ and Fmax on tridimensional graphs. In Figure 6, the results related to the calibration position are shown. It appears evident a net similarity with the DRS of the 1DoF system: for integer values of τ we obtain the minimum values of the cost function int. Moreover, it is clear the proportionality to the applied force intensity Fmax.

Estimation of Vibration Level at Nove End for a 1DoF System
After the results of the static analysis on the manipulator model, the investigation was extended towards the dynamic analysis releasing the bushings and actuating "large" motions.
We started from a Single input/Single output system described by a second-order ODE. It consists of a 1DoF elastic system with natural frequency equal to the one related to the first mode of vibration of the 9DoF system, to which a rotation u(t) having time functionality equal to the kinematic of the actuated joint is applied as input. According to the notation assumed in this paper, u(t) is the input and x(t) is the output, where these variables represent rotations expressed in rad. It is clear that in the ideal domain of rigid rotation or in quasi-static running, the expected result is x(t) = u(t).
A possible waveform of the positive and negative acceleration phases is obtained from portions of the function sin 2 (ω t). The angular velocity and displacement profiles of the actuated joint are obtained by further integrations of the acceleration.
The following second-order ODE, despite its compactness and simplicity: • describes the oscillatory dynamic of a 1DoF system; Figure 6. Cost function int with respect to τ and F max in calibration position for the 9DoF system.

Estimation of Vibration Level at Nove End for a 1DoF System
After the results of the static analysis on the manipulator model, the investigation was extended towards the dynamic analysis releasing the bushings and actuating "large" motions.
We started from a Single input/Single output system described by a second-order ODE. It consists of a 1DoF elastic system with natural frequency equal to the one related to the first mode of vibration of the 9DoF system, to which a rotation u(t) having time functionality equal to the kinematic of the actuated joint is applied as input. According to the notation assumed in this paper, u(t) is the input and x(t) is the output, where these variables represent rotations expressed in rad. It is clear that in the ideal domain of rigid rotation or in quasi-static running, the expected result is x(t) = u(t).
A possible waveform of the positive and negative acceleration phases is obtained from portions of the function sin 2 (ω t). The angular velocity and displacement profiles of the actuated joint are obtained by further integrations of the acceleration.
The following second-order ODE, despite its compactness and simplicity: • describes the oscillatory dynamic of a 1DoF system; • represents the operative tool in the simulation of the end effector vibrational behavior in 4, 5, 6 and 7 phases moves. ..
As it stands, the differential equation only requires the imposed motion u(t) as input and the modal parameter f, i.e., the natural oscillation frequency, due to an algebraic manipulation of the equation that allows to combine the torsional stiffness and the moment of inertia. In this way, a separate evaluation of stiffness and inertia is not required. Referring to Figure 1, the following hypothesis were adopted: f: frequency of the 1 st mode of vibration of the 9DoF system; • |a'|=|a"|, i.e., the absolute value of the maximum acceleration is equal to that of the maximum deceleration.
The solution, assuming t = 0 at the end of the move is: x fin (t) = a 8f 2 π 2 (−1+4f 2 t 2 1 ) (−16f 2 π 2 t 2 1 +64f 4 π 2 t 4 1 −8f 2 π 2 t 1 t 4 +32f 4 π 2 t 3 1 t 4 − cos(2fπt)+ cos(2fπ(t + 2t 1 ))+ + cos(2fπ(t + 2t 1 +t 4 ))− cos(2fπ(t + 4t 1 +t 4 )). (3) Equation (3) can be rearranged in the following formulation: A closed-form solution is now available. It is made of a combination of harmonics, which means it reflects an oscillatory behavior with a maximum amplitude of oscillation. The maximum amplitude or overshoot, can be chosen as quantitative index to quantify the vibration amount at move end. For the calculation of this index, for example, it is possible to sample the function x fin (t) in 30 ÷ 40 points of the period T = 1/f and extract the maximum value from this dataset. However, we can also calculate the maximum amplitude in closed-form way. The simulations show reduction of vibrations at move end at t 1 = k (1/2f) with k ∈ N (type 2 short moves with duration k (2/f) = k (2 t n )), except for k = 1, when Equation (3) presents a singularity. This case represents a short limitation in the adoption of this formulation for trajectory planning optimization. In type 1 short moves, together with the previous result, we have that for t 4 = k (1/f) with k ∈ N the vibrations are reduced. This calculation routine returns as output this kind of formulation: x max = F (a', t 1 , t 4 , f).
Still referring to a 1DoF system, the maximum oscillation amplitude x max was evaluated according to the previously identified routine as function of jerk time t j fixing the maximum acceleration value that can be reach and the natural frequency (depending on the manipulator in a given position). The results are shown in Figure 7: It is worthwhile to point out that there is always a condition where the end effector oscillation is null (x max = 0); the dashed black curve represents the delimitation line between the region where the given displacement ∆θ and the maximum acceleration can be both reach (right) and the region where it's necessary to relax the acceleration constraint to obtain the given ∆θ (left).
Electronics 2020, 9, x FOR PEER REVIEW 11 of 17 can be both reach (right) and the region where it's necessary to relax the acceleration constraint to obtain the given ∆θ (left).

Short Moves Simulation for 9DoF System
As observed in the previous section, the 1DoF elastic system dynamic response provides an analytic closed-form formulation to calculate the maximum oscillation amplitude at move end.
At this point, the investigation has moved to the analysis of the vibrational behavior of a 9DoF system, simulating type 2 short moves with the aim of detecting analogies with the 1DoF system. In this case, we decided for a numerical approach because finding a closed-form solution that describes the dynamic behavior of a complex 9DoF system is not trivial and it is not always possible. We have assumed as KPI the cost function int1s, that represents the cost function previously described by Equation (1) evaluated in the second immediately after move ending.
An assigned motion profile at bushing 1 was assumed for the first set of simulations. The rotation profile of bushing 1 was obtained integrating twice in time domain the antisymmetric acceleration profile built as a function sin 2 (π/(2t1) t) with amplitude equal to accmax. Simulations with total rotation Δθ equal to 5°, 10°, 15° and 20° and with several values of execution time were run. Hence, we obtain the value assumed by the KPI in different operative conditions, graphically represented as function of the duration tfin = 4t1 and the maximum acceleration accmax (Figure 8). In this way, we can proceed to a rapid identification of the favorable conditions from the vibrational point of view. It is worth highlighting that, fixing the move duration and the final position, the acceleration peak value is constrained. Thus, the KPI is showed on two different graphs: on one hand we have the vibrational index displayed as function of the move duration, on the other hand the same index is shown as function of the acceleration peak value. In the first case, the graph shows a vertical straight line, too, specified as 4tn. It represents the reference temporal condition for the planning method usually adopted by COMAU with tjerk = tn. The value of tn is given by the natural oscillation period corresponding to the first mode of vibration in each of the 7 given starting position (Tables 1-3). Independently of the entity of the rotation assigned to bushing 1, we find a local minimum of the index int1s for a move duration tfin ≈ 2tn, except for the starting position 7, where minimum is not visible. Consequently, we can obtain a local minimum of int1s for a maximum acceleration value accmax depending on the assigned rotation.
Δθ= ∫∫ acc(t) dt = f (accmax, t1) Figure 7. Representation of x max as function of t j for 1DoF system, fixing maximum acceleration and natural frequency.

Short Moves Simulation for 9DoF System
As observed in the previous section, the 1DoF elastic system dynamic response provides an analytic closed-form formulation to calculate the maximum oscillation amplitude at move end.
At this point, the investigation has moved to the analysis of the vibrational behavior of a 9DoF system, simulating type 2 short moves with the aim of detecting analogies with the 1DoF system. In this case, we decided for a numerical approach because finding a closed-form solution that describes the dynamic behavior of a complex 9DoF system is not trivial and it is not always possible. We have assumed as KPI the cost function int 1s , that represents the cost function previously described by Equation (1) evaluated in the second immediately after move ending.
An assigned motion profile at bushing 1 was assumed for the first set of simulations. The rotation profile of bushing 1 was obtained integrating twice in time domain the antisymmetric acceleration profile built as a function sin 2 (π/(2t 1 ) t) with amplitude equal to acc max . Simulations with total rotation ∆θ equal to 5 • , 10 • , 15 • and 20 • and with several values of execution time were run. Hence, we obtain the value assumed by the KPI in different operative conditions, graphically represented as function of the duration t fin = 4t 1 and the maximum acceleration acc max (Figure 8). In this way, we can proceed to a rapid identification of the favorable conditions from the vibrational point of view. It is worth highlighting that, fixing the move duration and the final position, the acceleration peak value is constrained. Thus, the KPI is showed on two different graphs: on one hand we have the vibrational index displayed as function of the move duration, on the other hand the same index is shown as function of the acceleration peak value. In the first case, the graph shows a vertical straight line, too, specified as 4t n . It represents the reference temporal condition for the planning method usually adopted by COMAU with t jerk = t n . The value of t n is given by the natural oscillation period corresponding to the first mode of vibration in each of the 7 given starting position (Tables 1-3). Independently of the entity of the rotation assigned to bushing 1, we find a local minimum of the index int 1s for a move duration t fin ≈ 2t n , except for the starting position 7, where minimum is not visible. Consequently, we can obtain a local minimum of int 1s for a maximum acceleration value acc max depending on the assigned rotation. The same set of simulations was carried out exciting bushing 2, realizing a rotation of the manipulator in the vertical plane. We have obtained, again, the representation of the cost function int1s in two different graphs, as function of the move duration tfin and the acceleration peak value accmax ( Figure 9). We can note a strong analogy with the simulations graphically described in Figure 8 (short moves in the horizontal plane). The best condition from the vibrational point of view occurs for tfin ≈ The same set of simulations was carried out exciting bushing 2, realizing a rotation of the manipulator in the vertical plane. We have obtained, again, the representation of the cost function int 1s in two different graphs, as function of the move duration t fin and the acceleration peak value acc max (Figure 9). The same set of simulations was carried out exciting bushing 2, realizing a rotation of the manipulator in the vertical plane. We have obtained, again, the representation of the cost function int1s in two different graphs, as function of the move duration tfin and the acceleration peak value accmax ( Figure 9). We can note a strong analogy with the simulations graphically described in Figure 8 (short moves in the horizontal plane). The best condition from the vibrational point of view occurs for tfin ≈ Figure 9. Index int 1s as function of t fin , acc max and rotation amplitude (rot) actuating bushing 2 in calibration position.
We can note a strong analogy with the simulations graphically described in Figure 8 (short moves in the horizontal plane). The best condition from the vibrational point of view occurs for t fin ≈ Electronics 2020, 9, 581 13 of 17 2t n , where t n is the oscillation period related to the first mode of vibration of each starting position. Repeating the set of simulations applying the motion at bushing 3, the same result is achieved. It is therefore detected an optimizing condition for vibrations at move end, unrelated to the energy application point.
At this point, it is clear that the frequency related to the first mode of vibration evaluated in the starting position of the manipulator strongly affects the oscillation entity at move end, no matter the actuated joint. In addition, it is worth noting a similarity between the DRS of the 1DoF system (Paragraph 3.1.3) and the trend of the cost function int 1s for a 9DoF system. In summary, type 2 short moves applied to 9DoF elastic system exhibit an optimal duration by the vibrational point of view equal to twice the natural oscillation period t n , where t n = f n −1 and f n is the natural frequency related to the first mode of vibration in each starting position.

Experimental Results
The outcomes of this study were implemented in the trajectory planner of a real robot to verify its deployment and effectiveness.
In Figure 10 we can see that the new planning strategy (MOD stands for modified) allows to realize a move with a vibration reduction at the end effector along the move direction and a growth of the oscillation in the two orthogonal directions with respect to the original planning strategy (STD stands for standard). However, the oscillation decrease along the main direction is one order of magnitude higher than the growth in the other two directions, resulting in a globally positive behavior due to the novel planning strategy applied to a short move. It is worth highlighting that we have achieved two goals: in fact, in addition to the vibration decrease, the new planner reduces the move duration for the same mission.
Electronics 2020, 9, x FOR PEER REVIEW 13 of 17 2tn, where tn is the oscillation period related to the first mode of vibration of each starting position.
Repeating the set of simulations applying the motion at bushing 3, the same result is achieved. It is therefore detected an optimizing condition for vibrations at move end, unrelated to the energy application point. At this point, it is clear that the frequency related to the first mode of vibration evaluated in the starting position of the manipulator strongly affects the oscillation entity at move end, no matter the actuated joint. In addition, it is worth noting a similarity between the DRS of the 1DoF system (Paragraph 3.1.3) and the trend of the cost function int1s for a 9DoF system. In summary, type 2 short moves applied to 9DoF elastic system exhibit an optimal duration by the vibrational point of view equal to twice the natural oscillation period tn, where tn = fn −1 and fn is the natural frequency related to the first mode of vibration in each starting position.

Experimental Results
The outcomes of this study were implemented in the trajectory planner of a real robot to verify its deployment and effectiveness.
In Figure 10 we can see that the new planning strategy (MOD stands for modified) allows to realize a move with a vibration reduction at the end effector along the move direction and a growth of the oscillation in the two orthogonal directions with respect to the original planning strategy (STD stands for standard). However, the oscillation decrease along the main direction is one order of magnitude higher than the growth in the other two directions, resulting in a globally positive behavior due to the novel planning strategy applied to a short move. It is worth highlighting that we have achieved two goals: in fact, in addition to the vibration decrease, the new planner reduces the move duration for the same mission.

Generalization of the Results for Non-Short Moves
A non-short move is characterized, as previously seen, by an assigned actuated joint rotation shape based on fifth order polynomials. Seven distinct phases con be detected. In general, such a move does not have any kind of symmetry. For this reason, suggesting a similar procedure to identify an analytic formulation of the overshoot to estimate vibrations, the time ranges related to the seven phases are completely independent from each other, except for t 3 and t 7 , respectively equal to t 1 and t 5 . The input generation starts from parabolic jerk profiles. Generally, the arbitrariness of the absolute values of the maximum jerk in acceleration (J max1 ) and deceleration (J max2 ) phase and that of the duration of the several phases do not assure the zeroing of the velocity at move end. To have this condition verified, it must be: In addition, in the proposed conditions we have: The duration of a non-short move can be relatively long and the position of the manipulator can evolve from an initial configuration to a final one with spatial mass distribution significantly different. This means a variation of the first natural frequency related to the final position achieved. For this reason, we chose to find a closed-form solution of an undamped elastic 1DoF system with natural frequency equal to f a until the end of part 4 of the move and then equal to f b , where f a and f b are respectively the first natural frequency in the initial and final posing. The solution of an undamped elastic 1DoF system excited by a joint actuated according to a non-short move profile with a natural frequency equal to f a until the end of part 4 and, after, f b is shown in Appendix A.
x max = F(J max1 , J max2 , t 1 , t 2 , t 4 , t 5 , t 6 , f a , f b ) As well as the previous case, a quantitative index to evaluate the vibrations entity is the overshoot, that can be calculated with a sampling of the function x fin (t) in 20 ÷ 50 points on a period equal to T = 1/f with f the lower value between f a and f b and extracting the maximum value from the array thus obtained. Compared to the previous one, this solution is numerically more stable because the denominator cannot be zero and t 1 can be equal to 1/(2f). Furthermore, introducing the appropriate symmetries, we can be led back to the previous case (Equation (3)) with a negligible deviation, confirming that the two formulations with third order polynomial acceleration of the actuated joint and sin 2 (ωt) shaped acceleration do not exhibit appreciable differences. The described procedure leads to the formulation: In particular, adopting the simplifying hypothesis of f b = f a = f and antisymmetric type 1 short move, i.e., t 2 = t 6 = 0 and t 1 = t 3 = t 5 = t 7 , we obtain no vibrations when: • t 1 = k 1/(2f), with k odd number, and t 4 = h (1/f) with h ∈ N; • t 1 = k 1/(2f), with k even number, and any value of t 4 ; • 2t 1 + t 4 = k (1/f) with k ∈ N, in general.

Conclusions
The Racer 7-1.4 robot multibody model, validated through the comparison of the first natural frequency in seven reference positions calculated through two different models of the same manipulator in two different software environments, provided an effective simulation tool to carry out virtual tests and numerical analysis. Specifically, the model allowed to explore the influence of characteristic parameters of short and non-short moves on the end effector vibrational behavior.
For the symmetric type 2 short moves, the 9DoF model shows vibrational behavior close to the simpler 1DoF system, with local minimum of the cost function occurring for a move duration close to 2t n , where t n is the natural oscillation period about the first vibrational mode. According to this similarity, the closed-form solution of a 1DoF undamped elastic system with natural frequency equal to the first mode vibration frequency of the 9DoF system can provide the equation to calculate an overall index able to quantify the oscillation amount at end effector position. In particular, the solution allows calculating the overshoot as key performance indicator of the planned trajectory to compare the vibration levels expected at move planning stage and to choose the best combination of characteristic parameters.
The adoption of the new strategy for the trajectory planning of a real manipulator validated the strong analogy between 1DoF and 9DoF systems, resulting in an overall reduction in end effector oscillations together with a significant move execution time decrease.
The study was extended, then, to the most general case of non-short move, with an input obtained by parabolic jerk waves without symmetry, considering the variation of the natural frequency of the system during the move due to the posing changing.
It is possible to lead back to specific cases and to short moves through adequate constraints on time durations t i and maximum jerk values and zeroing the segments with constant acceleration. The overshoot calculation at planning stage becomes the quantitative tool of estimation of the vibrations at the end effector, allowing the prediction of the oscillatory behavior of a complex system by the way of a closed-form equation useful to the online trajectories optimization. In this way, the trajectory planner of the manipulator is able to generate faster and higher-precision moves, enhancing both the productivity and the task quality.