Motion Planning of Nonholonomic Mobile Manipulators with Manipulability Maximization Considering Joints Physical Constraints and Self-Collision Avoidance

: The motion of nonholonomic mobile manipulators (NMMs) is inherently constrained by joint limits, joint velocity limits, self-collisions and singularities. Most motion planning algorithms consider some of the aforementioned constraints, however, a uniﬁed framework to deal with all of them is lacking. This paper proposes a motion planning solution for the kinematic trajectory tracking of redundant NMMs that include all the constraints needed for practical implementation, which improves the manipulability of both the entire system and the manipulator to prevent singularities. Experiments using a 10-DOF NMM demonstrate the good performance of the proposed method for executing 6-DOF trajectories while satisfying all the constraints and simultaneously maximizing manipulability.


Introduction
A mobile manipulator is a robotic system that consists of a standard robot manipulator mounted on a mobile platform. This system integrates the dexterity provided by the manipulator with the extended workspace provided by the platform. Additionally, the combination of both subsystems usually introduces kinematic redundancy, which increases flexibility and dexterity. Therefore, mobile manipulators are suitable to perform delicate tasks over a large space, such as welding large parts or painting large, curvy surfaces.
The practical applications of robotic systems commonly define tasks by either a pointto-point movement or a continuous path of the end-effector in task space (also known as operational space). This paper aims to solve the latter, and in particular, focuses on the task space trajectory tracking problem. A trajectory is a path on which a timing law is specified, for instance in terms of velocities and/or accelerations [1]. In other words, not only is the end-effector's pose profile defined, but so is its velocity profile. To accomplish this, a motion planning algorithm that exploits the capabilities of both the manipulator and the mobile platform and that coordinates their movements is required. The redundancy of mobile manipulators can be used to perform additional subtasks or satisfy system constraints. These constraints include joint limits, joint velocity limits, joint velocity boundary constraints (i.e., the constraints on the initial and final joint velocities), and self-collision avoidance. Furthermore, for task space trajectory tracking to be achievable, it is important that the system is kept away from singularities. All these requirements make the motion planning for trajectory tracking a challenging problem.
There exist recent efforts in solving the motion planning of mobile manipulators [2,3]. Liao et al. [2] presented an optimization-based solution that not only handles constraints at the position level, but can also set a target joint configuration for the manipulator at intrinsically improves the manipulability of the robot arm as well as the manipulability of the whole system.
The solution to the motion planning of NMMs for trajectory tracking presented in this work includes joint constraints (range and maximum velocities), self-collision avoidance and manipulation capability preservation. Figure 1 summarizes how all these constraints are included in our solution. Both the particular and homogeneous solutions of our proposed scheme are weighted to avoid joint limits and self-collisions while the trajectory is tracked. The homogeneous solution is used to maximize the manipulability by exploiting the redundancy of the system. To satisfy joint velocity boundary constraints, the step size for searching the maximum manipulability is varied. The joint velocity limits are satisfied by restraining the maximum step size based on the allowable self-motion. Experimental results demonstrate that our method can successfully solve the motion planning problem of NMMs under all the mentioned constraints. This work focuses on task space trajectory tracking at kinematic level. In other words, the outputs of the motion planning algorithm are the joint positions and velocities that will be fed to a joint space dynamic controller for motion control. Then, the motion controller is responsible for suppressing the model uncertainties and external disturbances to guarantee that the actual joint positions follow the ones output by the motion planning algorithm. In the experiments shown in this paper, we use the built-in motion controller of the commercial manipulator and leave the design of our own motion controller for future research.
The contributions of this work are detailed as follows: • A motion planning solution for NMMs that allows to include joint physical constraints and the execution of a secondary task is presented. • Multiple constraints required for the practical implementation of task space trajectory tracking are included in a unified solution. These constraints include joint limits, joint velocity limits, joint velocity boundary constraints and self-collision avoidance. • A new manipulability measure for mobile manipulators is presented. It is demonstrated that the maximization of this measure simultaneously improves the manipulabilities of the whole system and the robot arm.
This paper is organized as follows. Section 2 describes the kinematic modeling of NMMs. Section 3 describes the motion planning problem for trajectory tracking and presents the proposed solution. In Section 4, the concepts that are employed to satisfy each of the mentioned constraints are described and included in the motion planning algorithm.
Experimental results using 6-DOF tasks are presented in Section 5 to validate the efficacy of our approach. Finally, Section 6 concludes the paper.

Kinematic Modeling
The kinematic modeling described here follows the procedure of De Luca et al. [7] and Bayle et al. [6]. For a general procedure of kinematic modeling of wheeled mobile manipulators, please refer to Bayle et al. [5].
Let r ∈ R m be the end-effector's position and/or pose in the task space. The configuration vector q, also known as the generalized coordinates of the mobile manipulator, is given by the combination of the platform configuration vector q p and the robot arm configuration vector q a . Figure 2 illustrates these configuration vectors. A frame x y is attached to the mobile platform at the center of the wheels' axle (x p , y p ), with respect to the world reference frame xy, with its x axis pointing in the forward direction and the y axis pointing in the direction parallel to the wheels' axle. The angle between the x axis of the world reference frame and x attached to the platform is denoted as θ p . Then, the platform configuration is given by q p = x p y p θ p T ∈ R 3 . The robot arm configuration is given by the position of its joints as q a = q a 1 q a 2 . . . q n a T ∈ R n a . Finally, the configuration vector of the mobile manipulator is q = q p T q a T T ∈ R n with n = 3 + n a . The end-effector's position/pose is a function of the configuration vector defined by the kinematic map: The wheeled platform movement is constrained under the rolling without a slipping assumption on both wheels, which can be expressed as the following nonholonomic constraint:q where u p = v p ω p T ∈ R 2 is the velocity input of the platform, consisting of the linear and angular velocities of the platform, which are known as pseudo velocities or quasivelocities. The columns of matrix G(q p ) span the admissible velocity space for a given platform configuration. The matrix G(q p ) for a differential-drive platform is defined as On the other hand, the robotic arm kinematics at the velocity level is not constrained for any configuration, and therefore, the generalized velocities of the arm are equal to the velocity inputs of the arm:q a = u a .
The velocity input vector for the entire NMM can be constructed as u = u p T u a T T ∈ R δ m , with δ m = 2 + n a . The dimension δ m is called the mobility degree of the mobile manipulator and indicates the space dimension of the admissible generalized velocities [6]. Using (2) and (3), the map from the velocity input vector u to the generalized velocitieṡ Tq a T T can be written asq with: where G(q p ) maps the platform's velocity input vector u p = v p ω p T to the platform's generalized velocitiesq p = ẋ pẏpθp T , and matrix I n a is an identity matrix of size n a that sets the arm's generalized velocities equal to its input velocitiesq a = u a . The differential kinematics of the NMM is obtained by differentiating the relation (1) with respect to time:ṙ where the matrices J p ∈ R m×3 and J a ∈ R m×n a are the Jacobians of the platform and the arm, respectively. The matrix J ∈ R m×n is the Jacobian of the mobile manipulator, andJ ∈ R m×δ m is a reduced version of the Jacobian in which the admissible velocities of the platform under the nonholonomic constraint have been included. Equation (5) follows the same form as the differential kinematics of standard manipulators; therefore, the well-known methodologies for the motion planning of redundant manipulators can be extended to NMMs in terms ofJ, including concepts for joint limits avoidance, self-collisions avoidance and manipulability maximization.

Motion Planning Method
The motion planning problem for task space trajectory tracking consists of finding the input velocities u, given the desired end-effector's position/pose r d (t) ∈ R m for t ∈ [t 0 , t f ], such that r(t) follows r d (t). If the task space velocity of the end-effector is set aṡ where K ∈ R m×m is a positive definite matrix, then the convergence of r(t) to r d (t) is guaranteed. Consequently, the motion planning problem is turned into solving the input velocities u from the underdetermined linear equations (5) whereṙ is set according to (6). During the execution of the trajectory, the movement of the joints that get closer to a position constraint (joint limit or self-collision) must be penalized (slowed down). This can be achieved by using a weighting matrix. In this work, we define the weighted input velocity and the reduced weighted Jacobian as follows: where W(q) ∈ R δ m ×δ m is a configuration-dependent positive semidefinite matrix. It is constructed so that its elements penalize the motion of the joints based upon the system constraints. We choose W(q) as a diagonal matrix in this paper; furthermore, we drop the notational dependency of W on q to simplify the expression in the subsequent derivation. Using (7), the system (5) can be rewritten aṡ Because of the kinematic redundancy, i.e., m < δ m , infinite solutions to u W exist. The solution to (8) can be decomposed as a sum of a particular solution and a non-zero homogeneous solution. The particular solution satisfies the end-effector's task, while the homogeneous solution changes the manipulator configuration without changing the end-effector's position/pose.
One way to solve the redundant system (8) is to formulate the problem as a constrained optimization problem [1], where the goal is to minimize a cost function and satisfy the constraint (8). We propose the minimization of the quadratic cost function: this choice is aimed to find a vector of velocities u W that is as close as possible to W 1 2 u 0 , where u 0 ∈ R δ m is a velocity vector that is used to satisfy an additional task such as maximizing the manipulability. The choice of u 0 will be presented shortly. Notice that the physical constraints are also imposed to u 0 by weighting it similarly to howJ W is obtained in (7). By using the method of Lagrange multipliers to minimize g(u W ) with the equality constraint (8), the following solution is obtained: whereJ † W is the Moore-Penrose pseudoinverse of the matrixJ W such thatJ WJ † W = I, I −J † WJ W is the orthogonal projection into the null space ofJ W . The first term of (9) is the particular solution, and the second term is the homogeneous solution. It is trivial to show that I −J † WJ W projects u 0 onto the null space ofJ W by multiplying both sides of (9) byJ W . The velocity input vector is then recovered using the second part of (7): The degree of kinematic redundancy at the velocity level is relevant for this solution to be feasible, since it defines the number of simultaneous constraints that can be satisfied in the differential kinematics inversion withoutJ W losing its rank. The degree of redundancy for NMMs is calculated as R = δ m − m. Whilst analyzing the matrix W 1 2 for three different cases, it is possible to understand the expected behavior of solution (10). Let z represent the number of elements that are zero in the diagonal of W 1 2 , i.e., the number of joints that are forced to stop due to z simultaneously activated constraints. When z < R, both the particular and homogeneous solutions exist, and therefore, the secondary task will still be considered. When z = R, the system is not redundant anymore and only the particular solution exists. Finally, when z > R, the system (8) has no solution. Therefore, for the solution (10) to exist, the condition z ≤ R is necessary.
Our proposed solution (10) has an advantage over the projected gradient solution used in [5,6], because it includes the physical constraints and penalizes both the particular and homogeneous solutions. Furthermore, in contrast with the weighted least-norm method [11], this solution takes advantage of the redundant joints that have not been penalized, due to physical constraints, in attempt to satisfy the task defined by u 0 .
In this paper, the matrix W is constructed to satisfy joint position constraints (joint limits and self-collision avoidance), which is discussed in Section 4.2. The vector u 0 is designated to locally maximize a proposed new manipulability measure for mobile manipulators, as discussed in Section 4.1. To satisfy the joint velocity constraints, the particular and homogeneous solutions are analyzed. Because the particular solution is in charge of following the end-effector's task, it cannot be modified. On the other hand, the homogeneous solution can be varied to satisfy the joint velocity constraints. Furthermore, the homogeneous solution must also consider the admissible velocity commands with respect to the nonholonomic constraints [7]. Taking these two points into account, we define vector u 0 as u 0 = ±α(t)β(t)S T (q)∇ q F(q), (11) where α(t) is a scalar coefficient that is adjusted online to satisfy joint velocity limits (Section 4.3.2), β(t) is a positive variation factor used to satisfy joint velocity boundary constraints (Section 4.3.1), S(q) is defined in (4), and ∇ q F(q) is the gradient of the objective function F(q) : R 3+n a → R, which is set to a manipulability measure. The product α(t)β(t) is the step size of the gradient step ascent/descent, and the sign of u 0 defines whether the objective function F(q) will be maximized (+) or minimized (−).
Replacing (11) in (10), the final form of the proposed solution for u is: where the positive sign for the gradient step descent has been used because we are interested in maximizing the manipulability. Notice that when the sign of α(t) is positive, the objective function F(q) is maximized. However, for some cases, α(t) could be negative for short periods of time to respect joint velocity limits, as explained in Section 4.3.

Manipulability Maximization and Constraints Satisfaction
In this section, the details for manipulability maximization are described, the weighting matrix W is defined so that it penalizes the joint movement to avoid both joint limits and self-collisions, and a scheme to satisfy joint velocity limits as well as joint velocity boundary constraints is presented.

Manipulability Maximization
The term manipulability, introduced by Yoshikawa [12], describes the ability of a robotic system to provide end-effector's velocities in any direction for a given configuration. The manipulability of wheeled mobile manipulators was studied in detail in [5]. There exist different algebraic measures to characterize the manipulability of a robotic system [5,[12][13][14]. The most widely used measure, known as the manipulability measure and noted here as Ω, is given by Ω = σ 1 σ 2 . . . σ m , where σ 1 , σ 2 , . . . , σ m are the singular values of the system Jacobian J(q). Therefore, Ω is the manipulability measure for system configuration q. The measure Ω is proportional to the volume of the manipulability ellipsoid [12]. Furthermore, it can be shown that Ω has the following property: By definition, Ω ≥ 0, and Ω = 0 if and only if rank(J(q)) < m. The elements of ∇ q Ω are given mathematically by As mentioned before, we are not only concerned with improving the manipulability of the whole system but also that of the arm alone. In this work, we present a new manipulability measure that better expresses the manipulability of a mobile manipulator, because both the manipulability of the arm and the whole system are intrinsically considered.
Let us define the manipulability of the whole system, denoted as Ω p+a , and the manipulability of the robot arm, denoted as Ω a , to be the measures obtained using Equation (13) with the JacobiansJ and J a from Equation (5), respectively. Notice that by using the reduced JacobianJ, the platform's nonholonomic constraints are included in the manipulability measure of the whole system. Our proposed manipulability measure for mobile manipulators is defined as whereΩ p+a andΩ a are the normalized manipulabilities computed by dividing them by their respective maximum values over all the possible configurations of the system: The normalized values are used in our proposed measure to make sure that all its components are in the same scale. This normalization also makes the measures invariant to units and chosen reference frame [14,15].
Notice thatJ = J p J a S(q) ∈ R m×δ m in (5) indicates the relation betweenJ and J a . Hence, the singular values ofJ are not necessarily the same as the singular values of J a . There are cases where J a is rank-deficient whileJ is not. For these cases, the whole system manipulability Ω p+a = 0 even though the arm is in a singular configuration. Therefore, the measure Ω p+a fails to express the ability of the arm alone to provide end-effector's velocities in any direction. As a result, maximizing the measure Ω p+a might result in the poor manipulability of the arm, even though the manipulability of the whole system is preserved or improved. This is an issue that has been discussed in the literature [6,10].
On the other hand, our proposed measure Ω MM is defined by the product of the singular values ofJ and J a , i.e., the product of the manipulability ellipsoids ofJ and J a . Therefore, it encapsulates the abilities of the whole system and the arm to provide end-effector's velocities in any direction. If any of the singular values of J a orJ is zero, e.g., the arm is in a singularity, the measure Ω MM = 0. Hence, the measure Ω MM is a better representation of the manipulability of mobile manipulators because it intrinsically considers the manipulability of the arm and the whole system.
Using the product rule, the gradient of Ω MM is calculated as Analyzing the right hand side of (16), it is clear that the manipulability of the arm subsystem affects the gradient of the whole system. Likewise, the manipulability of the whole system affects the gradient of the arm subsystem. Therefore, maximizing Ω MM , i.e., setting ∇ q F(q) = ∇ q Ω M M in our solution (12) will simultaneously improve the manipulabilities of the arm and the whole system. In the experiments section (Section 5), a comparison is performed among the different objective functions for manipulability maximization to demonstrate the advantages of the proposed manipulability measure.
The reader may have noticed that because our proposed motion planning solution (12) multiplies W 1 2 with the gradient ∇ q F(q), the direction of the original gradient that maximizes the manipulability is changed. However, as our goal is to push the system away from singularities by gradually improving the manipulability rather than finding the shortest path towards its maximum value, it is sufficient to prove that the weighted vector W 1 2 u 0 also points in a direction that increases the manipulability as the unweighted gradient u 0 does. Since W 1 2 is positive semidefinite, we have u T 0 W 1 2 u 0 ≥ 0 for all u 0 . This implies that the weighted gradient also points in a direction that increases the manipulability.

Joint Position Constraints
Joint position constraints are common requirements in the motion planning of robotic systems. In this work, joint limits and self-collision avoidance are considered. These constraints are included in the proposed solution (12) using matrix W, which we define as the product of three terms as follows: where W Jlim is the weighting matrix for joint limits constraints, W Col is the weighting matrix for self-collision avoidance, and the matrix T q (known as the Jacobian normalization matrix), is used to normalize the velocity commands u as follows [12]: where u i max is the ith maximum velocity command. The weighting of this matrix handles the differences in units and scales among all the joints. For mobile manipulators, matrices W Jlim and W Col should be constructed only considering the arm. This is because there is no limit to the movement of the mobile platform along each degree of freedom. Likewise, the movement of the mobile platform cannot be used to avoid self-collisions because the platform moves the base of the arm. Therefore, in this work, the structure of both weighting matrices was designed to take this into account. To simplify the presentation, let W g represent either W Jlim or W Col . W g is a diagonal δ m × δ m matrix with the following form: where I 2 is an identity matrix of size 2 (for the case of the differential drive) and W a is a diagonal n a × n a matrix given by where each element of the diagonal matrix W a is defined using a performance criterion H(q).
In the next two subsections, two separate criterion functions for joint limits and selfcollision avoidance are defined. These two criterion functions have common properties. The values of H(q) and |∂H(q)/∂q i | become very large as the constraints are violated. When constructing matrix W a using these criterion functions, the joint movement towards or away from a constraint must be contemplated [11,16]. Under this consideration, the elements of W a are given by where ∆|∂H(q)/∂q i | is the change rate of |∂H(q)/∂q i | with respect to time, and is numerically calculated during implementation. With this choice, a value of one is assigned to w i , indicating that no penalty is imposed on the ith joint, if the ith joint is not moving (∆|∂H(q)/∂q i | = 0), or it moves away from a constraint (∆|∂H(q)/∂q i | < 0). On the other hand, w i tends towards zero if the movement of the ith joint gets closer to a constraint.
Hence, the element w i penalizes the movement of the ith joint by means of (12) if it moves towards a constraint and stops the joint if it is too close to it.

Joint Limits Avoidance
To construct the weighting matrix for joint limits avoidance W Jlim , a well-known criterion function [17] is used: where q i is the ith joint angle, q − i and q + i are its lower and upper limits, respectively, and γ is a scalar constant that adjusts the rate of change of H Jlim (q). This criterion function increases as the joint gets closer to its limits and goes to infinity at the joint bounds. The elements of the gradient of this function are given by Each element ∂H Jlim (q)/∂q i is equal to zero if q i is at the middle of its joint range and goes to infinity at either limit.

Self-Collision Avoidance
To construct the weighting matrix for self-collision avoidance W Col , an exponential function of the distance between two collision points is used as the criterion function [16]: where ρ > 0 controls the amplitude of H Col (q), and c 1 , c 2 > 0 control the rate of decay. This function has a maximum value when the distance d between two links is zero, and exponentially decreases as this distance increases. The distance between two collision points is defined as d(q) = p l1 − p l2 2 , where p l1 and p l2 represent the position vectors, referred to a common frame, of the collision points on two nonadjacent links. p l1 and p l2 can be calculated from the configuration vector q through forward kinematics.
The elements of the gradient of H Col (q) are given by where each of the partial derivatives is given by ∂d ∂q where J 1 and J 2 are the associated Jacobian matrices of p l1 and p l2 , respectively. The collision points are chosen from the surface of the links for which the collision distance is computed. For the case of potential collisions between the arm and the mobile platform, p l2 is picked as a point fixed on the surface of the mobile platform (p l2 does not move with respect to the arm). By using a frame attached to the mobile platform as the common frame, and selecting p l2 as a fixed point, (21) reduces to: When constructing matrix W a , the partial derivative ∂d/∂q in (19) was calculated using (21) or (22) depending on whether the collision is evaluated between two moving links or a moving link and a fixed one, respectively. Each element ∂H(q)/∂q i tends toward zero as the distance between two evaluated collision points increases, and tends towards infinity as the distance gets closer to zero. Because d = 0 is not admissible, i.e., the two links are in contact, the chosen points must contemplate a threshold.
Multiple pairs of potential collision points might exist in a robotic system. Let N c be the number of potential collision point pairs that are evaluated. A self-collision matrix W Col j is constructed for each distance d j , with j = 1, 2, . . . , N c . The weighting matrix that includes the contribution of each evaluated pair is finally obtained by With this combination, the ith diagonal element of the matrix W Col penalizes the movement of the ith joint. In contrast with the combination of collision weighting matrices proposed in [16], the combination (23) guarantees that the joints are stopped if they attempt to decrease any of the collision distances d j to a value of zero. Since some potential collisions are taken care of by the joint limits, the number of points considered in (23) is small, and therefore the computation cost can be relieved.

Joint Velocity Constraints
In a task space trajectory tracking problem, joint constraints at the velocity level are as important as joint limits and self-collision avoidance. These constraints include joint velocity boundary constraints and joint velocity limits. Satisfying joint velocity boundary constraints is a requirement for the end-effector to stop at the beginning and end of the task, i.e., the initial and final joint velocities must be zero. Joint velocity limits must be considered because the motion planning algorithm might generate joint velocity commands that are out of bounds in order to follow the task space velocity profile. These unfeasible velocities cannot be achieved by the real joints. Therefore, the trajectory tracking will fail if the joint velocity limits are not respected.

Joint Velocity Boundary Constraints
To satisfy joint velocity boundary constraints, the initial and final joint velocities must be zero. u p in (12) is directly associated with the end-effector's task; therefore, it implicitly satisfies the boundary constraints. However, u h in (12) does not necessarily satisfy them, because the manipulability maximization task is dependent on the mobile manipulator configuration at the initial and final poses. In this work, we propose using a time-varying self-motion magnitude to handle these constraints. The objective is to avoid non-zero values of u h only at the beginning and end of the trajectory. With this in mind, it is proposed to set the variation factor β(t) to increase from zero to one at the beginning of the trajectory, maintain a value of one for most of the trajectory, and then decrease from one to zero at the end of the trajectory.
Let t b be the blending time for β(t) to increase from zero to one, and t f be the trajectory execution time. To achieve a smooth transition, a 5th order polynomial β 1 (t) = a 0 + a 1 t + a 2 t 2 + a 3 t 3 + a 4 t 4 + a 5 t 5 is used when t < t b . The decrement in β(t) at the end of the trajectory, when t > t f − t b , is the complement of the polynomial β 1 (t) and is defined by After finding the polynomial coefficients, the self-motion magnitude variation factor is given by with: The blending time t b is chosen by the user. Figure 3 shows the shape of β(t) for t f = 15(s), t b = 2.5(s).

Joint Velocity Limits
To satisfy joint velocity limits, the maximum magnitude of self-motion is determined such that the velocity limit for each joint is not violated [18], namely: where u + i is the ith joint velocity limit with i ∈ 1, 2, . . . , δ m . Therefore, to avoid exceeding the joint velocity limits, α(t) must be selected such that it satisfies (25). The upper and lower limits of α(t) can be found using this equation. For each joint, it can be shown that the maximum and minimum values of α(t), denoted by α i max (t) and α i min (t), respectively, are given by Then, α max (t) and α min (t) for all the joints are: α min (t) = max α 1 min (t), α 2 min (t), . . . , α δ m min (t) , where α max (t) is the self-motion magnitude limit. In [18], the maximum value of α max (t) and minimum value of α min (t) are calculated for the whole trajectory, and the upper bound of α max (t) or lower bound of α min (t) is used as the step size to take advantage of the maximum admissible velocities. In our approach, a suitable initial value α s is selected through experimentation and the upper and lower limits of α are calculated for each time t. The value of α at time t is then given by This technique, in contrast with using the maximum/minimum self-motion for the entire trajectory, prevents sudden joint accelerations due to the use of maximum velocities in every step. Note that for some cases α(t) may be negative if α max (t) is negative, which means that the joint velocity limits cannot be satisfied without decreasing the manipulability of the system. The task can be executed as long as the system is still away from any singularity. The reader may have noticed that β(t) is in the denominator of (26) and (27), and as shown in Section 4.3.1, β(t) is zero at the beginning and the end of the trajectory. For these two cases, α min = −∞ and α max = ∞, but this does not cause instability in the system because by following (30), the value of α(t) is set to α s .
It is also possible to detect whether a task can be accomplished or not by using the limits of α(t). If the inequality α max (t) < α min (t) is true, then a suitable value of α(t) which avoids the joint velocity limits does not exist, because even the particular solution will violate them. In other words, the given task space trajectory cannot be accomplished without violating at least one of the joint velocity limits. For these situations, the task space trajectory must be replanned with a longer execution time t f or with lower end-effector's maximum velocities.

Experiments
Experiments were carried out to verify the efficacy of our scheme to solve the motion planning problem for trajectory tracking at the kinematic level. In this section, the results for the tracking of two trajectories, a Lissajous trajectory (see Section 5.4), and an elliptic trajectory (see Section 5.5), are analyzed. These trajectories were picked to demonstrate the ability of our approach to comply with the different constraints introduced in this manuscript while improving the manipulabilities of the whole system and the robot arm. For the Lissajous trajectory, a comparison of different objective functions for manipulability maximization is made to highlight the advantages of the proposed manipulability measure for mobile manipulators.
The experiments were performed using a 10-DOF NMM developed by the Industrial Technology Research Institute (ITRI), as shown in Figure 4. This NMM is composed of a differential-drive wheeled mobile platform, a prismatic joint mounted on top of the platform, and a Universal Robots UR5 6-DOF robotic arm attached to the prismatic joint. From this point forward, the UR5 arm's joints are denoted as q ai with i = 1, 2, . . . , 6 and the prismatic joint as z pj . The respective Denavit-Hartenberg (DH) parameters are shown in Table 1. The joint limits and joint velocity limits are shown in Table 2.
The Jacobian of the arm J a used for calculation ofΩ a in (15) was constructed only using the UR5 arm without the prismatic joint. If the prismatic joint were included, the arm would be allowed to stretch for most tasks since it would always be able to move the end-effector vertically using the prismatic joint, i.e., the manipulability is not affected when the arm is horizontally stretched. Stretching the arm for most tasks is an undesirable behavior, and therefore, the manipulability maximization of the UR5 arm without the prismatic joint is a suitable selection. The Jacobian J a is calculated with respect to the frame X a Y a Z a shown in Figure 4.
Due to the lack of a reliable positioning system, the odometry of the wheels was used in the experiments to compute the position and orientation of the mobile platform, and the forward kinematics were used to compute the end-effector's pose, which in turn was fed back to the motion planning algorithm. Therefore, the errors presented in the experiments are not the real-world errors, but the errors in the trajectory in which it is assumed that a reliable positioning system exists. Furthermore, as mentioned in the introduction section, the scope of this work is the motion planning of NMM for trajectory tracking at the kinematic level; therefore, the objective of the presented experiments is to validate the proposed algorithm at the kinematic level. Problems inherent to the dynamic behavior of the system, including vibrations of the mechanical structure, friction from the ground, etc., have an impact on the real end-effector tracking error, but they are out of scope of this paper and not considered in the algorithm. For these reasons, the simulation results of the position and orientation errors along the trajectory are also shown in this manuscript to highlight the performance of our motion planning algorithm.
An important remark is that the system vibrations due to the NMM mechanical structure greatly affected the performance of the motion planning algorithm when the system was close to its joint limits or self-collision constraints. This behavior was not observed when performing the simulations. To address this issue, a moving average filter with a window size of five was used to filter the gradients of the criterion functions for joint limits avoidance and self-collision avoidance. This filter diminished the impact of such vibrations on the motion planning algorithm.

Orientation Error for 6-DOF Tasks
In all the experiments, 6-DOF tasks were performed. For a 6-DOF task (m = 6), where the position and orientation of the end-effector are considered, the expression r d − r (mentioned in Section 3) has a specific definition that depends on the orientation representation, i.e.,r = r d − r does not hold for all orientation representations [1]. In this work, unit quaternions are used to describe the end-effector's orientation because of their efficiency and nonsingular representation for all orientations [1,[19][20][21].
A unit quaternion Q = [w + xi + yj + zk] is represented in scalar-vector form by Q = {s, v} with s ∈ R and v ∈ R 3 , where s and v are called the scalar and vector elements of Q, respectively. The desired and current pose are defined using unit quaternions for orientation as r d = P d Q d T and r c = P c Q c T , where P d = x d , y d , z d and If the desired orientation and current orientation are aligned, i.e., with zero orientation error, then ∆Q = {1, 0}. Consequently, it would be sufficient to define the orientation error to be ∆v. It is also important to follow a convention for the sign of the quaternion because Q = {s, v} and −Q = {−s, −v} represent the same orientation. A common convention is to keep the scalar quaternion element positive. Take this into account and the orientation error is defined as follows: Separating the position and orientation errors, the motion planning control law (6) is rewritten asṙ where K P and K O are positive definite diagonal 3 × 3 matrices. Note that e O in (32) is not an angle difference but its size represents the size of the orientation error that, as shown in [1,22], can achieve the convergence of the orientation error.

Evaluated Self-Collision Pairs
The types of self-collision can be significantly reduced by setting the joint limits properly. In addition, collisions among the first three links and the last three links can be taken care of by maximizing the arm manipulability, because the 6-DOF manipulator has poor manipulability if the arm is retracted to the point where the wrist is close to the base of the arm. Therefore, only the self-collision between the arm and the mobile platform needs consideration. Such a type of collision takes place when the elbow of the manipulator collides with the top of the platform, or when the wrist collides with the front of the platform. The distances associated with these potential self-collisions are depicted in Figure 5. Given that the platform is fixed with respect to the arm, setting frame X p Y p Z p (located on the center of the wheels' axle, as shown in Figure 4) as the reference frame allows to use (22), where p l1 is a point in the arm's elbow or wrist, correspondingly, and p l2 is a fixed point with respect to X p Y p Z p : As illustrated in Figure 5, to prevent the elbow collision, a collision pair between the z component of p elbow with respect to frame X p Y p Z p and a point located at 0.5 (m) from the origin of frame X p Y p Z p in the Z p direction is selected. The distance for this pair is named d elbow . To prevent wrist collision, a collision pair between the x component of p wrist with respect to frame X p Y p Z p and a point located at 0.37 (m) from the origin of frame X p Y p Z p in the X p direction is selected. The distance for this pair is named d wrist . The wrist collision is only evaluated if the wrist height h wrist , the z component of p wrist with respect to frame X p Y p Z p , is lower than 0.5 (m), otherwise, its associated weighting collision matrix is assigned to identity. This was done to avoid restricting the movement of the joints that reduce d wrist when the wrist is above the top of the mobile platform. All these parameters were chosen based on the physical dimensions of the NMM, and the second point in each of the collision pairs was selected so that when d elbow and d wrist are zero, a gap still exists between the potentially colliding links.

Common Parameters of the Simulations and Experiments
For all the experiments, the selected feedback gain matrices in (33) are K P = 10I 3×3 and K O = 20I 3×3 , where I is an identity matrix. The initial value of α(t) in (30) is set to α s = 3, and the blending time is set to t b = 0.2t f for the variation factor β(t) in (24). The parameters for self-collision avoidance in (20) are ρ = 1 × 10 −3 , c 1 = 50 and c 2 = 1. The starting configuration of the robot arm is q a = 0 −80 110 −120 −90 0 ( • ). Furthermore, a sampling time t s = 0.02(s) was used for the simulations and experiments.

Lissajous Trajectory
The Lissajous trajectory was picked to demonstrate the ability of the proposed approach to track a trajectory for which the nonholonomic constraints greatly affect the movement of the system. The proposed Lissajous trajectory is defined by where [x 0 , y 0 , z 0 ] T is the end-effector's starting position, s(t) is the trajectory timing variable, and A, B and C define the size of the path. The orientation is set to be the same for the entire trajectory, i.e., Q(t) = Q 0 . Notice that this is a 6-DOF trajectory because the orientation is constrained for the entire trajectory. The starting configurations of the mobile platform and prismatic joint for this trajectory were set to x p = −0.1 (m), y p = −0.13 (m), θ p = −90( • ) and z pj = 0.2 (m). With this configuration, the end-effector's initial pose is given by P 0 = −0.009 −0.6491 0.9884 (m) and Q 0 = {0, 0i + 1j + 0k}. A trapezoidal velocity profile [1] was used for the derivative of the timing variableṡ(t). The parameters for the path size were set to A = 1.3 (m), B = 1.3 (m) and C = 0.27 (m), and an execution time t f = 64(s) was chosen. Figure 6a,b show snapshots of the NMM's movement along the Lissajous trajectory, in simulations and experiments, respectively. Figure 7 compares the trajectory tracking results between the simulation and the experiment. In the simulation, the position and orientation errors are kept small during the whole trajectory, as shown in Figure 7b,c. This demonstrates the good tracking performance at the kinematic level of our proposed motion planning algorithm. The observed larger errors in the experiments, as exhibited in Figure 7e,f, are due to the vibrations of the system during the experiments. Furthermore, as mentioned before, the control of the dynamic behavior of the system is beyond the scope of this work. Nevertheless, the position errors in the experiments are kept below 2 × 10 −3 (m) and the orientation errors below 1.5 × 10 −3 . Likewise, we observe that the obtained trajectories in the simulation (Figure 7a) and the experiment (Figure 7d) are fairly similar. Figure 8 illustrates the experiment results. Note that the negative joint velocity limit −u + i is denoted as u − i in the pertinent figures. The manipulabilities of both the arm and the complete system are kept high during the execution of the trajectory, and their final values are higher than at the start of the trajectory, as shown in Figure 8b. It is important to remark that there are no potential self-collisions during the execution of this trajectory. Even though d wrist is negative at the beginning of the trajectory, as shown in Figure 8c, the wrist is above the top of the mobile platform and therefore the joint movements were not restricted, as explained in Section 5.2.   Figure 8f, is restricted in the time interval t = (29, 32)(s) to respect its lower limit. This time interval corresponds to the trajectory section between snapshots number three and four in both Figure 6a,b. For this section of the trajectory, the movement of the end-effector in the XY plane is taken care of by the platform due to the restriction imposed to q a1 . We observe in Figure 8e that the velocity limits of the prismatic joint are respected. The remaining joints do not reach their respective limits as shown in Figure 8f-h. Furthermore, as seen in Figure 8d,e,i, the velocity profiles for all the joints are smooth and satisfy the boundary constraints, i.e., the initial and final joint velocities are equal to zero.
To demonstrate the advantages of the proposed manipulability measure Ω MM , a comparison of the manipulability maximization results using different objective functions in simulations of the Lissajous trajectory tracking is shown in Figure 9. Here, a task is considered as failed when none of the constraints are satisfied, and thus the simulation is stopped. These results were obtained by using the same parameters shown in Section 5.3 with the only change being the objective function.
The Lissajous trajectory tracking fails when the objective function is set to maximize the manipulability of the arm, i.e., F(q) =Ω a , as shown in Figure 9a. The manipulabilities of both the arm and the whole system start a fast decay after t = 44s. This decrement in the manipulabilities is the result of restricting the homogeneous solution to respect the joint velocity limits. Consequently, the arm manipulability is not maximized anymore and the system moves towards singularity, ultimately failing the task. The results of setting the objective function to the manipulability of the whole system, i.e., F(q) =Ω p+a , are shown in Figure 9b. In this case, the task succeeds, however, the manipulability of the arm deteriorates, and at the end of the trajectory has a value close to zero. This is the behavior discussed in previous studies [6,10].  Figure 9c shows the results of maximizing a linear combination of the manipulabilities of the arm and the whole system, i.e., F(q) = 0.5Ω p+a + 0.5Ω a . Figure 9d shows the results of maximizing the proposed manipulability measure for mobile manipulators, i.e., F(q) = Ω MM from Equation (15). Both objective functions preserve the manipulability of the arm. However, it is clear that maximizing the proposed measure has better results overall. First, the manipulability of the arm is higher along the trajectory in Figure 9d compared to the arm manipulability obtained with the linear combination in Figure 9c. Secondly, the obtained final manipulabilities have improved with respect to the initial ones in Figure 9d. On the other hand, the linear combination improves the manipulability of the whole system, but not that of the arm with respect to the values at the start of the trajectory, as shown in Figure 9c.

Elliptic Trajectory
The elliptic trajectory was picked to demonstrate the ability of the proposed approach to comply with joint velocity limits and to prevent self-collisions. This trajectory consists of moving the end-effector from an initial pose r 0 = P 0 Q 0 T to a desired final pose r d = P d Q d T using an elliptic path. The trajectory position is defined by where A, B, c x , c y and s 0 are calculated using the XY coordinates of P 0 and P d , such that the elliptic path is centered at the closest point to the origin, while it covers a 90 • angle between the start and end points. The trajectory's Z coordinate follows a straight line between the points (z 0 , s 0 ) and (z d , s d ), where z 0 and z d are the Z coordinates of P 0 and P d , respectively; s 0 and s d are the starting and ending angles of the elliptic path, respectively; and m in the equation above is the slope of this line. A fifth-order polynomial profile [1] was used for the timing variable s(t). The orientations and rotational velocities along the trajectory were computed using quaternion polynomials [23]. This technique has two main advantages: a smooth velocity profile is obtained, and the rotational velocities and accelerations in the task space are included as boundary constraints.
The starting configuration of the mobile platform and prismatic joint for this trajectory was set to x p = −1.3 (m), y p = 0.56 (m), θ p = 0( • ), and z pj = 0.24 (m). With this configuration, the end-effector's initial pose is given by P 0 = −0.781 0.669 1.028 (m) and Q 0 = {0, 0.7071i − 0.7071j + 0k}. The final pose was selected to have position P d = 1.55 −1.0 0.26 (m) and orientation Q d = {0.2706, 0.6533i + 0.6533j − 0.2706k}, and an execution time t f = 20(s) was chosen. Figure 10a,b show snapshots of the NMM's movement along the elliptic trajectory, in simulations and experiments, respectively. Figure 11 compares the trajectory tracking results between the simulation and the experiment. The small position and orientation errors obtained in the simulation (Figure 11b,c) again demonstrate the good tracking performance of the proposed approach. In the experiment, the position errors are kept below 1.5 × 10 −3 (m) and the orientation errors are kept below 1 × 10 −3 as depicted in Figure 11e,f. Once again, the vibrations of the system played a role in the larger errors seen in the experiment.  Figure 12 illustrates the experiment results. The manipulabilities of both the arm and the complete system are again improved at the end of the trajectory compared to the initial configuration, as shown in Figure 12b. However, we notice how both manipulabilities decreased during the time interval t = (8, 15)(s). This behavior is due to the value of α(t) (manipulability maximization step size) being negative for this span of time in order to respect the joint velocity limits, as discussed in Section 4.3.2. Notice that v p (Figure 12d) andż pj (Figure 12e) are kept at their maximum values during this span of time. It would not be possible to track this particular trajectory while respecting the joint velocity limits without stretching the arm; hence, the manipulabilities of the complete system and the arm decreased.
We observe in Figure 12c how the self-collision distance of the elbow gets close to zero. In this trajectory, the elbow needs to get close to the platform in order for the end-effector to track the desired trajectory. The motion planning algorithm gradually stopped the elbow to prevent the collision with the platform, as shown by the slow decay of d elbow . To achieve this, the movement of the prismatic joint was restricted, as shown in Figure 12e. The wrist also moves towards the front of the platform, as shown by the decrement in d wrist in Figure 12c. This potential self-collision is also prevented by limiting the movement of q a2 and q a3 , as shown in Figure 12f,g. Figure 12d-i show that the joint limits, joint velocity limits and joint velocity boundary constraints are also satisfied. All the joints were kept within their respective limits, as shown in Figure 12e-h. As depicted in Figure 12d,e, the platform's linear velocity and the prismatic joint's velocity are kept within their limits. Once more, the velocity profiles are smooth and the initial and final velocities for all the joints are zero. Notice that all the joint limits and velocity limits are respected (panels (d-i)). The potential self-collisions described in Section 5.2 are prevented (panel (c)). Furthermore, the manipulabilities of both the arm and the whole system are improved (panel (b)). See Section 5.5 for a detailed discussion of this figure.

Discussion
A scheme was proposed to solve the motion planning of NMMs considering joint limits, joint velocity limits, joint velocity boundary constraints, and self-collision avoidance while maximizing the manipulabilities of both the robot arm and the whole system.
The proposed solution uses a weighted input velocity vector and a weighted Jacobian to penalize the movement of joints that get close to a position constraint. A proposed quadratic cost function is minimized when solving the motion planning problem for redundant NMMs. This cost function includes a secondary task to be satisfied that is also weighted to comply with the position constraints. The maximization of an introduced manipulability measure for mobile manipulators is used as the secondary task to push the system away from singularities. In the experiments section, it is demonstrated that maximizing this new measure simultaneously improves the manipulability of the whole system and that of the manipulator alone.
This work focuses on the motion planning for trajectory tracking at the kinematic level, which must not only comply with joint position constraints, but must also respect joint velocity constraints and joint velocity boundary constraints. Joint velocity boundary constraints are satisfied by varying the magnitude of the homogeneous solution at the start and end of the trajectory. The manipulability maximization at these points is not necessarily zero; hence, using an adequate variation in the magnitude of self-motion is needed. Joint velocity limits are satisfied by evaluating the maximum allowable self-motion for each joint. Using this information, the step size of the gradient ascent/descent is limited when required, and consequently, the joint velocity limits are not exceeded.
The experiments for 6-DOF trajectories were conducted to verify the efficacy of the proposed scheme. The results demonstrate that the proposed approach can solve the motion planning problem for NMMs to perform trajectory tracking at the kinematic level while considering the constraints required for real implementation including manipulation capability preservation or improvement.
The experiments designed in Section 5 consider an open environment without obstacles, because this is the scope of our manuscript. However, the proposed solution can be extended to prevent collisions with obstacles by including collision pairs between the robot arm and these obstacles, using the same definitions as in Section 4.2.2. This will penalize the movement of the arm's joints that get closer to an obstacle in the environment. Nevertheless, in the case of the platform, stopping it is not an efficient approach. In this case, an additional task can be added to the solution to push the platform away from the obstacles. One way to achieve this is by using a task priority scheme [24] using the nullspace of the motion planning algorithm.
Even though our work focuses on NMMs, our approach can be effortlessly adapted for use with holonomic mobile manipulators. Future work will focus on dynamic modeling and controller design to deal with system vibrations and tire slip such that the tracking errors of the system can be reduced.