Abstract
A new methodology using a direct method for obtaining the best found trajectory planning and maximum dynamic load-carrying capacity (DLCC) is presented for a 5-degree of freedom (DOF) hybrid robot manipulator. A nonlinear constrained multiobjective optimization problem is formulated with four objective functions, namely, travel time, total energy involved in the motion, joint jerks, and joint acceleration. The vector of decision variables is defined by the sequence of the time-interval lengths associated with each two consecutive via-points on the desired trajectory of the 5-DOF robot generalized coordinates. Then this vector of decision variables is computed in order to minimize the cost function (which is the weighted sum of these four objective functions) subject to constraints on joint positions, velocities, acceleration, jerks, forces/torques, and payload mass. Two separate approaches are proposed to deal with the trajectory planning problem and the maximum DLCC calculation for the 5-DOF robot manipulator using an evolutionary optimization technique. The adopted evolutionary algorithm is the elitist nondominated sorting genetic algorithm (NSGA-II). A numerical application is performed for obtaining best found solutions of trajectory planning and maximum DLCC calculation for the 5-DOF hybrid robot manipulator.
1. Introduction
The most primary consideration in designing and accomplishing a fundamental required task of a robot manipulator is to plan the trajectory for this robot manipulator moving from its current placement to a series of specified placements while carrying a payload. During these operations of the robot manipulator, the problem of finding such a path is referred to as the trajectory planning problem. Usually, solving an optimization problem is defined to deal with the trajectory planning problem. There are a number of popular methods for transcribing the trajectory planning problem into an optimization problem. Essentially, these transcriptions fall into two categories: direct (global) methods (such as [1–4]) and indirect (decoupled) methods (e.g., [5–8]).
Solving trajectory planning for the robot manipulator is difficult because of the highly coupled nonlinear robot manipulator dynamics and complex system constraints (which include limitations of actuator forces/torques and payload constraint). A feasible and efficient method must be able to find an appropriate trajectory such that it minimizes a multicriterion cost function and satisfies the limitations. This multicriterion cost function can include travel time, energy, actuator torques, and joint jerks. In order to achieve the greatly improved productivity in industrial situations, it is necessary to minimize the total traveling time of the robot manipulator. Many research works have concentrated on minimum time trajectories [9–13]. The robot manipulator trajectory planning using the energy criteria is reported in several literatures [14–18]. This kind of manipulator trajectory planning can produce smooth trajectories and reduce the stresses between the actuators and the robot manipulator structure. Moreover, it can save energy, which may be desirable in several applications, such as aerial, space, or submarine robots with a capacity limitation of energy source. Both optimal traveling time and the minimum energy of the actuators are taken together as objective functions in some literatures [19–22]. Several examples of minimizing joint jerks and acceleration in the objective function are studied in [23–27]. Minimization of the joint jerks and acceleration usually contributes to the reduction of errors during trajectory tracking and stresses between the actuators and structure of the robot manipulator. Also, excitation frequencies of the robot manipulator can be limited.
Many attempts to determine the optimal trajectory of robot manipulators are mainly conventional optimization techniques such as sequential unconstrained minimization technique, sequential quadratic programming, dynamic programming, and numerical iterative procedure. However, these conventional techniques have some drawbacks, such as easily falling into local minima or involving calculating the gradient and the Hessian of the objective function and constraints. Usually, traditional techniques cannot solve discontinuous terms (e.g., friction) in physical models. In order to overcome these drawbacks, the evolutionary methods can be applied to handle the optimization trajectory problems [23, 28]. There are several advantages of using evolutionary techniques, which include (a) population-based approach search, thus discovering the global optimal solution is possible, (b) complex and multimodal problems which can be solved for global optimality, (c) no need to provide any auxiliary relationships (like gradients, derivatives, etc.), and (d) problem independence; that is, they are suitable for solving all types of optimization problems.
The dynamic load-carrying capacity (DLCC) is another important term, which is an important characteristic of the robot manipulators. It is mandatory to consider the load-carrying capacity for a robot manipulator and the payload constraint for the trajectory planning. DLCC has received considerable attention in research literature. Wang and Ravani [29, 30] determined the DLCC of a given trajectory for a serial manipulator as an optimization problem according to the state space of the dynamics equation and the torque-speed features of the actuators. Korayem et al. proposed several methods in order to find DLCC of a wheeled mobile manipulator, such as a method which is a combination of open and closed loop optimal control [31], a method based on hierarchical optimal control [32], and a method using open-loop optimal control approach and considering stability on zero-moment point criterion [33].
Several classical optimization approaches have been used to deal with the trajectory planning and the DLCC calculation. Dynamic programming approach is applied to solve optimization problem with nonlinear constraints and the sequential quadratic programming formulation is used to increase the speed of convergence for both fixed and mobile base mechanical manipulators [34]. The iterative linear programming method and the discrete augmented Lagrangian technique are used to solve the nonlinear constrained optimization problem for cable driven robots [35] and parallel kinematic machines [21], respectively.
In this paper, a nonlinear constrained multiobjective optimization problem is formulated with four objective functions (namely, travel time, total energy involved in the motion, joint jerks, and joint acceleration), considering constraints of joint positions, velocities, acceleration, jerks, forces/torques, and payload. The sequence of the time-interval lengths associated with each two consecutive via-points is defined as the vector of decision variables. Then the cost function (which is the weighted sum of these four objective functions) is minimized by optimizing this vector of decision variables. Using a direct method, two approaches are proposed to deal with the trajectory planning problem and the calculation of maximum dynamic load-carrying capacity for a 5-DOF robotic manipulator, respectively. An evolutionary optimization technique, that is, an elitist nondominated sorting genetic algorithm (NSGA-II), is adopted to handle these optimization problems, taking into account these objective functions and constraints. The NSGA-II algorithm was originally proposed by Deb et al. [36]. It has been proved to be a highly effective and generally successful choice for optimization problems and becomes one of the mainstream solutions to multiobjective optimization questions [37–39].
The primary advantage of the proposed approaches in this paper is that the vector of decision variables (i.e., the sequence of the time-interval lengths associated with each two consecutive via-points) is optimized step by step during the evolution process (where the unfit individuals are replaced by the fit individuals through the crossover and mutation operator in each generation). The best found solutions of the joint trajectory displacement, velocity, acceleration, jerk, and force/torque can be obtained by these time-interval lengths. Additionally, the time intervals with variable lengths are more advantageous than the ones with same lengths (which are adopted by some classical optimal control approaches) because the former is more relatively likely to generate the smooth trajectory and avoid the suddenly change of the trajectory. Other advantages of the proposed approaches include the following: (a) both kinematics and dynamics aspects of the robot manipulator are considered, (b) they take several important decision criteria into account for trajectory planning (which are travel time, total energy involved in the motion, joint jerks, and joint acceleration), (c) adopting the different weights for the decision criteria, the required and suitable solutions that will satisfy various purposes can be obtained, (d) they consider payload mass constraint, (e) they ensure that the resulting best found trajectory is smooth and leads to minimum actuator power and energy, (f) they are easy to program and implement, and (g) they can also be extended to obtain the best found trajectory and the maximum DLCC for other types of robots.
The rest of this paper is organized as follows. The kinematics and dynamics of the studied 5-DOF robotic manipulator are described in Section 2. The problem formulation is presented in Section 3. An approach is proposed for the trajectory planning problem of the robot manipulator in Section 4. The maximum DLCC of the robot manipulator is calculated by a proposed approach in Section 5. The effectiveness of two proposed approaches for the trajectory planning and the maximum DLCC calculation is verified via a numerical example of this 5-DOF hybrid robot manipulator in Section 6. Finally, conclusions are drawn in Section 7.
2. Kinematic and Dynamic Modelling of the 5-DOF Robotic Manipulator
In this section, the kinematic and dynamics modelling of a 5-DOF robot manipulator are described. The structure of this robotic manipulator is shown in Figure 1. This robotic manipulator, consisting of a 3-DOF (3T) parallel module and a 2-DOF (2R) serial module, is a new 3T2R hybrid robot manipulator. Translational motion axes , , and (consisting of ball screws) are driven by motors , , and , respectively. Rotational motion axes and are driven by motors and , respectively. The axis is the end-effector which contains a high-speed motor spindle. More details of the structure description for this robot manipulator are described in [40].
(a)
(b)
(c)
2.1. Kinematic Modelling
Based on the Denavit-Hartenberg method [41], the reference frames and the kinematic modelling for this robot manipulator are presented.
2.1.1. Forward Kinematics
The reference frames for this robot manipulator are assigned as shown in Figure 2(a). The vector (shown in Figure 2(b)) is the location of the center of mass of a link relative to the frame representing the link. The rotation angle of , , , , and is represented by , , , , and , respectively. The stroke of the nut along , , and is represented by , , and , respectively. The length of the rod is denoted by . The distance between and in the horizontal direction is denoted as . The rotation angle of and link 1 is denoted as . represents the distance between point and point in the vertical direction when the nut is located at the top starting point along .
(a)
(b)
The detailed structure parameters of this robot manipulator are listed in Table 1. Considering the reduction ratio of the speed reducer, the pitch, the nominal torque of the AC servomotor, and the suitable safety coefficients, the constraints of displacement, velocity, acceleration, jerk, and force/torque (i.e., generalized force) are given in Table 2. Some structure parameter variables are related bywhere the constant represents the pitch of and , the constant represents the pitch of , and the available range of the rotation angle is because of the mechanical constraints.
Link 1 is driven in parallel by and , as shown in Figures 1 and 2, and a coupled movement is generated. While link 1 is rotating about and , its moving direction along the axis is always parallel to . The base reference frame 0 is defined to be the inertial reference frame. Based on this movement relationship, the homogeneous transformation matrix which refer to the inertial reference frame is directly given aswhere
Let be defined as the tool compensation, which is the value of the tool frame point in axis given as
Based on the D-H reference frame, using the attached reference frames to each link, the four parameters that locate one frame relative to another are defined as follows: is the distance from to along , is the angle from to about , is the distance from to along , and is the angle from to about . These geometric parameters of this manipulator for the links 2, 3, 4, and 5 (tool) shown in Figures 1 and 2 are listed in Table 3.
One vector is represented by in the frame and in the frame , respectively; the relationship between and can be expressed by . The homogeneous transformation matrix is conducted to provide a representation of reference frame relative to reference frame according to the D-H reference frames and the geometric parameters. It is given aswhere .
By successive multiplications of the transformation matrices, the equivalent homogeneous transformation matrix can be obtained. Let a vector be represented as in the end-effector reference frame; then this vector related to the inertial reference frame can be expressed as . Here, relates the inertial reference frame. The matrix is obtained aswhere the vector is the position vector of the origin point relative to the inertial reference frame. The vectors , and are the approach vector, orientation vector, and normal vector of the end-effector reference frame 5 relative to the inertial reference frame, respectively. The rotation matrix of the end-effector reference frame 5 relative to the inertial reference frame can be expressed as . The details are as follows:
Let the reference point of the end-effector be represented by the vector relative to the end-effector reference frame 5; the inertial position coordinates of the end-effector can be obtained bynamely,The results of the inertial position coordinates can be represented by the generalized coordinates through substituting the above related equations into (28).
2.1.2. Inverse Kinematics
The structure of this 5-DOF robot manipulator satisfies the Peiper criterion [42], so that closed-form solutions can be obtained through an algebraic method. The inverse kinematics and its concise expressions are given as follows:where the function is the four-quadrant arctangent function; for example, ; that is, it returns the argument of the complex number with real part and imaginary part .whereSince the available range of the rotation angle is as mentioned above,
2.1.3. Representation of Desired End-Effector Inertial Position Trajectory
Assume that the coordinates of the desired inertial position trajectory of the 5-DOF robot end-effector denoted by are specified in the following two forms (see example in Section 6).
(a) Implicit Formwhere is a given smooth function and the initial and final positions and are given.
The implicit form is used in the formulation of the basic multiobjective optimization problem (see Section 3.1).
(b) Parametric Formwhere is a parameter, , are given smooth functions, and the initial and final positions and are given.
The parametric form is used in the formulation of the finite dimensional multiobjective optimization problem (see Section 3.2).
Given a 5-DOF robot motion trajectory , then the resulting actual inertial position of the end-effector, denoted by , is computed by using (28).
Assume that the actual inertial position of the end-effector matches the desired inertial position trajectory specified in implicit form; see (33). Then, it follows that satisfies the following constraint equation:
2.2. Dynamic Modelling Using the RNEA
Because of the need for computational efficiency, the robotics community has focused on this problem of the robotic manipulator dynamics. For inverse dynamics, the recursive Newton–Euler algorithm (RNEA) proposed by Luh et al. [43] remains the most important algorithm [44]. The RNEA uses a Newton–Euler formulation of the problem, and this formulation is especially amenable to the development of efficient recursive algorithms for dynamics computations.
The RNEA is an algorithm with calculation complexity of (where is the number of degrees of freedom in the mechanism), which is aimed to calculate recursive dynamics of fixed base open-loop robot. The calculation divides into four steps shown as follows, with two steps in each of the two recursions.
Step 1. Calculate the velocity and acceleration of each link in sequence, beginning with the known velocity and acceleration of the fixed base and calculating towards the end link of robot. The velocity of each link in the robot is obtained by the recursive formulawhere is the velocity of link , is the motion matrix of joint , and is the vector of joint velocity variables for joint .
The equivalent formula for acceleration is given by differentiating (36):where is the acceleration of link and is the vector of joint acceleration variables.
Through initializing to instead of zero (where is the gravitational acceleration vector and the gravity is simulated by a fictitious base acceleration), the effect of a uniform gravitational field on robot can be simulated. In this case, is not the true acceleration of link , but the sum of its true acceleration and .
Step 2. Calculate the equation of motion for each link. This step computes the forces required to cause the acceleration calculated in Step . The equation of motion for link is shown aswhere is the spatial inertia of link i and is the net force acting on link i.
Step 3. Calculate the spatial force on each joint. According to Figure 3, the force acting on link is as follows:where is the force transmitted across joint i, is the sum of all relevant external forces acting on link , and is the connected set of link . Rearranging (38) gives the recursive formula for calculating the joint forces:where iterates from to 1.
Step 4. Calculate the joint force variables, . They are given by the equation
In this section, the dynamics problem (i.e., inverse dynamics problem) of this 5-DOF robot manipulator is solved based on this very efficient RNEA. The velocity and acceleration of each link and the resultant forces on each link are first computed through an outward recursion in turn, starting with the known velocity and acceleration of the fixed base and working towards the tips. A second, inward recursion uses the force balance equations at each link to compute the spatial force across each joint and the value of each joint torque/force variable. Note that, without loss of generality, only the case of the rigid robot manipulator is stressed here; other contributing factors in the dynamic description of the robot manipulator (which may include the dynamics of the actuators, joint and link flexibility, friction, noise, and disturbances) are not considered in this paper.
Define that , , , , , and are the angular velocity, the linear velocity, the angular acceleration, the linear acceleration, the force, and the moment of the reference frame with respect to the inertial reference frame and expressed in the reference frame , respectively. Also define that , , and are the linear acceleration of the center of mass of link , the inertial force, and the inertial moment acting on the center of mass of link in the reference frame with respect to the inertial reference frame and expressed in the reference frame , respectively.
The vector is defined as the axis of the reference frame expressed in the reference frame . The vector is defined as the origin of the reference frame with respect to the reference frame and expressed in the reference frame . The rotation matrix is defined as the reference frame with respect to the reference frame and expressed in the reference frame .
The mass (), center of mass (), and moment of inertia about the center of mass () of link expressed in reference frame are listed in Table 4.
2.2.1. Outward Recursion
(a) Link 1. As shown in Figures 1 and 2, link 1 is driven in parallel by and , and a coupled movement is generated. The angular and linear velocity and acceleration of reference frame 1 representing link 1 with respect to the inertial reference frame and expressed in the inertial reference frame can be directly given as follows:where m/s2, according to (1), yielding
The angular and linear velocity and acceleration of the reference frame 1 and the linear acceleration of the center of mass of link 1 and the inertial force and the inertial moment acting on the center of mass of link 1 in the reference frame 1, with respect to the inertial reference frame and expressed in the reference frame 1, are derived aswhere the rotation matrix is given by
(b) Link 2. The angular and linear velocity and acceleration of the reference frame 2 and the linear acceleration of the center of mass of link 2 and the inertial force and the inertial moment acting on the center of mass of link 2 in the reference frame 2, with respect to the inertial reference frame and expressed in the reference frame 2, are derived aswhere
(c) Link 3. The angular and linear velocity and acceleration of the reference frame 3 and the linear acceleration of the center of mass of link 3 and the inertial force and the inertial moment acting on the center of mass of link 3 in the reference frame 3, with respect to the inertial reference frame and expressed in the reference frame 3, are calculated as follows:where
(d) Link 4. The angular and linear velocity and acceleration of the reference frame 4 and the linear acceleration of the center of mass of link 4 and the inertial force and the inertial moment acting on the center of mass of link 4 in the reference frame 4, with respect to the inertial reference frame and expressed in the reference frame 4, are given as follows:where
2.2.2. Inward Recursion
For each link of robot manipulator, except the gravitational force, the forces and moments acting one link belong to internal forces and moments acting between the connected links; that is, the external forces and moments relationships between connected links are forces/moments and reaction forces/moments. A diagram is given for displaying the external forces and resultant forces (shown in Figure 4). External forces and moments acting on link are shown in Figure 4(a). Here, the force and moment of link acting on link are resultant and , respectively; also, the force and moment of link acting on link are resultant and , respectively; is the gravitational force of link (where is the mass of link and is gravitational acceleration vector). The resultant force and moment at the center of mass of link caused by these external forces and moments acting on link based on static balancing conditions are and (shown in Figure 4(b)), respectively.
(a)
(b)
(a) Link 4. When a robot performs a manipulation task (e.g., polishing a surface), the robot establishes physical contact with external environment; the relevant external forces and/or moments will act on the end-effector. In order to express this influence in the dynamic equations of motion, suppose that the resultant of all relevant external forces and moments applied on the origin of the end-effector reference frame 5 is and , respectively. Thus, the force and moment of link 4 applied on the origin of the reference frame 5 are and , respectively. The payload mass can be entered as an external force .
The force and the moment of the reference frame 4 with respect to the inertial reference frame and expressed in the reference frame 4 and the required joint actuator torque of joint 5 are derived aswhere
Note that, for link 4, the inertial force and the inertial moment (which are expressed in (70) and (71), resp.) are represented as the term in (74) in order to apply the influence of the center of mass of link 4 on the origin of the reference frame 4. Other influences of the center of mass of link have the same representation forms.
(b) Link 3. The force and the moment of the reference frame 3 with respect to the inertial reference frame and expressed in the reference frame 3 and the required joint actuator torque of joint 4 are obtained as
(c) Link 2. The force and the moment of the reference frame 2 with respect to the inertial reference frame and expressed in the reference frame 2 and the required joint actuator torque of joint 3 are expressed as follows:
(d) Link 1. The force and the moment of the reference frame 1 with respect to the inertial reference frame and expressed in the reference frame 1 are given as follows:
The interaction forces between link 1 and screw joints (consisting of ball screws) are shown in Figure 5. Combining Newton’s and Euler’s equations for the link 1 leads towhere and .
(a)
(b)
The required joint actuator forces of joint 1 and joint 2 are derived by (85); that is,
Ignoring transmission power loss, the resulting five applied motor torques (while the five reaction torques are equal in magnitude and opposite in direction) can be obtained as
Thus, given a 5-DOF robot trajectory , then the applied forces and torques required to produce the above-mentioned trajectory are computed as follows.
Step 1. Calculate the velocity and acceleration of each link in turn, starting with the known velocity and acceleration of the fixed base, and working towards the end. They are computed by (44)–(47), (52)–(54), (59)–(61), and (66)–(68).
Step 2. Calculate the equation of motion for each link. They are determined by (50), (57), (64), and (71).
Step 3. Calculate the spatial force on each joint. They are given by (74), (78), (81), and (84).
Step 4. Calculate the joint force variables (i.e., the applied forces and torques ). They are obtained by (75), (79), (82), (86), and (87).
2.3. Dynamic Modelling Using the Lagrange Equations
Besides the RNEA, the robot dynamic equations of motion can also be derived by the Lagrange equations. The dynamic equations can be deduced by calculating the kinetic and potential energies of the links and by defining the Lagrangian and by differentiating the Lagrangian equation with respect to the joint variables [45].
2.3.1. Kinetic Energy
As mentioned above, the link 1 is driven in parallel by and , and a coupled movement is generated. In order to separately express this coupled movement, an addition virtual link with zero mass and zero length is introduced between the base and link 1 to produce a translational motion. Thus, a translational motion of the virtual link along the axis and a rotation motion of the link 1 about the axis can be decoupled from this coupled movement. The corresponding transformation matrix can be represented aswhere
Other transformation matrices for the link are represented by , respectively; that is,where the joint variables are given as follows:
The transformation between the reference frame and the inertial reference frame can be written as
Using local generalized coordinates to represent these joint variables ( and for prismatic joints and , and for revolute joints), the derivative of matrix for a revolute joint or a prismatic joint can be expressed aswhere
Extending the same differentiation principle to the matrix of (93) with multiple joint variables differentiated with respect to only one variable gives
Higher-order derivatives can be formulated similarly from
Using to represent a point on any link of the robot relative to frame , the position of the point can expressed by premultiplying the vector with the transformation matrix representing its frame as follows:
Differentiating (98) with respect to all the local variables yields the velocity of the point:
The kinetic energy of an element of mass on a link iswhere the mathematical symbol “Trace” is the trace of a matrix. Integrating this equation and rearranging terms yield the total kinetic energy:
The Pseudo Inertia Matrix, representing the terms, can be written aswhere is the location of the center of mass of the link relative to the frame representing the link. For the local generalized variables , the corresponding Pseudo Inertia Matrices are , , , , and , respectively. The Pseudo Inertia Matrices and the location of the center of mass of the studied robot manipulator are listed in Table 5.
Substituting (102) into (101) gives the final form for kinetic energy of the studied robot manipulator (here, the kinetic energy of the actuators is not stressed in this paper):
2.3.2. Potential Energy
The potential energy of the system is the sum of the potential energies of each link and can be written aswhere is the gravity matrix. Notice that the values in the gravity matrix are dependent on the orientation of the reference frame.
2.3.3. The Lagrangian
The Lagrangian of the robot manipulator is then expressed as
2.3.4. The Lagrangian Dynamics for the Robot’s Equations of Motion
The Lagrangian of the robot manipulator can be differentiated for local generalized coordinates in order to form the dynamic equations of motion:
The final equations of motion for all local generalized variables can be derived as follows:where
The following derivations of transformation relationship, between equations of motion (represented by the generalized coordinates ) and equations of motion (represented by the local generalized variables ), can be yielded as follows.
The relationship between forces and , respectively, produced by the generalized coordinates and and and , respectively, produced by the virtual generalized coordinate and can be derived using the principle of virtual work and (90) as
Based on (6) and (109), it yields
Similarly, based on (92), it gives
From the above derivations in this section, the dynamic equations of motion represented by the local generalized variables and the generalized coordinates for the studied 5-DOF robot manipulator are obtained, respectively.
Silver compared the computational efficiency between Newton–Euler and Lagrangian formulations for manipulator dynamics and proved that they are indeed equivalent to each other [46]. There is in fact no fundamental difference in computational efficiency between two. The former is significantly more efficient than the latter due to two factors: the recursive structure of the computation and the representation chosen for the rotational dynamics. Namely, the former works directly with Newton’s and Euler’s equations for a rigid body, which are contained within the spatial equation of motion. This formulation is especially amenable to the development of efficient recursive algorithms for dynamics computations. The latter only considers the total robot energy and proceeds via the Lagrangian of the robot mechanism; then the dynamic equations of motion are developed using Lagrange’s equation for each defined generalized coordinate.
3. Problem Formulation
For a given prespecified trajectory in the task space of this 5-DOF robot manipulator in this paper, the aim is to optimize a trajectory and calculate the maximum DLCC. This can be formulated as an optimization problem which optimizes four objective functions considering the robot manipulator physical constraints, payload constraint, and actuator limitations.
3.1. Cost Function and Constraints
In this paper, the cost function is the weighted sum of four objective functions representing traveling time, total energy involved in the motion of the robot manipulator, joint jerks, and joint acceleration, respectively. The cost function is considered as the evaluation criterion for the optimal trajectory planning problem.
In practice, the results of the multiobjective optimization problems are influenced by the values of the weighting coefficients. Hence, appropriate values are required to characterize the importance of the corresponding decision criteria for the trajectory planning applications. Multiple objectives can be combined into a scalar cost function via a weight vector. Weights may be assigned through the direct assignment, the eigenvector method, the empty method, the minimal information method, the random determination, or the adaptation determination. In this paper, it is considered that minimizing operation time (minimum traveling time) and saving energy (minimizing total energy involved in the motion) take the highest priority over others. The next priority is applied for reducing vibrations of robot manipulator joints (minimizing joint jerks and acceleration). Thus, a corresponding vector of normalized weights, , , and , is assigned for the objective functions in the optimization problem.
Significant differences may be made in the magnitude of each function when the objective functions are simply weighted and used to produce a single cost function, resulting in a dominating effort that is mainly caused by the function corresponding to the largest magnitude. In particular, the larger the magnitude of an objective function, the higher the influence on the cost function. Hence, a feasible and efficient means will be needed to deal with this problem effected by the large magnitude of the objective function. One possibility of avoiding this problem is to make an adequate normalization of the objective functions and formulate the same magnitude of objective functions. It can have a roughly equivalent influence on the objective functions for the cost function. Also normalizing parameters make all objective functions as unitless functions. In this paper, the normalized weighted objective functions method is used to balance the influence on four objective functions and select the best found solution. More details on this method can be referred to [47, 48]. The values of , , , and are taken as the normalizing parameters of the objective functions , , , and .
As mentioned above, joint variables of this robot manipulator form a vector of generalized coordinates. Likewise, the variables , , , and are generalized velocities, acceleration, jerks, and forces.
Given the following date, (a) the initial time ; the final time which is unknown and must be determined (see below); (b) the coordinates of the desired inertial position trajectory of the 5-DOF robot end-effector specified in two forms (see Section 2.1.3), the resulting desired initial and final conditions of the generalized coordinates and are calculated by using (29)–(32); set and ; (c) the initial and final conditions are , , , , , and ; thus, the basic multiobjective optimization problem is stated as follows.
Compute the final time , the payload mass , and a smooth 5-DOF robot motion trajectory satisfying the above-mentioned initial and final conditions so as to minimize the cost function given byWhere one has the following: Total traveling time between initial and final configurations: Total energy involved in the motion of the robot manipulator: Integral of joint jerks squared: Integral of joint acceleration squared:subject to the following constraints: The final time constraint: The displacement constraint: The velocity constraint: The acceleration constraint: The jerk constraint: The force/torque constraint:that is, . The payload constraint: The kinematic constraints: that is, (29)–(32) (which are specific specified in Section 2.1). The dynamic constraints: that is, (75), (79), (82), and (86)–(88) (which are specific specified in Section 2.2 are also considered for the multiobjective optimization problem).
Here, and are the jth joint lower and upper bounds of displacement, respectively. , , , and are the jth joint upper bounds of velocity, acceleration, jerk, and force/torque, respectively. and are the minimum and maximum allowable load corresponding to the DLCC of the robot manipulator, respectively.
3.2. Finite Dimensional 5-DOF Robot Motion Trajectory Representation
The cubic spline in trajectory planning is very useful, because the generated trajectories have continuous values of the acceleration; also, unlike higher-order polynomials, cubic splines would not generate problems such as excessive oscillations and overshoot between any pair of reference points.
In this paper, cubic splines are applied to represent the trajectory of the joint space, given that the desired inertial position trajectory of the 5-DOF robot end-effector is specified in parametric form (see Section 2.1.3). Also, given a positive integer , select an increasing sequence of parameter values lying in the parameter interval such that
Thus, via-points points are obtained on the desired inertial position trajectory of the 5-DOF robot end-effector with coordinates
By using the inverse kinematics of the 5-DOF robot (see (29)–(32)), via-points are obtained on the desired trajectory of the 5-DOF robot generalized coordinates at a sequence of increasing time instants indexed as follows:where the initial time is given and .
An additional two via-points and , corresponding to times and , are added in order to accommodate the specified initial and final conditions , , , and (see below).
Thus, a total of via-points on the desired trajectory of the 5-DOF robot generalized coordinates are known and are denoted as follows:
What is not known and must be determined is the sequence of increasing time instants,
The sequence of time instants , defines consecutive time intervals given byLet denote the length of the time interval ; that is, .
Given the above-mentioned via-points and the time-interval lengths , then the actual motion trajectory of the 5-DOF robot is constructed by using a cubic spline segment to connect every pair of consecutive via-points; see (127).
In this section, a methodology is proposed for computing the unknown time-interval lengths , such that the resulting actual motion trajectory of the 5-DOF robot minimizes a cost function subject to constraints (see later).
Let be the cubic polynomial for the jth joint defined on the interval . Solving the interpolation problem means to find cubic polynomials satisfying the imposed conditions. In order to avoid considering cumbersome calculations with unknowns, it can be adopted that the second derivative of is a linear function of in the interval ; namely,
Integrating (130) and setting and , the following interpolating functions are obtained:
The joint displacements of these two extra via-points are given as
Since the derivative of is continuous at , the inner via-points can be expressed as follows:
Combining and (134) giveswhere
Then, a linear system with unknowns (the acceleration at the inner via-points) is obtained aswhere and are -dimensional vector and the coefficient matrix (the same for all joints) is a nonsingular and band-diagonal matrix as follows:
Given the acceleration and velocities for the first and the last via-point, this linear system with unknowns can be obtained. Thus, the acceleration at the inner via-points, namely,with , can be obtained by solving the linear system (137).
Thus, the actual motion trajectory of the 5-DOF robot is constructed by joining together the above-mentioned spline functions; that is, if then , , , , ,,, and (which is a new symbol for force defined for ).
In order to solve the trajectory planning problem, considering that the trajectory is formed by cubic splines, an explicit expression for objective functions and constraints needs to be formulated.
As for the displacement, the constraints are denoted by
As for the velocity, the constraints are given by
Here we note that these constraints are of the semi-infinite type; that is, they must hold for any value of the continuous variable in the indicated time interval. It would be much easier to deal with them if they can be transformed into a finite form. This is feasible noticing that the velocity is parabolic in the interval ; hence it reaches its maximum value either at one of the interval ends or at a certain intermediate instant when the acceleration becomes zero. This intermediate instant can be expressed asThe value of the velocity at the instant is obtained asThe values of the velocity at the instant ends are given byHence, the semi-infinite constraints (26) can be expressed much more accurately aswhich represent a finite number of relationships that must be satisfied for the planned trajectory.
The kinematic constraints on acceleration and jerk can be similarly reduced to a finite form. We note that the acceleration is linear in the interval and therefore it must reach its maximum value at one of two interval ends. Hence, the semi-infinite constraints on the acceleration,can be expressed in a finite form as
The jerk can be obtained as follows:
As for the jerk, the semi-infinite constraint is given bySince the jerk is constant on any interval , the finite form of the jerk yields
As for the force/torque, the semi-infinite constraint givesthat is, , , , , and , .
Moreover, in the optimization of the trajectory, the time intervals which are the true optimization parameters (i.e., the decision variables) must be particularly considered. Since in any interval the condition must be satisfied. Hence, the lower bound of the intervals must be subjected to
Thus, the finite dimensional multiobjective optimization problem is formulated as follows. Compute so as to minimize the cost function :wheresubject to the constraints (which are (140), (141), (147), (150)–(152), (29)–(32), (75), (79), (82), (86), and (87)).
The vector of decision variables is defined by the sequence of the time intervals associated with each two consecutive via-points; thus . When the intervals and joint displacements (i.e., these two extra via-points and joint displacements obtained by solving the inverse kinematics of the 5-DOF robot manipulator) are determined, other functions , , , , and can also be determined by these time intervals and joint displacements. The initial and final velocity ( and ) and acceleration ( and ) for all joints are set to zero corresponding to the initial and final time instants. The payload mass is fixed to a feasible value satisfying (123).
A genetic optimization algorithm NSGA-II is applied in order to try and solve the above-mentioned optimization problem. The algorithm NSGA-II runs for a number of iterations and then stops after meeting some termination criterion, thus producing a vector of decision variables in solving trajectory planning and DLCC of the 5-DOF robot. These obtained decision variables and other above-mentioned functions can be represented as the best found solution that have been able to find so far.
Actually, the real advantages of NSGA-II used in this application compared to the conventional constrained optimization algorithms (e.g., sequential quadratic programming) are that the former is population-based approach search performing the genetic operations like crossover and mutation for the elements of the decision variables (where the unfit individuals are replaced by the fit individuals through the crossover and mutation operator in population), thus discovering the global optimal solution is possible; also, the former is no need to provide any auxiliary relationships (e.g., gradients and Hessian matrix); moreover, the latter is easy to fall into local minima as the susceptibility to initial values, and it needs to calculate the gradients and the Hessian matrix of the objective functions and constraints resulting in more computational time. In Section 6, a classical nonlinear constrained optimization function (i.e., the function fmincon in Matlab Optimization Toolbox) based on sequential quadratic programming algorithm is also applied to solve the above-mentioned optimization problem and compared with the best found solution of the genetic optimization algorithm NSGA-II.
4. The Proposed Approach for Trajectory Planning
In this section, an approach is proposed to deal with the multiobjective optimization problem for a given prespecified trajectory taking into account a feasible DLCC and obtain time intervals corresponding to the via-points, joint displacement, velocity, acceleration, jerk, and force/torque, considering the objective functions and two categories of constraints (inequality linear constraints and equality nonlinear constraints). This feasible DLCC is marked as and must fall within the DLCC range of reachable values. The detailed calculation for maximum DLCC will be presented in the next section.
Note that the cost function results would converge after a certain number of generations by using the calculation mechanism of the NSGA-II algorithm. In this paper, this stabilized condition is made as the termination criterion, and the corresponding calculation results are taken as the best found solution for the trajectory planning problem and maximum DLCC calculation of this 5-DOF robot manipulator.
The desired joint displacement (i.e., the via-points in the joint space) can be determined by solving the robot manipulator inverse kinematics problem. When the optimization process is performed, the equality nonlinear constraints and inequality linear constraints described in the previous section are considered in the objective functions. The suitable population size (pop), total number of generations (gen), crossover probability (pc), mutation probability (pm), crossover distribution index (dc), and mutation distribution index (dm) in NSGA-II algorithm are needed to select so that the algorithm leads the robot manipulator to have a smooth trajectory and minimum chattering in profiles of velocity, acceleration, jerk, and force/torque.
The detailed description of the proposed approach for trajectory planning is presented in the following steps.
Step 1. Calculate desired joint displacement (i.e., the via-points in the joint space) according to the given desired trajectory in the task space by solving the robot manipulator inverse kinematics problem using (29)–(32).
Step 2. Select , pop, gen, pc, pm, dc, and dm, set the initial and final velocity ( and ) and acceleration ( and ) for all joints , set the constraints of displacement, velocity, acceleration, jerk, and force/torque, and calculate the lower bound of the time intervals using (152).
Step 3. A set of small time-interval lengths (except the time intervals of two extra via-points), are randomly generated using (152) and the constraints , , , , , and need to be satisfied, then, respectively, interpolating randomly at the first and last time-interval lengths (i.e., and ) for the time intervals of two extra via-points generates new time-interval lengths , , , and , and finally assembling , , , , , and constructs the initialized time-interval lengths ; generate two extra via-points using (133), and arrange the displacement in the NSGA-II algorithm.
Step 4. If lower bound constraints of the time intervals are satisfied (which are judged by (152)), go back to Step . Otherwise, go back to Step .
Step 5. Calculate , , , and using (131), (137), and (148) and the dynamic equations (75), (79), (82), (86), and (87), respectively. If the constraints (displacement, velocity, acceleration, jerk, and force/torque) are satisfied (which are judged by (140), (145), (147), (150), and (151), resp.), go back to Step . Otherwise, go back to Step .
Step 6. Optimize the time intervals (i.e., the decision variables) and objective functions (which are formulated as (154)) by using the calculation mechanism of the NSGA-II algorithm, compare the corresponding cost function result (i.e., (153)) of each population for the current generation, and find out the minimum one (, where i denotes the current th generation and is the number of the total population) and the corresponding best found calculation results (decision variables, i.e., the time intervals, and the corresponding results of objective functions, displacement, velocity, acceleration, jerk, and force/torque). Then, save and report these selected best found calculation results for each generation. This step continues until all generations are completed.
Step 7. Compare the cost function results of each generation, and find out the minimum one (, where is the total number of generations) and the corresponding other calculation results for the whole calculation process.
Through this detailed calculation process, the required calculation decision variables, that is, the time intervals, and the corresponding best found calculation values (objective functions, cost function, displacement, velocity, acceleration, jerk, and force/torque), taking equality nonlinear constraints and inequality linear constraints into account, are obtained by using the NSGA-II algorithm.
5. The Proposed Approach for DLCC Calculation
In many practical applications, the maximum DLCC must be a feasible one to satisfy certain performance requirements of the robot manipulator. It is appropriate to specify a robot manipulator in terms of useful payload as a function of performance rather than just in terms of maximum capacity. Therefore, a feasible maximum DLCC of the robot manipulator, from a safety and reliability point of view, plays a significant role in practical applications. The dominant constraint that limits the DLCC of a robot manipulator is the limitations on the applied forces/torques in the robot manipulator joints. A feasible maximum DLCC, considering this dominant constraint (i.e., the applied forces/torques), is addressed in this section.
Another approach is proposed to calculate the maximum DLCC for a given prespecified desired trajectory which is described in the task space considering the objective functions and two categories of constraints as mentioned above. Firstly, a feasible and suitable DLCC is introduced for the trajectory planning based on the proposed approach in the above section by solving the multiobjective optimization problem for this given trajectory, and the related best found results (its best found time intervals, joint displacement, velocity, acceleration, and jerk) can be obtained. Then these determined results are taken as the initial conditions to calculate the maximum DLCC and also the motor torque requirements using the iterative method. The process of the iterative solution is repeated until the termination condition (i.e., the motor torque requirements are beyond the rated torques) is reached. Here, at each iteration, the DLCC, that is, the allowable payload mass that is carried by the robot manipulator, is increased by an iteration procedure. and are the initial value and the increment of at each iteration, respectively.
The proposed approach for the maximum DLCC calculation proceeds in nine steps as follows.
Step 1. Calculate desired joint displacement (i.e., the via-points in the joint space) according to the given trajectory in the task space by solving the robot manipulator inverse kinematics problem using (29)–(32).
Step 2. Select , , pop, gen, pc, pm, dc, and dm, set , set the initial and final velocity ( and ) and acceleration ( and ) for all joints , set the constraints of displacement, velocity, acceleration, jerk, and force/torque, and calculate the lower bound of the time intervals using (152).
Step 3. A set of small time-interval lengths (except the time intervals of two extra via-points), are randomly generated using (152) and the constraints , , , , , and need to be satisfied, then, respectively, interpolating randomly at the first and last time-interval lengths (i.e., and ) for the time intervals of two extra via-points generates new time-interval lengths , , , and , and finally assembling , , , , , and constructs the initialized time-interval lengths ; generate two extra via-points using (133), and arrange the displacement in the NSGA-II algorithm.
Step 4. If lower bound constraints of the time intervals are satisfied (which are judged by (152)), go back to Step . Otherwise, go back to Step .
Step 5. Calculate , , , and using (131), (137), (148) and the dynamic equations (75), (79), (82), (86), and (87), respectively. If the constraints (displacement, velocity, acceleration, jerk, and force/torque) are satisfied (which are judged by (140), (145), (147), (150) and (151), resp.), go back to Step . Otherwise, go back to Step .
Step 6. Optimize the time intervals (i.e., the decision variables) and objective functions (which are formulated as (154)) by using the calculation mechanism of the NSGA-II algorithm, compare the cost function results of each generation, and find out the minimum one (, where is the total number of generations) and the corresponding other calculation results (decision variables, i.e., the time intervals, and the corresponding results of objective functions, displacement, velocity, acceleration, jerk, and force/torque). Then, save and report these selected best found results.
Step 7. Set ; increase payload by .
Step 8. Calculate the required motor torques using the dynamic equations (75), (79), (82), (86), and (87) according to the selected best found results (i.e., the best found displacement , velocity , and acceleration ) and the current payload mass . If the constraints of the rated motor torques are satisfied (which are judged by (151)), go back to Step . Otherwise, go back to Step .
Step 9. Save and report the latest DLCC (i.e., the maximum DLCC).
Through this detailed calculation steps while taking into account equality nonlinear constraints and inequality linear constraints, according to the best found initial conditions (i.e., the selected best found results of displacement , velocity , and acceleration ), the maximum DLCC for the given trajectory are determined based on the NSGA-II algorithm and the iterative method.
6. Numerical Example
Based on the above proposed approaches, the trajectory planning and the maximum DLCC calculation of this robot manipulator for a given desired end-effector inertial position trajectory are presented in this section.
A saddle line segment is selected as the given desired inertial position trajectory of the 5-DOF robot manipulator end-effector for the trajectory planning and the maximum DLCC calculation. The saddle line segment is defined between initial coordinates (mm, mm, and mm) and final coordinates (mm, mm, and mm). In this paper, the equation of this saddle line segment is expressed in implicit form as follows:
The saddle line segment is also expressed in parametric form as follows:where the parameter . In this case via-points are chosen on the desired end-effector trajectory leading to the following values for the parameter : .
This saddle line segment (i.e., the prespecified desired trajectory) and each corresponding via-point and orientation in the task space are shown in Figure 6, where the orientation is determined as follows: the vector is defined along the normal vector of the via-point (cyan arrows as shown in Figure 6), is defined along the tangent vector of the via-pointand points the next via-point (blue arrows as shown in Figure 6), and is defined by . The generalized coordinates of this 5-DOF robot manipulator can be represented as . The robot manipulator is at rest initially and comes to a full stop at the end of the trajectory (i.e., , , , and for all joints ). In this paper, joint force/torque of this robot manipulator is generated by AC servomotors. The feasible DLCC (i.e., ) is taken as a fixed payload mass kg for solving the trajectory planning problem.
The via-points (i.e., the desired joint displacement) including two extra via-points in the joint space are determined by solving the inverse kinematics of the robot manipulator and given in Table 6. The real variable type is considered in the NSGA-II algorithm. The values of the parameter which is used in algorithm, population size (pop), total number of generations (gen), crossover probability (pc), mutation probability (pm), crossover distribution index (dc), and mutation distribution index (dm), are 100, 200, 0.9, 0.1, 10, and 100, respectively.
The proposed approaches are applied both to trajectory planning problem and maximum DLCC calculation for this 5-DOF robot manipulator, and the related best found solutions are obtained. Note that, without loss of generality, only the case of the rigid robot manipulator is stressed here; other contributions to the dynamic description of the robot manipulator (which may include the dynamics of the actuators, joint and link flexibility, friction, noise, and disturbances) are not considered in this paper. The relationship between the cost function and the total number of generations is shown in Figure 7. It can be found that the cost function value reached the convergence to a stabilized condition from 66th to 200th generation. The relationships between objective functions (total traveling time (), total energy involved in the motion (), integral of joint jerks squared (), and integral of joint acceleration squared ()) and the corresponding total number of generations (pop) are shown in Figures 8(a), 8(b), 8(c), and 8(d), respectively.
(a)
(b)
(c)
(d)
The best found variables (i.e., time intervals) obtained from the proposed approach for trajectory planning are given in Table 7. The best found joint displacement is given in Table 8. The best found joint displacement, velocity, acceleration, jerk, and force/torque are shown in Figures 9, 10, 11, 12, and 13, respectively. The actual end-effector inertial position is obtained as shown in Figure 14. The end-effector inertial position tracking error between the actual inertial motion trajectory and the desired inertial position trajectory is shown in Figure 15. The standard deviation tracking error, the variance tracking error, and the maximum tracking error are 0.0840 mm, 0.0071 mm, and 0.2933 mm (the last one shows the nearby point of the first extra point), respectively.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
The best found joint displacement, velocity, and acceleration are taken as the initial conditions (the same as shown in Figures 9, 10, and 11, resp.) to calculate the maximum DLCC and also the motor torque requirements (i.e., requirements of joint force/torque). Note that, without loss of generality, only the case of the rigid robot manipulator is stressed here; other contributions to the dynamic description of the robot manipulator (which may include the dynamics of the actuators, joint and link flexibility, friction, noise, and disturbances) are not considered in this paper. The proposed DLCC calculation approach is applied and the maximum DLCC, that is, the maximum allowable load, is found to be 12.65 kg.
In this paper, the cost function value converges to a stabilized status from 66th to 200th generation, so produced vector of decision variables (i.e., the time-interval lengths ) and other above-mentioned functions of the last generation are used to represent the solution of trajectory planning and DLCC of the 5-DOF robot. Strictly speaking, it cannot be rigorously proved that this found solution is the constrained global minimizer of the cost function and it only represents the best found solution that have been able to find so far.
The Matlab function of the proposed approach for trajectory planning implementing the optimization algorithm NSGA-II and an example multimodal cost function subject to some nonlinear constraints showing how to call/use the function implementing the NSGA-II algorithm are given in Appendix and Appendix (both in Supplementary Materials available online at http://dx.doi.org/10.1155/2016/1302537), respectively.
A classical nonlinear programming solver, that is, fmincon in Matlab Optimization Toolbox (which is based on sequential quadratic programming algorithm), is applied to solve the above-mentioned minimum of constrained nonlinear objective functions optimization problem and compared with the best found solution of the genetic optimization algorithm NSGA-II.
Initializing randomly the time-interval lengths , just like the initialized way of Step in Section 4 generates the initial value of the fmincon function. The same normalized weights () and normalizing parameters () as mentioned in Section 3 are used in fmincon solver. The same processes of initializing way and solving are performed 100 times and generate 100 solving results with 1032.78 seconds. The mean value of calculation results is listed in Table 9. The best found solution of the algorithm NSGA-II for the above-mentioned optimization problem (which is computed with is 74.55 seconds) is also listed in Table 9. The results of the cost function and the objective functions in fmincon function solver are given in Appendix in Supplementary Materials.
Comparing these results, it can be observed that the found solution of the genetic optimization algorithm NSGA-II is better than the fmincon function solver with respect to the cost function (), total traveling time (), integral of joint jerks squared (), and integral of joint acceleration squared (). This is because the former is population-based approach which can perform the genetic crossover and mutation operations for the elements of the decision variables and desert the unfit individuals in each generation. Also, these genetic operations can extend and discover the possibly global solution. The computational time of the former is less than the latter, mainly because the algorithm NSGA-II does not need to provide any auxiliary relationships such as the gradients and Hessian matrix of the objective functions and constraints.
7. Conclusions
A new methodology using a direct method is presented for obtaining the best found trajectory planning and maximum dynamic load-carrying capacity (DLCC) for a given prespecified trajectory based on an evolutionary algorithm (i.e., NSGA-II). A cubic spline curve is used to represent the trajectory. This problem is transcribed into a nonlinear constrained multiobjective optimization problem. To summarize, the major contributions of this paper are as follows:(i)The nonlinear constrained multiobjective optimization problem is formulated with four objective functions (i.e., total traveling time, total energy involved in the motion, joint jerks, and joint acceleration) and a cost function (which is a weighted sum of these objective functions). These decision criteria can ensure the best found trajectory is smooth and lead to minimum actuator power and energy.(ii)Two separate approaches are proposed to deal with trajectory planning problem and maximum DLCC calculation for a 5-DOF robot manipulator, respectively. These approaches are easy to program and implement and can also be extended to handle similar problems for other types of the robots.(iii)A numerical example verified the effectiveness of the two proposed approaches and obtained best found solutions for trajectory planning and maximum DLCC calculation of a 5-DOF robot manipulator.
In order to get more accurate and smoother trajectory, the nonuniform rational B-spline (NURBS) curve with more via-points will be used to represent the trajectory in future work. The presented new methodology based on the evolutionary optimization technique and two proposed approaches are expected to provide some insight into the foundational aspects of trajectory planning and maximum DLCC calculation.
Competing Interests
The authors declare that there are no competing interests regarding the publication of this paper.
Acknowledgments
This work was supported by Self-Planned Task (no. SKLRS201609B) of State Key Laboratory of Robotics and System (HIT).
Supplementary Materials
Supplementary Materials provided the Matlab function of the proposed approach for trajectory planning implementing the optimization algorithm NSGA-II, an example multimodal cost function subject to some nonlinear constraints showing how to call/use the function implementing the optimization algorithm NSGA-II, and the results of the cost function and the objective functions in fmincon function solver.