A Multicriteria Motion Planning Approach for Combining Smoothness and Speed in Collaborative Assembly Systems

: Human–robot interaction is an important aspect of Industry 4.0, and the extended use of robotics in industrial environments will not be possible without enabling them to safely interact with humans. This imposes relevant constraints in the qualitative characterization of the motions of robots when sharing their workspace with humans. In this paper, we address the trade-off between two such constraints, namely the smoothness, which is related to the cognitive stress that a person undergoes when interacting with a robot, and the speed, which is related to normative safety requirements. Given an execution time, such an approach will allow us to plan safe trajectories without neglecting cognitive ergonomics and production efﬁciency aspects. We ﬁrst present the methodology able to express the balance between these qualities in the form of a composite objective function. Thanks to the variational formalism, we identify the related set of optimal trajectories with respect to the given criterion and give a suitable parametrization to them. Then, we are able to formulate the safety requirements in terms of a reparametrization of the motion. Finally, numerical and experimental results are provided. This allows the identiﬁcation of the preferable sets of the possible motions that satisfy the operator’s psychological well-being and the assembly process performance by complying with the safety requirements in terms of mechanical risk prevention.


Introduction
Collaborative robots (cobots) are robots designed to enable the combination of the strength, force, and repeatability of automatic machines with the dexterity and cognitive capacities of humans, allowing them to work together in a shared workspace. It is expected that they will boost the Industry 4.0 revolution and, in particular, the digitalization and the human-centered design of an industrial workspace [1].
However, the physical human-robot interaction (pHRI) required in such a shared workspace introduces novel mechanical risks with respect to "traditional" industrial robots [2]. In this regard, the International Organisation for Standardization (ISO) has introduced four methods of operation [2]: i.e., hand guiding, safety-rated monitored stop, speed and separation monitoring, and power and force limiting. Such modalities constitute a framework to define a cobot, proposing proper implementations and evaluations of safety requirements in terms of mechanical risk.
Commercially available solutions provide out-of-the-box user interfaces and tunable parameters to allow the integrator to implement these collaborative operation methods [3][4][5].
In addition to safety, other important topics that have to be considered when designing and implementing a collaborative workspace are physical and cognitive ergonomics. For the former, we can utilize industry standards, guidelines, and best practices already developed for manual work and modify and apply them to collaborative scenarios [6,7]. On the contrary, the latter has not yet been fully studied, and guidelines and best practices are not available [8].
It has been demonstrated that the operator's subjective acceptance of pHRI can be increased by means of smooth, i.e., predictable, motions [9][10][11][12]. We therefore formulate a robot motion planning approach to achieve a psychologically less-stressful working environment.
Among all motion planning techniques that are available in the literature, trajectory planning through via-points is of particular importance for collaborative robotics. The Technical Specification ISO TS-15066 [13] foresees the hand guiding modality as collaborative operation: even unskilled workers can intuitively program a path by simply moving the robot manually in its workspace and saving step-by-step the needed via-points through an interface. The worker can thus directly define the robot's configuration and avoid the obstacles that may be placed in the reachable workspace. The computation of a safe trajectory along the via-points is then made by the robot controller according to the ISO specifications.
Starting from the outcome of the hand-guiding operation and path definition, this work focuses on the problem of designing trajectories that allow guaranteeing both safe and psychologically less stressful movements without overlimiting the robot performance.
In the literature, there is a plethora of different techniques for planning a robotic trajectory [14] that focus on the optimization/minimization of, e.g., time, smoothness, energy, and torque. Among these, minimum jerk motion profiles have drawn particular interest in the robotics community for mainly two reasons:

•
There is experimental evidence that suggests that humans adopt minimum jerk movements [15][16][17], and thus, these can be perceived as predictable and familiar by a person that interacts with the robot. In [18], the impact of different industrial robot motion profiles was investigated in a cooperative human-robot interaction (HRI) task showing that minimum jerk trajectories significantly reduce the heart rate in humans, thus suggesting an improvement of robot acceptance and a reduction of the psychological stress; • Minimum-jerk trajectories reduce the induced mechanical vibrations and, thus, wear in the elements of an automatic machine, e.g., robot, as well as the error during the trajectory tracking phase [19][20][21][22][23][24].
Thus, we conclude that minimum jerk trajectories allow psychologically less stressful movements with respect to well-known ones like trapezoidal speed profiles. However, as pointed out by the authors in [25], minimum jerk motions could have two potential drawbacks that degrade the robot's behavior: • For a given execution time, minimum jerk motions could require higher speeds than trapezoidal speed profiles. To make the motion feasible for the robot and for the safety requirement, an upper bound value of velocity has to be set, and a linear scaling of the trajectory penalizes this kind of trajectory in terms of execution time much more than it does for trapezoidal speed profiles; • Large deviations from the straight lines that join consecutive via-points may occur. The magnitude of such deviations is not predictable and may create unexpected deviations from the intended path between via-points, forcing the user to define additional via-points to possibly limit this problem.
In this paper, we overcome these issues by formulating a multi-criteria optimization problem in order to take into account both the jerk and the speed along the path that joins a sequence of via-points. Following the approach developed in [25], we directly apply the variational formalism to derive a family of trajectories that are parametrized with respect to a parameter that expresses the required balance. As a result, we obtain a human-friendly motion planning methodology, the input of which is given by the via-points introduced during the hand-guiding programming phase, and the output is a reference trajectory for the robot's controller.
In Section 2, we formulate the problem of trajectory planning through via-points and demonstrate that the previously mentioned drawbacks of minimum jerk motions are consequences of a naturally large arc-length. Then, we formulate a multi-criteria (multi-objective) optimization problem that accounts for balance between arc-length and jerk. This section finishes with the application of the variational formalist to the proposed problem. In Section 3, we provide the means to implement this kind of motion in a collaborative environment, and we leverage a Pareto front to provide the user the necessary information to choose the balance parameter for the given execution time. Finally, in Section 4, the experimental results of the implemented trajectories are given.

General Problem Formulation
In order to exploit the benefits of a smooth, i.e., minimum jerk, trajectory, and limit the possible overshoots in the path, a trade-off method between the two objectives is investigated and proposed.
To obtain a uniform description of the motion regardless of the space in which the motion is planned, e.g., the joint space of a specific task, we describe it as a curve in an open subset Ω of an n-dimensional manifold M such that there exists a chart (ψ, Ω) such that: This allows identifying the motion with a curve in R n without considering the specific details of the space where the planning is done.
Let us consider a partial description of a desired robot motion given as the execution time t f and an ordered sequence of N + 1 via-points in the manifold.
Utilizing the previously introduced local coordinate map ψ, Equation (1), the problem can be rewritten in R n , and the ordered sequence of desired via-points can be described as This operation translates our problem into that of finding a reasonable curve joining the via-points VP in such a way that the via-point w i and the velocityẇ i are attained at the time t i with t i < t i+1 in which every instant t i is within the execution time interval T = [0, t f ]. Then, the position, velocity, and acceleration constraints can be written as: with t 0 = 0 and t N = t f .
We define the term reasonable curve as a curve that optimizes a desired cost function and is "physically" feasible. As the cost functions of this work include the integrals of the derivatives of q, it is required that the k-th first derivatives are square integrable. Such a set of functions constitutes a Hilbert space, denoted by W 2,k (T, R n ). The second requirement demands that q connects the via-points VP with a trajectory that is continuous to the k-th order. Furthermore, such variations must be attainable by the dynamics of the systems.
Thus, we introduce the set of admissible motions as: The problem can then be seen as the one of finding the curve among the admissible motions Q such that a given performance index is minimized and additional possible constraints are satisfied.
We define the performance index as a map I : Q −→ R that describes a desired quality of the motion under analysis by associating it with a real number. A simple way to construct such a performance index is to take the average value of a given scalar function F, which depends on the desired i-th derivative. The cost function I can then be expressed as: We can then express the optimal trajectory planning problem in the following form: to which further constraints can be added. The design problem is therefore often formulated in terms of minimum velocity F = q , minimum acceleration F = q , and/or minimum jerk F = ... q . Therefore, this is a multi-criteria (or multiple objective) design problem. There are a number of techniques to handle multiple objectives (or cost functions; see, e.g., [26]). In engineering cases where computational effort is high or where the computational time should be kept to a minimum, the problem is transformed into a single-objective problem by either constraining all objectives except one, i.e.: or, as in this work, by creating a composite objective function using weighting factors: As such, the problem can then be solved with a wide range of optimization techniques. In our specific case, a min-jerk trajectory could lead to a very large arc-length. In fact, since the curvature is proportional to the second derivative, the jerk minimization creates a smooth trajectory with low curvature values.
As the curvature is the rate of change of the direction of the curve with respect to the arc-length, it follows that a minimum jerk trajectory will require a curve with a significant arc-length to join a sequence of non-collinear via-points. We can conclude that minimum jerk motion has a larger arc-length than other kinds of trajectories such as the trapezoidal velocity profile. It follows that a limit on the velocity will penalize a minimum jerk motion more than other motions in terms of execution time. To overcome these issues, we formulate a problem that grasps the balance between smoothness and arc-length.
Our approach consists of recognizing that the velocity cost function is proportional to the arc-length. It follows that we can construct a cost function that integrates the smoothness and arc-length with a convex combination of the velocity and the jerk: where α varies from zero, i.e., only considering the jerk-related objective, to unity, i.e., only considering the arc-length related objective.

The Variational Formalism
The variational formalism of the optimization of (4) requires the formal definition of a small variation of a curve q ∈ Q. Following a commonly adopted approach [27,28], we construct a small variation of a curve q ∈ Q: by using an arbitrarily smooth map η ∈ C ∞ (T, R n ) called variation and with an arbitrarily small number ε > 0. Both curves p and q satisfy the constraints (2) at their respective set of times t 0 , . . . , t N and r 0 , . . . , r N such that q(t i ) = q i and p(r i ) = q i . By definition, the difference between t i and r i is arbitrarily small and depends on ε through the relation [28] r i = t i + εδt i , where the scalar δt i is an arbitrary variation of t i . Furthermore, due to the small differences in the values of velocities and acceleration at the prescribed points, consequently, we have to consider the arbitrary variations . Such variations have a concrete relation with the variations η(t) given by [28], The variations δt i , η(t i ), δq i , and δq i vanish for i = 0 and i = N due to the constraints in (2). Thanks to the formal definition of small variations stated above, we can disclose the necessary conditions that make an element of Q a stationary point of (4). Such elements are called extremals of the functional (4) and are characterized by the relation: where O(ε 2 ) represents a negligible quantity proportional to ε 2 called big-O notation. Such a relation may be expanded at each interval as: By considering the relations between t i and r i , we have: By integrating by parts and neglecting the terms proportional to ε 2 , this relation can be expanded to: After substituting (10) in (14), we can rearrange such an equation as a linear combination of the variations η, δt i , δq i , δq i . Such variations being arbitrary, they are linearly independent [27].
It follows that the necessary stationarity conditions of I at q are defined in the intervals where t − i and t + i represent the limits of t when approaching t i from the left and the right, respectively, and (15) are the Euler-Lagrange equations at each interval J i .

Final Formulation
The necessary optimality conditions (15)- (17) state that the optimal of (4) is constituted by a set of optimal arcs that join the via-points and satisfy (15) in each interval. An explicit formula of these optimal arcs may be easily written by introducing the set J 0 = [−1, 1], the vector τ ∈ R N with components τ i = t i+1 − t i , and the family of linear dilations s i : By applying the transformation (18) to (15), we obtain that at each interval J i , its solution is given by the solution of: If follows that the optimal arcs that join the via-points and satisfy the Euler-Lagrange equations at the interior of each interval J i can be written in the nominal interval J 0 as a linear combination of the general solutions of (19). Following the traditional approach for solving ordinary linear differential equations [29], we consider a guest solution of the form: where a is an arbitrary column vector and β an unknown scalar to be found. We can confirm that the general solution of (19) is constituted by a linear combination of functions of the type (20) for which: After substituting the solutions of (21) into (20) and selecting the linear combinations that have only the real part [29], we have that the general solution of (19) is: where: and: Using the general expression (22), we can uniquely represent the i-th component at the interval J j of any extremal of (4) using a vector y j i ∈ R 6 and the length of the interval τ j as: and (25) represents the general formula. It is therefore possible to achieve a uniform representation of the extremals of (4) in the whole execution interval T with a vector y = [(y 1 1 ) , . . . , (y N n ) ] of dimension n p = 6nN and the vector τ. Using this notation, we rewrite (4) as: where Q(τ) is a 6nN × 6nN block diagonal matrix with 6 × 6 blocks Q (τ) of the form: where k = (j mod 6n) + 1. However, we can prove that the representation of the extremals of (4) using the pair (y, τ) is redundant. In fact, Constraints (2) and (16) may be written as: where b is a vector containing the via-points and the boundary conditions in (2). Thanks to (28), we can describe the problem of the optimization of (4) as an optimization problem in N real variables contained in the vector τ. In fact, the problem of minimizing (4) subject to (2) can be formulated as: where 1 ∈ R N is the all-ones column vector. This formulation has a convex feasible domain, and due to the inversion of (28), its computational complexity is proportional to O(n 3 N 3 ).

Implementation
The problem (29) and its resulting family of minimum jerk and velocity trajectories requires the definition of the execution time t f and of a weighting factor better α. This is thus different with respect to the minimum jerk problem and its solution developed in [25], which provides a family of trajectories independent of the execution time. Threemain steps are then to be addressed for choosing the optimal trajectory: 1.
calculation of the family of trajectories dependent on the weighting factor α; 2.
multi-criteria optimization finding the optimal time interval τ i and selection of the weighting factor α; 3.
trajectory feasibility verification and possible scaling.

Multi-Criteria Optimization and Pareto Front
As introduced in Section 2.1, we use a composite objective (cost) function (8). This can be also solved in a series of optimizations with varying weighting factors α, resulting in a set of optimal designs based on these weighting factor values. Each optimal design belongs to the Pareto optimal set, which is also referred to as a Pareto front. The Pareto front allows the worker or engineer to decide on the trade-off between the different objectives represented in the composite function, here smoothness and arc-length. Therefore, by exploiting this method, we avoid having to choose a value for the weighting factor α a priori.
The Pareto front is shown here as discretely approximated for each value of the weighting factor α. Due to differences in magnitude between the objective functions of smoothness and arc-length, a logarithmic distribution can be chosen to ensure proper spreading of the points and, therefore, a well-approximated Pareto front. The set of constrained optimization problems leading to the Pareto front then becomes: where where F is the cost function in (29) and τ i are the time intervals between via-points on the trajectory found via the variational problem. It should be noted that the equality constraint on time, requiring that the operation be completed in the desired amount of time, is implemented as two inequality constraints, The efficient second-order algorithm NLPQLP [30] is utilized to solve this optimization problem posed in (30)- (32). This optimization formulation, though, is general, and therefore, any gradient-based optimization algorithm can be applied.
In Section 4, it will be shown that with a relatively low number of points, the Pareto front is well defined and can be used effectively by the worker or engineer to make design decisions. It should be further pointed out that this compromise decision is to be made by a worker, and therefore, more than two or three objectives are typically not able to be handled (unable to be readily shown graphically) and therefore not advisable.

Feasibility of the Motion
Given an execution time t f fixed a priori, the trajectory resulting from an optimization from certain value α may result in a non-feasible trajectory. Therefore, in order to find a feasible trajectory in terms of robot limits and mechanical safety, it may have to be scaled in order to have a feasible design. However, since the arc-length does not change after a re-parametrization, we can state that for any α > 0, the solutions of (29) may be linearly scaled without losing their property of being minimum jerk trajectories.
If the planned motion violates robot limits or mechanical safety standards, we apply a linear dilation σ : C 0 (R, R n ) −→ C 0 (R, R n ) to the computed trajectory given by: where σ is a function that takes the continuous vector-valued trajectory q and returns the same trajectory scaled by the scaling factor σ 0 . The mechanical safety is considered by using the concept of a danger field [31], which is a map DF : R 3 × R n × R n −→ R such that DF(x, q,q) represents a measure of the hazard at the position x induced by a robot, which has its joints at positions q moving with speedsq. A robotic application is considered safe if the value of the danger field is below a given threshold DF max determined during the application risk assessment. In particular, we consider danger fields that are homogeneous with respect toq: DF(x, q, σ 0q ) = σ 0 DF(x, q,q).
The limits of the robots are considered by taking into account the maximum values of velocity V max and acceleration A max achievable by the robot and retrieved by the data sheets and manuals provided by the robot manufacturer. If follows that, to make a motion feasible, we have to compare the given feasibility thresholds with respect to the following values: Then, if the planned movement violates the limits of the robot or mechanical safety, the linear scaling to assure the feasibility of the motions is:

Implemented Procedure
The final procedure implemented for running the developed methodology is as follows:

1.
User: definition of the trajectory via-points and desired execution time.

2.
Robot controller: computation of the optimal solutions and generation of the Pareto front chart.

3.
User: selection of the desired trajectory choosing the proper value of α from the Pareto front chart.

4.
Robot controller: evaluation of the trajectory feasibility and, if needed, linear scaling of it.

5.
Robot controller: execution of the trajectory according to the starting command given by the user.

Experimental Setup
The developed method and procedure was tested considering a six degree-of-freedom UR3 collaborative robot from Universal Robots available at the Smart Mini Factory laboratory of the Free University of Bozen-Bolzano [32]. The robot is an anthropomorphic robot with a non-spherical wrist, and the main characteristics are reported in Figure 1.
This robot is part of a collaborative work-cell developed by making a manual assembly operation for assembling pneumatic cylinders collaborative, as described in [33]. The operation workflow of the station consists of a pick-and-pass operation where the robot passes the required components of the assembly to the user. These components are organized a priori in a blister, and the robot trajectories are computed with our approach using the user-provided via-points. Following ISO TS-15066 [13], a mechanical risk assessment was performed [33], and the computed velocity limits are stated in Table 1. In order to simplify the expression of the danger field, we suppose the source of the hazard as concentrated at the robot's gripper. Then, the safety limits are encoded in a danger field as: whereê o ∈ R 3 is the direction pointing towards the operator, x is the Cartesian position of the gripper, k : R n −→ R 3 is the direct kinematics of the gripper, A is the volume of the collaborative zone, and J is the Jacobian matrix of the gripper's position. In our specific case, we haveê o = [0, −1, 0] . To follow the trajectories, we implement a PDvelocity controller of the form: where u represents the control command, q (t) the computed trajectory, and (q(t),q(t)) the robot state.

Upper arms and elbow joints 330
Lower arms and wrist joints 360

Trajectories Chosen for Numerical and Experimental Comparison
In order to evaluate and validate the proposed method, two benchmark examples representing two qualitatively different motions, as well as a pick-and-pass movement for a collaborative activity have been chosen.

Benchmark Trajectories
The benchmark trajectories represent different shapes drawn in the horizontal plane with the end-effector; their via-points in the operational space are reported in Figure 2 and correspond to a zigzag pattern (eight via points) and to a figure-eight symbol (19 via-points), respectively. The shape of these trajectories may be appreciated in Figure 2. The relative joint configurations, obtained via a kinematic inversion on the given path in the operational space, are reported in radians in Table 2.

Collaborative Trajectory
The collaborative operation consists of a pick-and-pass task composed by two consecutive movements. In the former, the robot starts from a position close to the operator and moves close to the component to be picked, approaching it for a proper grasping. Once the gripper has performed its grasping action, the latter movement is executed: the robot moves back the starting position in order to pass the object to the operator. The via-points were introduced through a hand-guiding programming phase and are reported in Table 3. The required execution time given by the particular production process was 14 s.

Benchmark Trajectories
The zigzag and figure-eight considered benchmark trajectories were evaluated for execution times of t f = 10.5 and t f = 19.0 s, respectively. Given the path and time of execution, the Pareto optimal sets were computed. Eleven discrete values of the weighting function α and therefore eleven separate optimizations were carried out considering the composite objective function and a logarithmic distribution used for a well-approximated Pareto front, For the two sets of via-points, the optimization procedure utilized as a first guess τ 0 = t f /N1. To reduce computational effort, each subsequent optimization started at the previous optimum. In this study, the number of evaluations was hereby able to be reduced by approximately 30% for both cases. Further, as both sets of starting points led to the exact same results, we assumed that gradient-based optimization was effective (without commenting on the convexity of the problem). The eleven optimization runs resulting in the Pareto front for the zigzag trajectory found in Figure 3 (left) needed 98 evaluations and 72 analytical sensitivity evaluations, requiring 7.37 s on a four core Intel i5-7200U CPU with 2.50 GHz.
Physically, this solution lengthened the time on the short ends of the zigzag trajectory and shortened the time spent on the longer ends. The L 2 norm of the jerk was reduced, while increasing the L 2 norm of the velocity.
An 11 point Pareto front was also calculated for the figure-eight trajectory, and this can be seen in Figure 3 (right).Due to the higher number of via-points and therefore design variables, the complexity of the optimization was increased. Here, four-hundred thirty-eight system evaluations and 179 analytical sensitivity evaluations were required and, in turn, 65 s on the same four core Intel i5-7200U CPU with 2.50 GHz.
Being that the figure-eight path had no edges, the optimization had less effect on the arc-length and jerk values with respect to the previous benchmark trajectory. However, the jerk values can be greatly modified according to the selection of the weighting parameter.  To complement these results, in Figures 4 and 5, the numerical results obtained for the two trajectories with the weighting factor α → 0, i.e., min-jerk trajectory in red, and with α → 1, i.e., arc-length minimization in red, are reported. As one might have expected, the red curve shows less overshoots with respect to the red one, in particular in the case of an "edgy" path.

Collaborative Trajectory
Considering the collaborative operation and the fact that it can be seen as a sequence of two movements, two Pareto fronts were computed, and the results are presented in Figure 6. The obtained trajectories (in particular, the min-jerk one) were verified taking into account the speed limits, i.e., Table 1, and a danger field (37) given by the volume with y < −10 cm. Thanks to these, the user can select the proper weight according to the needs and kind of movement and give to the robot controller the command for starting the execution of the chosen optimal trajectory.  To better show the developed method's capabilities, a comparative picture was realized by computing and running three types of trajectories for the given sequence of via-points: • a trapezoidal velocity profile generated by the robot's controller; • a min-jerk trajectory computed considering α → 0; • a minimum arc-length trajectory computed considering α → 1.
A comparison of the three results is reported in Figures 7 and 8 where the robot's trapezoidal velocity profile trajectory is depicted in blue, the proposed approach is depicted with red, α → 0 with a dashed line and α → 1 with a solid line. Figure 7 shows the comparison in the joint space for the first three joints of the robot, while Figure 8 reports a comparison in the operational space in the directory of the operator (y-axis; see Figure 1). We underline that the constraints in Table 1 are respected for the given trajectories in Figure 8.    The proposed approach is inherently smoother than the robot's trapezoidal speed profile and possibly results in being less stressful for the operator.
In Table 4, we present a comparison between the maximum attained velocities. The trajectory obtained with the proposed approach with α → 1 showed a shorter path with respect to the one with α → 0 and the lowest maximum velocity. As a consequence, from the mechanical-risk point of view, the trajectory was safer with respect to the others, confirming the validity of the approach.

Conclusions
In this study, we address the trade-off between two constraints, namely the smoothness and speed, via a multicriteria trajectory planning approach. Smoothness is related to the cognitive stress that a person undergoes when interacting with a robot, while speed is related to safety requirements. Given an execution time, we developed the proposed method to plan safe trajectories without neglecting cognitive ergonomics and production efficiency aspects.
The balance between these two qualities was achieved through a composite objective function, and thanks to the variational formalism, a family of optimal trajectories, as well as a compact formulation for the optimization problem were presented. Thanks to the Pareto front, the operator was able to choose a trade-off between smoothness and the drawbacks of minimum jerk and obtain the trajectory that minimized the chosen cost function. Numerical and experimental results showed the feasibility of the approach, as well as the improvement of the performances both with respect to the native planner, i.e., trapezoidal speed profile, and min-jerk trajectories. Future work will cover: (i) further and extensive comparison and evaluation of the proposed method considering also the dilation effect, as well as (ii) its integration as a path option in the robot software and/or as an online service for the calculation of the optimal path and subsequent transfer to the robot.
Author Contributions: Conceptualization, R.A.R.; methodology, all authors; software, R.A.R. and E.W.; validation, R.A.R. and E.W.; writing, original draft preparation, all the authors; writing, review and editing, all authors; supervision, R.V. All authors read and agreed to the published version of the manuscript.