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

Abstract

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 first 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 first 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 first 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 confirmed by experimental tests, allowing either residual vibration decrease and execution time reduction as well.

1. 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].

1.1. 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 the trajectory of the joints (or of the end effector), expressed as time-history of position, velocity and acceleration values. Usually, the geometric path is specified in the operating space, i.e., with reference to the end effector of the robot, because both the tasks and the physical obstacles to be skipped can be naturally described in such a space. However, the trajectory is normally planned in the joint space of the robot, after a mathematical inversion of its kinematics. Then, joint trajectories are obtained by way of interpolating functions which meet the kinematic and dynamic constraints. The main advantage of trajectory planning in the joint space comes from the control system whose action is naturally linked with joints actuators; thus, it would be easier to adjust the trajectory according to the design requirements once the planning is provided in the joint space. Moreover, trajectory planning in the joint space would allow to avoid the drawbacks arising with kinematic singularities and manipulator redundancy. The main disadvantage of trajectory planning in the joint space consists of the not easy prediction of robot end effector position, due to the non-linearities introduced by mathematical calculation from joints trajectories to of end effector trajectories through direct kinematics [8].

In the joint space, trajectories are mostly depicted by polynomial, spline, Bezier interpolation functions, etc. e.g., Cook and Ho [9] developed an algorithm based on third and fifth order polynomial functions, while Petrinec and Kovacic [10], as like as Boscariol et al. [11], combine fourth and fifth order ones. Xu et al. [12] use Bezier curve to obtain the trajectory. Further algorithms are based on third [13], fifth [8,14] and seventh order [15] B-splines.

1.2. 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.

1.3. 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:

a’ = acc_max,

a″ = dec_max,

v’ < v_max,

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:

v’ = v_max.

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₂ and t₆ 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 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].

The investigation was carried out according to the four following phases:

• multibody modeling of the industrial manipulator Racer 7-1.4 in Simscape Multibody environment;
• 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”;
• experimental validation of analytic and numerical results obtained for short moves simulations;
• extension of point 2 to “non-short move” case.

2. 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.

3. 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. 3.1. 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.

3.1.1. 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.

3.1.2. 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:

int= _∫0tfin^{(x(t)−_xob(t))2}+^{(y(t)−_yob(t))2}+^{(z(t)−_zob(t))2} dt ,

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:

_intx=_∫0tfin|x(t)_−xob(t)| dt ; _inty=_∫0tfin|y(t)_−yob(t)| dt ; _intz=_∫0tfin|z(t)_−zob(t)| dt .

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.

3.1.3. 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₁/T_n where t₁ 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₁ 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₁; T_sim) where t₁ 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₁ 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 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₁ 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₁ 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.

3.1.4. 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₁/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₁; 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.

3.1.5. 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² (ω 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.

x¨(t)^(2πf)2+x(t)=u(t).

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:

• t₁ = t₃ = t₅ = t₇;
• t₄ ≥ 0;

• v′=_∫02t1acc(t)dt;

• 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.

For the above-mentioned hypothesis, we have:

v’ = a’t₁;

∆θ =2a’t₁² + a’t₁t₄.

The solution, assuming t = 0 at the end of the move is:

_xfin(t)=a′^{8f2 π2}(^−1+4f2 _t12)^{(−16f2 π2} _t12 ^{+64f4 π2} _t14 ^{−8f2 π2} _{t1 t4} ^{+32f4 π2} _{t13 t4}−cos(2fπt)+cos(2fπ(_t+2t1))++cos(2fπ(_{t+2t1 +t4}))−cos(2fπ(_{t+4t1 +t4})).

Equation (3) can be rearranged in the following formulation:

_xfin(t)=1E^a′[A+B+C+D+αcos(2πft+β)]=F+Gcos(2πft+β)=F(1+GFcos(2πft+β)).

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₁ = 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₄ = 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₁, t₄, 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).

3.1.6. 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²(π/(2t₁) 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₁ 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 (Table 1, Table 2 and Table 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.

Δθ= ∫∫ acc(t) dt = f (acc_max, t₁)

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).

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 ≈ 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⁻¹ and f_n is the natural frequency related to the first mode of vibration in each starting position. 3.2. 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.

3.3. 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₃ and t₇, respectively equal to t₁ and t₅. 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:

_Jmax2=_{Jmax1 t12 +Jmax1 t1 t2t5}(_{t5 +t6}).

In addition, in the proposed conditions we have:

^a′=23 _{Jmax1 t1}; |^a″|=23|_Jmax2|_t5; ^v′=_{2 Jmax1 t12}3+_{2 Jmax1 t1 t2}3; Δθ=_{2 Jmax1 t13}3_{+Jmax1 t12 t2}+13 _{Jmax1 t1 t22}+(_{2 Jmax1 t12}3+_{2 Jmax1 t1 t2}3) _t4+43 _{Jmax1 t12 t5}+43 _{Jmax1 t1 t2 t5}+−_{2 Jmax2 t53}3+23 _{Jmax1 t12 t6}+23 _{Jmax1 t1 t2 t6 −Jmax2 t52 t6}−13 _{Jmax2 t5 t62}.

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₁, t₂, t₄, t₅, t₆, 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₁ 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²(ω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₂ = t₆ = 0 and t₁ = t₃ = t₅ = t₇, we obtain no vibrations when:

• t₁ = k 1/(2f), with k odd number, and t₄ = h (1/f) with h ∈ N;
• t₁ = k 1/(2f), with k even number, and any value of t₄;
• 2t₁ + t₄ = k (1/f) with k ∈ N, in general.

In antisymmetric non-short moves with f_b = f_a = f and t₁ = t₂ = t₃ = t₅ = t₆ = t₇, we have no vibrations when:

• t₁ = k 1/(2f), with k ∈ N, and any value of t₄;
• 3t₁ + t₄ = k (1/f) with k ∈ N, in general.

4. 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.

Position	First Natural Frequency Expected (Hz)	First Natural Frequency Recorded (Hz)	Err %
Position 1	14.4972	14.664	+1.15
Position 2	14.7856	14.969	+1.24
Position 3	13.6087	13.993	+2.82
Position 4	13.6305	13.993	+2.66
Position 5	13.8658	13.993	+0.92
Position 6	20.4669	20.827	+1.76
Position 7	20.5024	20.827	+1.58

Position	First Natural Frequency Expected (Hz)	First Natural Frequency Recorded (Hz)	Err %
Position 1	14.4972	14.644	+1.01
Position 2	14.7856	14.969	+1.24
Position 3	13.6087	13.993	+2.82
Position 4	13.6305	13.993	+2.66
Position 5	13.8658	13.993	+0.92
Position 6	20.4669	20.827	+1.76
Position 7	20.5024	21.152	+3.17

Position	First Natural Frequency Expected (Hz)	First Natural Frequency Recorded (Hz)	Err %
Position 1	14.4972	14.644	+1.01
Position 2	14.7856	14.969	+1.24
Position 3	13.6087	13.993	+2.82
Position 4	13.6305	13.993	+2.66
Position 5	13.8658	13.993	+0.92
Position 6	20.4669	20.827	+1.76
Position 7	20.5024	21.152	+3.17

Author Contributions

Conceptualization, A.A., V.P. and F.N.; methodology, A.A., A.S, V.P. and F.F.; software, A.A., F.F. and G.A.; validation, S.P., A.G.; investigation, A.A., V.P., F.F. and G.A.; writing-original draft preparation, A.A.; writing-review and editing, A.S., V.P. and R.S.; supervision, V.P., F.N.; project administration, R.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research received funding in the framework of research contract between Department of Industrial Engineering at University of Salerno and COMAU SpA, "Elastic models of robots and trajectory planning", Year 2018.

Acknowledgments

This paper collects a part of the results achieved in the framework of the above stated research contract.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

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 fa until the end of part 4 and, after, fb is the following:

_xfin(t)=1_{24fa5 fb5} ^π5 _{t12 t52[16fa5 fb5 Jmax1} ^π5 _{tt14 t52 +16fa5 fb5 Jmax1} ^π5 _{t15 t52 +16fa5 fb5 Jmax1} ^π5 _{tt13 t2 t52}+_{+24fa5 fb5 Jmax1} ^π5 _{t14 t2 t52 +8fa5 fb5 Jmax1} ^π5 _{t13 t22 t52 +16fa5 fb5 Jmax1} ^π5 _{t14 t4 t52 +16fa5 fb5 Jmax1} ^π5 _{t13 t2 t4 t52}+_{+32fa5 fb5 Jmax1} ^π5 _{t14 t53 +32fa5 fb5 Jmax1} ^π5 _{t13 t2 t53-16fa5 fb5 Jmax2} ^π5 _{tt12 t54-16fa5 fb5 Jmax2} ^π5 _{t12 t55}+_{+16fa5 fb5 Jmax1} ^π5 _{t14 t52 t6 +16fa5 fb5 Jmax1} ^π5 _{t13 t2 t52 t6-16fa5 fb5 Jmax2} ^π5 _{tt12 t53 t6-24fa5 fb5 Jmax2} ^π5 _{t12 t54 t6}+_{-8fa5 fb5 Jmax2} ^π5 _{t12 t53 t62-6fa5 fb Jmax2 πt12 t5 cos(2fb πt)-6fa5 fb Jmax2 πt12 t5 cos(2fb}π(_t+t5))++_{6fa5 fb Jmax2 πt12 t5 cos(2fb}π(_{t+t5 +t6})_{)+6fa5 fb Jmax2 πt12 t5 cos(2fb}π(_{t+2t5 +t6}))+_{-3fa2 fb4 Jmax1 πt1 t52}cos(2π(_{-fa t4 +fb}(_{t+2t5 +t6})))+_{+3fa fb5 Jmax1 πt1 t52}cos(2π(_{-fa t4 +fb}(_{t+2t5 +t6})))+_{+3fa2 fb4 Jmax1 πt1 t52}cos(2π(_{fa t4 +fb}(_{t+2t5 +t6})))+_{+3fa fb5 Jmax1 πt1 t52}cos(2π(_{fa t4 +fb}(_{t+2t5 +t6})))+_{-3fa2 fb4 Jmax1 πt1 t52}cos(2π(_-fa(_{t1 +t4})_+fb(_{t+2t5 +t6})))+_{+3fa fb5 Jmax1 πt1 t52}cos(2π(_-fa(_{t1 +t4})_+fb(_{t+2t5 +t6})))+_{+3fa2 fb4 Jmax1 πt1 t52}cos(2π(_fa(_{t1 +t4})_+fb(_{t+2t5 +t6})))+_{+3fa fb5 Jmax1 πt1 t52}cos(2π(_fa(_{t1 +t4})_+fb(_{t+2t5 +t6})))+_{+3fa2 fb4 Jmax1 πt1 t52}cos(2π(_-fa(_{t1 +t2 +t4})_+fb(_{t+2t5 +t6})))+_{-3fa fb5 Jmax1 πt1 t52}cos(2π(_-fa(_{t1 +t2 +t4})_+fb(_{t+2t5 +t6})))+_{-3fa2 fb4 Jmax1 πt1 t52}cos(2π(_fa(_{t1 +t2 +t4})_+fb(_{t+2t5 +t6})))+_{-3fa fb5 Jmax1 πt1 t52}cos(2π(_fa(_{t1 +t2 +t4})_+fb(_{t+2t5 +t6})))+_{+3fa2 fb4 Jmax1 πt1 t52}cos(2π(_-fa(_{2t1 +t2 +t4})_+fb(_{t+2t5 +t6})))+_{-3fa fb5 Jmax1 πt1 t52}cos(2π(_-fa(_{2t1 +t2 +t4})_+fb(_{t+2t5 +t6})))+_{-3fa2 fb4 Jmax1 πt1 t52}cos(2π(_fa(_{2t1 +t2 +t4})_+fb(_{t+2t5 +t6})))+_{-3fa fb5 Jmax1 πt1 t52}cos(2π(_fa(_{2t1 +t2 +t4})_+fb(_{t+2t5 +t6})))+_{-6fa5 Jmax2 t12}sin(_2fbπt)_{+6fa5 Jmax2 t12}sin(_2fbπ(_t+t5))_{+6fa5 Jmax2 t12}sin(_2fbπ(_{t+t5 +t6}))+_{-6fa5 Jmax2 t12}sin(_2fbπ(_{t+2t5 +t6}))_{+3fa fb4 Jmax1 t52}sin(2π(_{-fa t4 +fb}(_{t+2t5 +t6})))+_{-3fb5 Jmax1 t52}sin(2π(_{-fa t4 +fb}(_{t+2t5 +t6})))_{+3fa fb4 Jmax1 t52}sin(2π(_{fa t4 +fb}(_{t+2t5 +t6})))+_{+3fb5 Jmax1 t52}sin(2π(_{fa t4 +fb}(_{t+2t5 +t6})))_{-3fa fb4 Jmax1 t52}sin(2π(_-fa(_{t1 +t4})_+fb(_{t+2t5 +t6})))+_{+3fb5 Jmax1 t52}sin(2π(_-fa(_{t1 +t4})_+fb(_{t+2t5 +t6})))_{-3fa fb4 Jmax1 t52}sin(2π(_fa(_{t1 +t4})_+fb(_{t+2t5 +t6})))+_{-3fb5 Jmax1 t52}sin(2π(_fa(_{t1 +t4})_+fb(_{t+2t5 +t6})))+_{-3fa fb4 Jmax1 t52}sin(2π(_-fa(_{t1 +t2 +t4})_+fb(_{t+2t5 +t6})))+_{+3fb5 Jmax1 t52}sin(2π(_-fa(_{t1 +t2 +t4})_+fb(_{t+2t5 +t6})))+_{-3fa fb4 Jmax1 t52}sin(2π(_fa(_{t1 +t2 +t4})_+fb(_{t+2t5 +t6})))+_{-3fb5 Jmax1 t52}sin(2π(_fa(_{t1 +t2 +t4})_+fb(_{t+2t5 +t6})))+_{+3fa fb4 Jmax1 t52}sin(2π(_-fa(_{2t1 +t2 +t4})_+fb(_{t+2t5 +t6})))+_{-3fb5 Jmax1 t52}sin(2π(_-fa(_{2t1 +t2 +t4})_+fb(_{t+2t5 +t6})))+_{+3fa fb4 Jmax1 t52}sin(2π(_fa(_{2t1 +t2 +t4})_+fb(_{t+2t5 +t6})))+_{+3fb5 Jmax1 t52}sin(2π(_fa(_{2t1 +t2 +t4})_+fb(_{t+2t5 +t6})))]