Rigid-body inverse dynamics of a spatial redundantly actuated parallel mechanism constrained by two point contact higher kinematic pairs

. This paper presents a comparative study of the rigid-body inverse dynamics of a spatial redundantly actuated parallel mechanism constrained by two point contact higher kinematic pairs (HKPs). Firstly, its constrained motions are analysed comprehensively, then four different models are formulated by the generalized momentum approach and the Lagrange-D ’ Alembert formulation to explore its inverse dynamics. In each method, the ﬁ rst model is built by employing the method directly to the mechanism. In the second model, the dynamic model of its non-redundantly actuated counterpart free of HKPs is built by this approach ﬁ rst, then the constraints from HKPs are modelled, to ﬁ nally reach the model of the redundantly actuated parallel mechanism (RAPM) where that of its counterpart is utilised as the core. The four models give rise to equivalent numerical results, and the second model in both methods of the RAPM can alleviate the strong coupling between the parasitic motion variables and degrees of freedom (DOFs), boosting the computational speed as fast as that of its non-redundantly actuated counterpart without simpli ﬁ cation or loss of accuracy. The comparisons between the mechanism and its counterpart validate that the HKP constraints greatly increase the computational complexity, and the torques required by the parasitic motions of the end effector are signi ﬁ cantly smaller than those by the corresponding DOFs.


Introduction
In the food industry, there is a great interest in evaluating the time-varying dynamics of newly developed food textures during the chewing process. To this end, a machine that can accurately replicate the human-like chewing behaviours in a biomimetic fashion can come into play a significant role in terms of the three-dimensional (3D) chewing motions and bite forces. Inspired by the masticatory system of human beings, a spatial redundantly actuated parallel mechanism (RAPM) constrained by two point contact higher kinematic pairs (HKPs) has been developed [1,2]: the base is the skull, the six RSS (revolutespherical-spherical) kinematic chains are the primary chewing muscles, the end effector is the mandible, and the two point contact HKPs are the left and right temporomandibular joints (TMJs), respectively. The underlined letter means the joint is active in the chain. At the early stage of our study in the bio-inspired masticatory robot, the six RSS kinematic chains were designed to play the role of primary chewing muscles in the human's masticatory system, as in Chapter 4 of [3]. However, during the motions of the lower mandible, it is not only driven by chewing muscles, but also guided at TMJs at its left and right condyles. In this regard, to better replicate the chewing system and produce human-like chewing behaviours, TMJs were modelled using two HKPs as in [1,2]. Readers interested in the human masticatory system and chewing robotics can refer to [3] for a comprehensive description.
Before the applications of the RAPM in practice, developing a precise and computationally efficient inverse dynamic model is a prime requirement. In fact, in our previous work [2], an initial attempt at the rigid-body inverse dynamics of the RAPM of interest has been made via the hybrid of the Newton-Euler's law and the Lagrangian equations. However, due to the complicated constraints from HKPs at the end effector in the RAPM, the model was quite sophisticated and error-prone. For instance, the forces at HKPs without friction effects are ideal constraints on the end effector which cannot be included in the Lagrangian formulation, whilst they must be considered in the Newton-Euler's law. Meanwhile, the Lagrange formulation introduces many unknown Lagrange multipliers, greatly increasing the number of variables in the equations of motion (EOMs). In these regards, finding efficient solutions to its inverse dynamics is a strong motivation for the study in this paper.
From the literature review, the conventional methodologies in dynamics of serial mechanisms and PMs without redundant actuation have been successfully applied in a wide range of RAPMs. For instance, the Newton-Euler's law was used in a planar 3-DOF RAPM to find its maximum dynamic load-carrying capacity and the impact of actuation redundancy [4]. The Euler's equation was applied in a 4RRR planar RAPM to achieve its shaking force/moment balancing [5]. In [6], this law was employed to solve the inverse dynamics of the spatial RAPM with parallelograms inside a serial-parallel hybrid machine. One can find that under this law, the free body diagram must be drawn to each body, and there are a number of intermediate variables associated with the constraint wrenches in the joints. Therefore, its computational complexity and efficiency would be negatively influenced by the increased number of actuators and moving links in RAPMs. In [7], the Lagrange formulation with unknown Lagrangian multipliers was used to address the dynamics of a planar RAPM with two translational DOFs. The multipliers representing the magnitude of constraint forces were eliminated using the null space of the differential matrix of the closed-loop constrained equation. The identical procedure was also employed in [8,9]. The utilisation of the Lagrange equations in the modelling process is accompanied by the introduction and elimination of Lagrangian multipliers, and it involves a number of partial derivatives. As a result, the computational demands are not light. The principle of virtual work was used in evaluating the superiority of redundant actuations in the dynamic performance of an 8PSS spatial RAPM in [10]. In [11], this principle was applied to compare the dynamic performance of planar PMs with both actuation and kinematic redundancies. It was employed in [12] to calculate the joint torque to compensate for disturbances to the mobile platform of a dual-arm underwater robot. It also has been used in [13][14][15] where the EOM was used for model-based control. In this principle, the ideal constraint wrenches that do not generate virtual work are not taken into account, and it does not require as many partial deviations as in the Lagrange formulation. As such, compared to the foregoing two methods, its computational cost is lower.
Apart from them, another two methods in the literature attract our attention. The first is the newly developed generalized momentum approach in [16,17]: it has been successfully applied in both the conventional Gough-Stewart platform and a 6PUS (prismatic-universal-spherical) PM. In this approach, the kinetic and the gravitational components of generalized forces acting on a moving body are computed by the differentiation of its generalized momentum and potential energy, respectively, and they can be expressed by the generalized coordinates. In addition, the contributions to the total actuating torques by each moving body can be found separately. As such, compared to the classical Newton-Euler' law, it is characterised by being well-structured, straightforward, and systematic. The second method is the classical Lagrange-D'Alembert formulation, by which the dynamic models of the planar 2-DOF RAPM from [8] and 3-DOF purely translational PM were built in [18,19], respectively. In this approach, the dynamic model of each rigid body is firstly built by the Lagrange formulation, then following the D'Alembert principle, the EOM of the mechanism is reached. Compared to the procedure by the Lagrange's equations, it is not indispensable to introduce and then eliminate the unknown Lagrangian multipliers. In addition to the two methods' respective advantages, they have two common features: the constraints imposed by the closed loops in the mechanism are embedded via the velocity relationships between the end effector and the chains, obviating the use of explicit constraint equations; they do not take constraint wrenches in joints that do not generate virtual work into consideration. As such, a minimal set of variables and equations is generated and it is straightforwardly beneficial to the model-based motion control.
Even though the successful applications of these two methods in RAPMs are not quite extensive, due to these foregoing merits, an attempt is going to be made to extend them to the RAPM with HKPs which has a sophisticated topology, and this is the first reason to use them in this paper. The second is, the computational cost of these two methods in RAPMs has not been compared to date, despite their common advantages. Thus, it is necessary to find which is faster. The third reason is, the software at hand i.e., Matlab does not support redundant actuation mode in simulations of inverse dynamics. Thereby, it is required to use two methods simultaneously to cross-check the correctness of the numerical results, and the foregoing two methods are adopted in this paper.
In the meantime, from the view point of mechanism, there are primarily two methods to implement actuation redundancy in PMs as stated in [20]: the first is adding identical kinematic chains as the existing ones in a PM, and the second is actuating some of the passive joints of PMs. By a sharp contrast, from the analysis in Section 3 in this paper, redundant actuations are realised by imposing two point contact HKPs at the end effector of a 6RSS PM. Evidently, by this process, both the RAPM and its nonredundantly actuated counterpart, namely, the 6RSS PM have identical numbers of moving bodies and actuators, respectively. Nonetheless, intuitively, the kinematic variables in the RAPM are more strongly coupled, leading to more complex kinematics and dynamics, and the computational demands would be higher than that of its counterpart. Furthermore, studies on this sort of RAPMs with HKP constraints are quite rare: from the literature, the publications are almost about RAPMs with lower kinematic pairs and their end effectors are not constrained by the base directly. Thus, under the adopted methods, finding a short-cut that can alleviate the coupling among variables to increase the computational speed of the dynamic model is the second strong motivation for conducting the study in the paper. In fact, the discovery of this short-cut would indicate the nature of HKP constraints in the inverse dynamics of the RAPM is explored deeply.
In this paper, the rigid-body inverse dynamics is studied. It is assumed that all bodies and joints are rigid, frictionless, and free of clearances in the mechanism. The inertia of the spherical joints is rather small and then ignored. The remainder of the paper is organized as follows: a detailed description of the mechanism is given in Section 2. In Section 3, the kinematics of the mechanism including the constrained end effector and the kinematic chains is analysed comprehensively. The inverse dynamics of the RAPM is analysed by the generalized momentum approach in Section 4, in Model 1, the approach is directly employed to the RAPM, while in Model 2, firstly the dynamic model of the 6RSS PM is built by this approach, then the two constraints from HKPs are modelled, finally, the model of the RAPM is built by incorporating the model of the 6RSS PM as its core. In Section 5, an identical procedure as in Section 4 is implemented by using the Lagrange-D'Alembert formulation. Section 6 is devoted to the numerical computations and discussions in the four models by the two methods, and the influences of the HKP constraints. Finally, some conclusions are made in Section 7.
The three contributions of this paper are: -By using the inverse dynamic model of the 6RSS PM as the core of that of the RAPM can alleviate the coupling in its parasitic motion variables and increase its computational speed as fast as that of the 6RSS PM. -The insight into the nature of HKP constraints on the inverse dynamics has been provided clearly in terms of the modelling process, the computational cost, and the numerical results. -The applications of the generalized momentum approach and the Lagrange-D'Alembert formulation have been extended to the RAPM constrained by HKPs.

Description of the mechanism
The scheme of the RAPM constrained by two HKPs is illustrated in Figure 1. The maxilla (i.e., the base) is fixed on the ground and the movable mandible (i.e., the end effector) is connected to the base by six independent kinematic chains. The maxilla, to which the inertia frame {S} is assigned, is not shown in the figure for a clear exhibition of moving bodies. This frame consists of a horizontal X S -Y S plane perpendicular to the vertical Z S axis. A frame {M} is established at the mass centre O M of the end effector. The origins and orientations of {S} and {M} overlap when the mechanism is at the home position, that is, the maxilla and the mandible are in the occlusal state. The origin O M is used as the reference point to describe the mandibular translations, and its orientations with respect to {S} are expressed by Bryant angles. In each chain, the crank G i S i (i =1,…,6) is driven by a basemounted rotary actuator, whose shaft connects it with a rotational joint at G i , and a coupler S i M i that joins the crank and the end effector via two spherical joints at its two ends S i and M i , respectively. The layout of six chains is in accordance with that of the six primary chewing muscles. The rotation of the ith actuator with respect to {S} is described by the actuator frame{C i }attached at G i . In it, the.X C i .axis is directed from G i to S i , theZ C i axis runs through the driving shaft of the actuator, and theY C i axis completes the frame, obeying the right-hand rule. A frame {N i } is attached at the mass centre E i of S i M i to describe its motions with respect to {S}. The X N i axis points from S i to M i , theY N i axis is parallel to the cross product of two unit vectors defined along theX N i and X S axes, and theZ N i axis is defined by the right-hand rule. Two HKPs modelling the left and right TMJs are formed by the two condyle balls being in contact with the articular surfaces. The two contacts between the left and right condyles of the mandible and the maxilla at TMJs, are realised by the two condyle balls being always constrained within a condylar socket in the mechanism prototype, as shown in Figure 2. The width of the socket equals the diameter of the condyle ball; thus, it can always guarantee the point contact during the movement of the mandible. By this design, the motion of the condyle ball centre is always constrained onto a surface, which is offset from the upper and lower surfaces of the socket by the ball radius. Thereupon, it is clear that the end effector is actuated by six chains and constrained by the base at the two HKPs simultaneously. Besides, the end effector is shown in Figure 2d. the green component is a manufactured mechanical part, while the red lower mandible and golden teeth are scanned from a human cadaver by computer tomography in the School of Dentistry at the University of Otago, Dunedin, New Zealand [3]. The prototype of the mechanism is provided in Figure 3.

Kinematics of the mechanism
A priori development of accurate kinematic analysis is fundamental to the inverse dynamics. In this section, the kinematics of the end effector and the chains in terms of their displacements, velocities, and accelerations is analysed comprehensively.

Constrained motions of the end effector
A second-order surface was used as the workspace of the centre of the condylar ball in [1,2,21,22]; its shape and size was designed by referring to [23][24][25]. However, on the one hand, when the mechanism tracks the real chewing trajectories of healthy human subjects, for instance, the one in Section 6 which serves as an illustrative example in this paper, the surface where the path of the condylar ball centre is situated can be approximated as a flat one. On the other hand, it is quite difficult to derive explicit analytical expressions of the parasitic motions in the RAPM under the second-order surface. In these regards, in this paper where the chewing system is explored from the viewpoint of mechanical dynamics and with an emphasis on the constrained dynamics of the end effector, the surfaces in {S} where the left and right condyle ball centres T L and T R slide on are designed as flat (unit: mm): From the Kutzbach-Grübler criterion, the mechanism now has four DOFs, but the information on which four to choose is not given. They are to be derived from a rigorous computation below: the coordinates of T i (i = L,R) in {S} can be expressed as andR Z (g)are three rotation matrices about the X M , Y M , and Z M axes by three Bryant anglesa, b, andg, respectively. It is worth noting that in this paper, a matrix/vector/scalar in local frames owns a leading superscript on its left to denote the specific frame it refers to, but those in {S} omit their superscripts for the sake of convenience and clarity. The six motion variables of the end effector are grouped and expressed as From equation (2), one can obtain where M S R i;: Because of the left-right symmetry ofMO M T L and-MO M T R in {M}, a summation and a subtraction of the two equations in equation (5) sidewise yield where MO M T L (1) and MO M T L (3) are the first and third terms ofMO M T L , respectively. From these computations, it is found Z and g are transferred from DOFs to parasitic motions and they are functions of q EE , which is a 4 Â 1 vector by grouping four DOFs as and it constitutes the task space of the mechanism. To characterise the instantaneous configuration of the mechanism, equation (3), or both equations (6) and (7) ad hoc are needed. In other words, the RAPM can still perform motions in six directions with four DOFs and two parasitic motion variables. Regarding this, redundant actuations in the mechanism are essentially caused by constraints from the base directly onto the end effector, which is completely different from the two methodologies mentioned in [20]. It is also worth noting that though the workspace of the centre of the condylar ball is simplified as a flat surface as in equation (1), a strongly nonlinear and sophisticated relationship between Z/gand q EE in equation (6) can be observed.  After figuring out the DOFs of the end effector, its motions can be defined below. The angular velocity is where 0 3 is 3 × 3 the zero matrix, and M J denotes the 6 × 4 Jacobian matrix between X EE and q EE , namely From it, the second time-derivative of X EE can be derived as The translational velocity of the end effector is expressed by where E 3 is the 3 Â 3 identity matrix.

Kinematics of the ith chain
The inverse kinematics of the mechanism, i.
consists of a system of six decoupled equations expressed byq EE , has already been derived in Section 3.2 of [22]. Nonetheless, the motions of the coupler S i M i (i = 1,…,6) are still needed for its rigid-body dynamics. Due to the two spherical joints at S i and M i , the coupler can rotate around the three orthogonal axes of {N i }. The rotation around X N i axis is a passive DOF for it is not controllable; this rotational range is quite small thanks to the physical restrictions from the used spherical joints in the mechanical design, however. In these regards, it is assumed that there is no axial rotation in the coupler. As such, two Euler angles b i and g i around theY N i and Z N i axes, respectively, are used to express the rotation of S i M i in {S} in terms of the two rotational matrices R Y (b i ) and R Z (g i ). Thereafter, the coordinate vector of the coupler can be expressed as where N i0S R is the orientation of S i M i in {S} at the initial configuration of the mechanism at hand, andkS i M i kis the length of S i M i . From the geometry of the mechanism in {S}, the vector of the coupler can also be found from the difference in the position vector of M i and S i : Substituting equation (12) and where n i1 n i2 n i3 , which consists of the joint space of the ith chain and completely specifies its configuration. It is highlighted thatq r i is the function of q EE , i.e.,q r i ¼ q r i q EE ð Þ. To derive the relationship between the first timederivatives ofq r i andq EE , equation (13) can be rewritten as The two sides can be expressed by q r i and q EE , respectively. The first time-derivative of equation (15) yields where Moreover, one can find that where M 3i ¼ M À1 1i ⋅M 2i , and the second time-derivative ofq r i is given as So far, _ q ri _riand __ q ri rihave been derived as quantities intimately associated with the motions of the ith chain.

The angular velocity of
The coordinate vector of the mass centre E i in {S} is The translational velocity of E i can then be found by the differentiation of equation (21) with respect to time as Upon substitution of equation (17) into equation (22), it yields As far as the crank G i S i is concerned, its mass centre G i is passing through the rotational shaft of the actuator, as such where C iS R :;3 ð Þ is the third column of C iS R. Heretofore, the motions of all moving bodies in the mechanism have been explicitly expressed. These are all the ingredients needed for the construction of the EOM based on the two dynamic methods in the following two sections. In each of them, two models are built by the generalized-momentum approach and the Lagrange-D'Alembert formulation, respectively. In Section 4, the generalized-momentum approach is directly used to the RAPM in Model 1; in Model 2, by following an intuitive procedure, the EOM of the 6RSS PM free of HKPs is firstly written, then the two constraints from HKPs are modelled, to finally reach the EOM of the RAPM where the model of the 6RSS PM is used as its core. An identical procedure is implemented in Section 5 by the Lagrange-D'Alembert formulation. The structure of dynamic modelling in the following two sections are given in Figure 4.

Singularity analysis
To analyse the singular configurations of the mechanism, a 6◊1 constraint vector is defined as where and their sizes are 6◊6 and 6◊4, respectively. It is also noted that the ith(i = 1,…,6) row of L (u, q EE ) is the function of u i and q EE , then the non-diagonal entries of J A must be zero, i.e., J A is a diagonal matrix.
From Section 5.3 of [26], there are three sorts of singularity conditions: 1: An inverse kinematic singularity occurs when J A is rank-deficient, i.e., Thus, at least one diagonal entry of J A is zero. In this case, a motion of the actuators does not lead to the displacement of the end effector.

A direct kinematic singularity occurs when
In this case, the end effector can possess uncontrollable motions in some directions while all actuators are locked.
3. Combined singularities occurs when both equations (27) and (28) are valid. In this case, the end effector can remain stationary while the actuators undergo some infinitesimal motions. On the other hand, it also can undergo infinitesimal displacements in some directions while all actuators are locked.

Generalized momentum approach
In writing the EOM by this approach, it must be kept in mind that the forces acting at the end effector via the two HKPs are ideal constraint forces rather than active forces, as such, they do not appear in the formulation. In this section, Model 1 is built by directly applying this approach to the RAPM, and in Model 2, the dynamic model of the non-redundantly actuated counterpart, namely, the 6RSS PM is firstly built, then the HKP constraints are modelled, to arrive at the dynamic model of the RAPM finally.

Model 1
Under this approach, the kinetic and gravitational components of the generalized forces of each moving body is obtained by differentiating its generalized momentum and potential energy with respect to time and q EE , respectively. The EOM of the RAPM is acquired by the law of power equivalence: the power generated by the generalized forces of all rigid bodies in the mechanism is equivalent to that by actuating torques from actuators and external wrenches from the environment.

The end effector
On the one hand, the generalized momentum of the end effector is computed as is the 6 Â 6 inertia dyad, m EE and I EE are the mass and the inertia matrix with respect to O M and expressed in {S}, respectively. By equations (8) and (11), it can be found that The kinetic component of the generalized force acting on the end effector can be computed as On the other hand, the potential energy of the end effector is where g = 9800 m/s 2 is the gravitational acceleration. The gravitational component of the generalized force is 4.1.2 The ith chain As far as the coupler S i M i of the ith kinematic chain is concerned, by virtue of equations (20) and (23), its generalized momentum is computed as where The kinetic component of the generalized force acting on S i M i is computed as The potential energy of S i M i is For the crank G i S i , its generalized momentum is where tor, and the scalar I G i S i is the rotational inertia of G i S i around theZ C i axis of frame{C i }. As such, its kinetic component of the generalized force is Because G i is a fixed point in {S}, its potential energy is constant and computed as in whichO S G i (3) is the third term of O S G i . Thereby, its gravitational component of the generalized force is

The overall mechanism
In the entire mechanism, the instantaneous power generated from the sum of the generalized forces on all the moving bodies is balanced by that from the input torques and bite force. As such, one can write wheretis the6 Â 1actuating torque vector by the six actuators, andW F B is the6 Â 1wrench vector by the bite force F B . Substituting equations (8), (11), and (17) into equation (41), and rearranging the term byW F B give rise to Jacobian matrix mapping q_ EE into u_, i.e., u_=J u1 ⋅q_ EE , and is the first row of. Because q EE is the independent generalized coordinate vector of the entire mechanism under study, q_ EE is free to vary, thus one can obtain which is the inverse dynamic model of the RAPM and it contains four equations and six unknowns.

Model 2
In this model, the EOM of the 6RSS PM without HKP constraints is firstly written, then the constraints from two HKPs are modelled, to finally reach the model of the RAPM, where the model of the 6RSS PM is utilised as its core. However, it must be kept in mind that now the kinetic and gravitational components of generalized forces are functions of X EE ; X EE ; X EE rather than q EE ; q EE ; q EE . In the 6RSS PM, the step about the instantaneous power equivalence among the sum of the generalized forces on all the moving bodies, the input torques, and the bite force is as follows: are the kinetic and gravitational components of the generalized force acting on the end effector, the ith coupler and the ith crank, respectively. They can be computed in an identical manner as their counterparts in the last section; thus, the process is omitted to save pages. After the two constraints from HKPs are exerted on the end effector, its DOFs are changed from X EE to q EE . Then using equations (8), (9), and (11), equation (44) can be rewritten as vector by the bite force F B , J u2 is the 6 Â 6 Jacobian matrix mapping X_ EE to u_, namely, u_=J u2 ⋅X_ EE . Because q_T EE is free to vary, it can be rewritten as which is the model of the RAPM with four equations and six unknowns. It must be emphasised that even though the generalized forces in the bracket are ultimately functions of q EE ; q EE ; q EE , they are directly expressed and computed by X EE ; X EE ; X EE . In this manner, the strong coupling in the parasitic motion variables Z and by q EE is completely alleviated, reducing the computational complexity to a huge degree, which will be exhibited in Section 6.2.1.
One can easily find that, by substituting only equations (8) and (11) into equation (44), and "deleting" the independent vector X_ EE produces which is the EOM of the 6RSS PM without HKPs. That is to say, from the comparison between equations (46) and (47), the EOM of the RAPM can be easily generated by left-multiplyingM T J at the two sides of that of the 6RSS PM, which is used as the core of the dynamic model of the RAPM.

Lagrange-D'Alembert formulation
In writing the EOM of the RAPM by this formulation, there are two models to achieve the EOM in this section: in Model 1, the EOM is directly formatted to the RAPM; while in Model 2, firstly, the dynamic model of the 6RSS PM is built, and next, the two HKP constraints are modelled, to finally obtain the model of the RAPM of interest.

Model 1
In this model, this formulation is directly applied to the RAPM: at first, the mechanism is virtually cut at M i , then the dynamics of the end effector and the chains is formulated in the task space and the joint space independently by the Lagrangian formulation; next, the complete dynamic model of the mechanism without constraints is built. Finally, the model of the RAPM is acquired via the D'Alembert formulation.

The end effector
The kinetic energy of the end effector is computed as whereM EE is the inertia dyad that has been given in equation (29), and M EE1 is its4 Â 4 mass matrix that can be expressed as Its potential energy has already been derived in equation (31). Regarding these, its Lagrange function is One can find that by virtue of which the Lagrangian formulation of the end effector is d dt SeeFormulainPDF are the 4 Â 4 Coriolis/centrifugal force matrix and the4 Â 1gravitational force vector of the end effector, respectively.F EE1 (q EE ) is the4 Â 1generalized force vector corresponding to the four DOFs of the end effector and has been computed in equation (42); evidently, it is a function of q EE . Now the dynamic model of the end effector has been formulated straightforwardly in the task space.

The ith chain
The EOM of a serial chain by the Lagrangian formulation is well developed and can be found in many textbooks on robotics. As such, the details are not given and only the final EOM is written as ð Þ are themass matrix, theCoriolis/centrifugal force matrix, thegravitational force vector, and thegeneralized force vector of the ith chain, respectively, and is the actuating torque provided by the ith actuator. Note that the model is built in the joint space, and

The overall mechanism
For the complete RAPM, the model without constraints can be written as where are the generalized coordinate vector, the mass matrix, the Coriolis/centrifugal force matrix, the gravitational force vector, and the generalized force vector of the six chains without constraints, respectively.
Following the procedure of the D'Alembert formulation, one can obtain that where L q r 1 q r ;q r ;q r ð Þ ¼ M q r q r ð Þ⋅q r q EE ;q EE ;q EE ð Þ þC q r q r ; q r ð Þ⋅q r q EE ;q EE ð ÞþG q r q r ð Þ; and dq r and dq EE are the virtual displacement vectors formed by q r and q EE , respectively. Because dq EE is an independent vector, it can be found that Substituting equations (17) and (18) into equation (56) and transposing the result, it can be equivalently written as In equation (57), there are four equations and six unknowns, indicating the mechanism is under actuation redundancy.

Model 2
In this model, the EOM of the 6RSS PM without HKP constraints is written, as the function of X EE ,X_ EE ,X EE . Then the constraints from HKPs are modelled, to finally reach the model of the RAPM, where the model of the 6RSS PM is incorporated as its core.
Likewise, the Lagrange equation of the 6RSS PM free of constraints among rigid bodies is written as M q r q r ð Þ; C q r q r ; q r ð Þ; G q r q r ð Þ; F q r t ð Þ have identical structures as their counterparts in equation (54). However, it must be remembered that they are functions of X EE ,X_ EE , X EE rather than q EE ,q_ EE ,q EE . Specifically, differentiating equation (15) with respect to time once and twice, and rearranging the results for the six chains prodocue By the D'Alembert formulation, one can obtain where L q r q r ; _ q r ; q r ð Þ¼M q r q r ð Þ⋅q r X EE ; Substituting equation (59) into equation (60) and transposing the result give rise to where and J u2 can be found in equation (46). After the two HKP constraints are imposed to the end effector, its DOFs have transferred from X EE to q EE . Thus, from equation (9), one can obtain Putting it into equation (61) and "deleting" the free termdq T EE produce This set of equations is the model of the RAPM with four equations and six unknowns.
From this process, there are two things should be emphasised: the first is, though in equation (61), the terms inside the bracket are functions of q EE ,q_ EE ,q EE ultimately, the numerical values of X EE ,X_ EE ,X EE computed from equations (6), (9), and (10) are directly used in them. That is to say, the model of the 6RSS PM is used as the core of the model of the RAPM. It can alleviate the coupling in the parasitic motion variables and boost the computational efficiency significantly, and this will be illustrated in Section 6.2.1. The second is, by "deleting" the term in ((61), one can find the following equation is the EOM of the 6RSS PM. In this regard, the model of the RAPM can be formatted by left-multiplying M T J to the two sides of the EOM of the 6RSS PM, which is an identical manipulation as in Section 4.2.

Numerical computations and discussions
To verify the methods and models in this paper, the mechanism is commanded to follow a real incisor trajectory in free chewing by a healthy human subject which lasts 5 seconds. The trajectory from different perspectives is given in the first four subplots in Figure 5. The experimental setup and procedure to obtain the mastication movements from human subjects can be found in Chapter 6 of [3]. The corresponding mandibular motions in terms of the four elements of q EE are computed, using the same method in equation (38) of [1]. Next, via equation (6), the numerical values of parasitic motion variables Z andgin terms of q EE can be obtained. Finally, the first and second time-derivatives of these six motion variables are also computed via equations (9) and (10). The values of X EE ,X_ EE ,X EE versus time are provided in the following six subplots.
Four interesting features in this chewing trajectory in Figure 5 can be observed. Firstly, the sagittal view in the X S -Z S plane is like a straight line, indicating the trajectory is nearly concentrated on one plane. Secondly, from the last six subplots, during the first one second, the motions of the mandible are quite small. After that period, it performs a rhythmic chewing motion with larger amplitudes of Z andbthan those of X, Y, a and g, and the first and second time-rates of Z andbare also larger than their counterparts, indicating the chewing motions are dominated by the mouth opening and closing movements. In the third place, the lateral movement indicated by the values of Y,Y_,Y are quite small. Finally, the trajectories of a and g are nearly coincident in the time history, so are their first and second time-rates. In following this predefined trajectory, numerical values of six diagonal entries of J A in equation (26) are exhibited in Figure 6. It clearly shows that none of them is zero and as a result, equation (27) is not satisfied and an inverse singularity cannot occur. Meanwhile, det(J W ) always equals four in the time history and then equation (28) is invalid. In these regards, the planned trajectory is not passing the three sorts of singular configurations as mentioned in Section 3.3.
Concerning the 6RSS PM, its six DOFs and their first and second time-derivatives as a function of time are identical as those in the final six subplots in Figure 5. As a consequence, in the two mechanisms, the end effectors undergo identical motions in six directions, and each of the six chains in them also undergoes identical motions in terms of the numerical values of € q r i ; € q r i ; € q r i . In tracking this prescribed chewing trajectory, the 6RSS PM is at or near singular configurations using the analysis of the RAPM as in Section 3.3. Correspondingly, an experimentally measured 3D bite force in {S} on peanuts by a healthy human subject on the molars in Figure 7, is exerted onto the right molar of the end effector. In this figure, there are in total five bursts in the three components of the bite force in the timeline, and each of the burst corresponds to the jaw closing portion for each chew. The size of the z component is larger than those in the other two directions, because the force in this direction plays a more important role in the food chewing process.

Generalized momentum approach
In both Models 1 and 2 by this approach, there are four equations and six unknowns in equations (43) and (46), indicating there are infinite solutions for the torque distribution theoretically. A closed-form solution for the actuating torques can be obtained simply by the pseudoinverse of J T u1 and M T J ⋅J T u2 , respectively, as ð65Þ which corresponds to the minimum of actuating torques. From model 1, the input torquestversus time along the chewing trajectory are exhibited in the first column of Figure 8, and the torque differences between Model 1 and Model 2 in this approach are given in the second column. From the first column, one can see that the five bursts produced by each actuator are tightly following those in the bite force in the time history in Figure 7, rather than those in the free chewing motion in Figure 5. It means the torque cost required by the bite force is larger than that by the inertia and gravitational forces of the mechanism itself. More importantly, in the second column the differences are quite minor, thus, the two models are equivalent in terms of the numerical results.

Lagrange-D'Alembert formulation
Via the two models by the Lagrange-D'Alembert principle in equations (57) and (64), the actuating torques can also be minimised using the pseudo-inverse solution, respectively, as The actuating torques from Model 1 and the differences between the actuating torques from the two models under this formulation are equivalent to those in Figure 8, respectively. As a consequence, they are not exhibited for the sake of brevity. Furthermore, the equivalent results demonstrate the correctness and accuracy of the developed four models by the two methods.

Computational demands
The computational efficiency of inverse dynamics is crucial to the real-time application of the RAPM, and it is an important factor in judging the performance of the model. To assess the suitability of the built models for real-time control, reliable quantitative measurement of the computational load is critical. The time consumption is summarised in Table 1, where GM and L-D are short for the generalized momentum approach and the Lagrange-D'Alembert formulation, respectively, and the number 1 and 2 are meaning the first and the second model in the corresponding method, respectively. The procedures under each model are all divided into symbolic and numeric computations, and have been implemented in programs written in Matlab, using an Intel(R) Core(TM) i7-8700K CPU@3.70GHz and 64GB of RAM. By the symbolic EOMs, it is rather convenient to obtain the mathematical functions of the Jacobian matrix, the kinetic and gravitational components of the generalized forces, etc. in these methods, and they can be called in the numerical computations.  In the dynamic models where the two methods are directly utilised to the RAPM, i.e., in the first two columns in Table 1, the model under the newly developed generalized momentum approach is more computationally efficient than that under the classical Lagrange-D'Alembert formulation by about 40% and 43.6% in the symbolic and numeric computations, respectively. The reason for the superiority of the generalized momentum approach roots in the high efficiency in deriving the Coriolis/centrifugal force matrix by the differentiation of the generalized momentum, as stated in [16,17]. In the two models where the EOM of the 6RSS PM is incorporated as the core of the model of the RAPM, namely, from the third and fourth columns, one can also reach the conclusion that the generalized momentum approach is faster than the Lagrange-D'Alembert formulation.
From the comparisons between the first four columns in Table 1, it can be seen using the dynamics of the 6RSS PM can boost the efficiency to a huge degree. The total time consumption by the two methods in Model 2 is only about 22.5% and 16.2% of that in Model 1, respectively. The main reason for this considerable improvement is that the resulting complexity of Model 1 by each method is a consequence of using the complex symbolic expressions of Z and g from equation (6), and the Jacobian matrix M J from equation (9), for q EE ,q_ EE ,q EE are used directly. By contrast, in Model 2, even though the matrices and vectors are ultimately functions of q EE ,q_ EE ,q EE , they can be computed by the numerical values of X EE ,X_ EE ,X EE directly, which can be available from q EE ,q_ EE ,q EE at the beginning of the computation, and then are fed into equations (46) and (63). In this manner, the strong coupling between the two parasitic motion variables Z and g and the four DOFs q EE is completely alleviated. In other words, Model 2 is free of the complex symbolic computations brought by Z,g, and M J . As such, the symbolic mathematical functions, which are called and used in numerical computations, are less complex than those in Model 1, offering simpler possible computational algorithms. In this regard, a shortcut that can greatly boost the computational efficiency is found, by incorporating the dynamics model of the 6RSS PM as the core of that of the RAPM. On this basis, a left-multiplication ofis simply performed to the two sides of EOMs of the 6RSS PM as in equations (47) and (64). Indeed, this shortcut is inspired intuitively by the comparison between the RAPM and the 6RSS PM: the former is generated by the two HKP constraints onto the end effector of the latter.

Influence of HKP constraints
It is remembered that the RAPM is built by imposing two HKP constraints onto the end effector of the 6RSS PM. Thereby, the dynamic model of the 6RSS PM free of HKPs is also built via the foregoing two methodologies, as far as the role of the HKP constraints in the numerical results and the computational cost is concerned. equations (47) and  (64) are the EOMs of the 6RSS PM by these two methods, and the variables inside them are functions of X EE ,X_ EE , X EE , which are identical as in the final six subplots of Figure 5. By these two approaches, its inverse dynamics problem is reduced to solving a system of six linear algebraic equations in six unknowns. The input torques can be uniquely determined and there is no optimisation as in the RAPM. Their profiles in the time history are similar to those in the first column of Figure 6, and the only difference is the magnitudes of the torques in the 6RSS PM are larger than those in the RAPM. Thus, they are not exhibited to save pages. Specifically, the contributions the total torques by the end effector, the six couplers, and the six cranks, respectively in the RAPM and the 6RSS PM are to be discovered and compared. in the RAPM, according to the EOM by the generalized momentum approach, the torques required by the abovementioned three parts are given in the three columns of Figure 9, respectively. Clearly, the bursts in each actuator are tightly tracking those in the free-chewing trajectory as shown in the final six subplots in Figure 5. In addition, in one subplot, the torques offered by each pair of actuators are more or less equivalent. A very large portion of the total torques is required by the end effector, but those by the couplers and the cranks are quite minor. From the comparison between Figure 9 and the first column of Figure 8, it is obvious that the torque cost by the motion of the RAPM itself is quite small, but the torque required to sustain the bite force occupies a large portion. In the 6RSS PM, the torque distributions by the three parts are experiencing a similar profile as those in Figure 9 but with larger magnitudes, thus they are not exhibited to save pages.
To quantitatively justify and compare the influence by HKPs to the actuating torques, in the RAPM and the 6RSS PM, four indices are set about the numerical values of total actuating torques, and torques required by the aforemen-tioned three parts as and N is the number of sampling points in the time history. Obviously, the total actuating torques t can also be obtained from Model 2 by the generalized momentum approach, or from the two models by the Lagrange-  ; t GS 1 6 can also be computed from Model 2 by the generalized momentum approach; however, by the Lagrange-D'Alembert formulation, the contributions by each of the three parts to the total actuating torques are not as easily available as from the generalized momentum approach.
The numerical values of these indices of the RAPM and the 6RSS PM are given in the first and the second line of Table 2, respectively. The value of F is always 0.0207N.m from the RAPM under the four models; while in the 6RSS PM, F equals 0.0723N.m from two methods and is larger than that by the RAPM over three times. It means that redundant actuation can minimise the input torques, which is a well-established and well-known opinion and has been proved in a number of publications, even though the way to generate redundant actuations in the mechanism is rather different from the two primary methods in [20]: they are caused by transferring two DOFs into parasitic motions under the HKP constraints. In other words, smaller torques are devoted to the remaining four DOFs, even though they can give rise to two parasitic motion variables. Besides, the torques required to maintain the motion of the aforementioned three parts increase by 72.32%, 52.94%, and 9.09% from the RAPM to the 6RSS PM, respectively. Thereby, both the proportion in the end effector and the torques required by it are the largest. In these regards, an important conclusion can be reached: even though the end effector undergoes identical motions in six directions in terms of displacements, velocities, and accelerations in both the RAPM and the 6RSS PM, the required torques are evidently smaller in the parasitic motion variables Z and gin the RAPM than those by the corresponding two DOFs in the 6RSS PM.
The computational time for the 6RSS PM under the two methods is given in the last two columns of Table 1. It is noted that the symbolic computational time is equivalent to that in Model 2 by each method in the RAPM, respectively, because one common file in Matlab for the symbolic computation is shared by the RAPM in Model 2 and the 6RSS PM. The numerical computational burden of the generalized momentum approach is also lower than that of the Lagrange-D'Alembert formulation, which is an identical conclusion found in the RAPM. Meanwhile, as expected, the time used in the two methods in the 6RSS PM is comparable to that in Model 2 by the two methods, respectively, because the core of the model of the 6RSS PM is used in Model 2, and the difference is only a numerical left-multiplication of M T J . By comparing the time consumption between the 6RSS PM and Model 1 from each method in the RAPM, where both the two methods are directly used to the two mechanisms, it is found that the two HKP constraints greatly increase the computational complexity. The reason why the RAPM consumes considerably more time is that the complex expressions of parasitic motion variables Z andg, and the Jacobian matrix M J , etc. caused by the HKP constraints greatly increase the demands in the symbolic computations. On this basis, more complicated numerical computations are processed. By contrast, in the 6RSS PM, Z andgare directly used as DOFs which are much more straightforward, and there is no Jacobian matrix M J .

Conclusion
The inverse dynamics of a spatial RAPM constrained by two HKPs was solved systematically via the generalized momentum approach and the Lagrange-D'Alembert formulation. The methods and models were thoroughly studied and compared. Generally, studies of inverse dynamics concentrate on RAPMs with only lower kinematic pairs. Thereby, the scientific contribution of this paper is the deep study of inverse dynamics in a spatial RAPM with both lower kinematic pairs and HKPs in terms of revealing the nature of HKP constraints. Specifically, the following conclusions can be drawn: -The two models in each method give rise to identical numerical results; however, by incorporating the dynamic model of the 6RSS PM as the core of that of the RAPM in Model 2 in each method, the computational efficiency has been enhanced remarkably as fast as that of the 6RSS PM without any simplification or loss of accuracy. -The two methods are equivalent in terms of the numerical results; nonetheless, the computational demands are greatly lower in the newly developed generalized momentum approach than in the classical Lagrange-D'Alembert formulation. -By comparing the computational requirements between the RAPM and the 6RSS PM, it is discovered that the HKP constraints greatly raise the computational complexity. -Compared to the 6RSS PM with six DOFs, in the RAPM smaller torques are devoted to its four DOFs and two parasitic motion variables. In other words, the torques required by the parasitic motion variables Z and gin the RAPM are smaller than those by the corresponding DOFs in the 6RSS PM.
The models in this paper can provide a great potential to facilitate the real-time motion and/or force control of the RAPM. It is also expected to straightforwardly extend the presented process to a class of RAPMs whose end effectors are constrained by the base directly.
In future work, the new form of Lagrange-D'Alembert formulation in [27] for dynamic modelling of PMs will be employed to this RAPM.

Nomenclature
To comprehend the paper conveniently, a list of symbols together with their definitions is provided in the following