© C. Cheng et al., Published by EDP Sciences 2023

1 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 (revolute-spherical-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–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 non-redundantly 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.

2 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_iS_i (i =1,…,6) is driven by a base-mounted rotary actuator, whose shaft connects it with a rotational joint at G_i, and a coupler S_iM_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. .axis is directed from G_i to S_i, the axis runs through the driving shaft of the actuator, and the axis completes the frame, obeying the right-hand rule. A frame {N_i} is attached at the mass centre E_i of S_iM_i to describe its motions with respect to {S}. The axis points from S_i to M_i, the axis is parallel to the cross product of two unit vectors defined along the and X_S axes, and the 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.

Fig. 1 A schematic view of the RAPM constrained by two point contact HKPs, where ① right condyle ball, ② left condyle ball, ③ articular surface of right TMJ, ④articular surface of left TMJ [21].

Fig. 2 The RAPM: (a) CAD model, (b) magnification of the right HKP, (c) prototype of the HKP [2], (d) end effector.

Fig. 3 Prototype of the mechanism [2].

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.

3.1 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–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):

(1)

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

(2)

where denotes the 3◊1 position vector of O_M in {S}, is the rotation matrix from {S} to {M}, R_X (α) , R_Y (β), andR_Z (γ)are three rotation matrices about the X_M, Y_M, and Z_M axes by three Bryant anglesα, β, andγ, 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

(3)

From equation (2), one can obtain

(4)

where is the ith(i = 1,3) row of . Putting equation (4) into equation (1) produces

(5)

Because of the left-right symmetry ofMO_MT_LandMO_MT_Rin {M}, a summation and a subtraction of the two equations in equation (5) sidewise yield

(6)

where MO_MT_L (1)and MO_MT_L (3)are the first and third terms ofMO_MT_L, respectively. From these computations, it is found Z and γ 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

(7)

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/γand 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

(8)

where

0₃ is 3 × 3 the zero matrix, and M_J denotes the 6 × 4 Jacobian matrix between X_EE and q_EE, namely

(9)

From it, the second time-derivative of X_EE can be derived as

(10)

The translational velocity of the end effector is expressed by

(11)

where E₃ is the 3 × 3 identity matrix.

3.2 Kinematics of the ith chain

The inverse kinematics of the mechanism, i.e., that 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_iM_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 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 β_i and γ_i around the and axes, respectively, are used to express the rotation of S_iM_i in {S} in terms of the two rotational matrices R_Y (β_i) and R_Z (γ_i). Thereafter, the coordinate vector of the coupler can be expressed as

(12)

where is the orientation of S_iM_i in {S} at the initial configuration of the mechanism at hand, and∥S_iM_i∥is the length of S_iM_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:

(13)

where and is the coordinate vector of M_i in {M}; O_SG_iis the constant position vector of G_i in {S}; and , in which is the orientation of G_iS_i in {S} at the initial configuration of the mechanism,R_Z (θ_i) is the rotation matrix about the axis by θ_i, and ∥G_iS_i∥is the length of G_iS_i.

Substituting equation (12) and θ_i = θ_i (q_EE) into equation (13) produces

(14)

where . For the sake of convenience, a3 × 1 generalized vector is defined as , which consists of the joint space of the ith chain and completely specifies its configuration. It is highlighted that is the function of q_EE, i.e., .

To derive the relationship between the first time-derivatives of andq_EE, equation (13) can be rewritten as

(15)

The two sides can be expressed by and q_EE, respectively. The first time-derivative of equation (15) yields

(16)

where

Moreover, one can find that

(17)

where , and the second time-derivative of is given as

(18)

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 S_iM_i is

(19)

where . Substituting equation (17) into equation(19) renders

(20)

The coordinate vector of the mass centre E_i in {S} is

(21)

The translational velocity of E_i can then be found by the differentiation of equation (21) with respect to time as

(22)

Upon substitution of equation (17) into equation (22), it yields

(23)

As far as the crank G_iS_i is concerned, its mass centre G_i is passing through the rotational shaft of the actuator, as such . Its rotational velocity in {S} is

(24)

where is the third column of .

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.

Fig. 4 The structure of dynamic modelling by the two methods in Sections 4 and 5.

3.3 Singularity analysis

To analyse the singular configurations of the mechanism, a 6◊1 constraint vector is defined as

(25)

Its first time-derivative is

(26)

where

and their sizes are 6◊6 and 6◊4, respectively. It is also noted that the ith(i = 1,…,6) row of L (θ, q_EE) is the function of θ_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.,

(27)

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.

2. A direct kinematic singularity occurs when J_W is rank-deficient, i.e.,

(28)

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.

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

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

4.1.1 The end effector

On the one hand, the generalized momentum of the end effector is computed as

(29)

where 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 , and the size of q_EE1 is 4 × 1.

The kinetic component of the generalized force acting on the end effector can be computed as

(30)

On the other hand, the potential energy of the end effector is

(31)

where g = 9800 m/s² is the gravitational acceleration. The gravitational component of the generalized force is

(32)

4.1.2 The ith chain

As far as the coupler S_iM_i of the ith kinematic chain is concerned, by virtue of equations (20) and (23), its generalized momentum is computed as

(33)

where is the 6 × 6 inertia dyad of S_iM_i, and are the mass and the inertia matrix with respect to O_M and expressed in {S}, respectively, and . The kinetic component of the generalized force acting on S_iM_i is computed as

(34)

The potential energy of S_iM_i is

See Formula in PDF

(35)

whereO_SE_i (3) is the third term of O_SE_i. The gravitational component of the generalized force is

(36)

For the crank G_iS_i, its generalized momentum is

(37)

where is a constant6 × 1vector, and the scalar is the rotational inertia of G_iS_i around the axis of frame{C_i}. As such, its kinetic component of the generalized force is

(38)

Because G_i is a fixed point in {S}, its potential energy is constant and computed as

(39)

in whichO_SG_i (3) is the third term of O_SG_i. Thereby, its gravitational component of the generalized force is

(40)

4.1.3 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

(41)

whereτis the6 × 1actuating torque vector by the six actuators, and is the6 × 1wrench vector by the bite force F_B. Substituting equations (8), (11), and (17) into equation (41), and rearranging the term by give rise to

(42)

4526+++++++where is the generalized force vector by F_B, is the 6 × 4 Jacobian matrix mapping q˙_EE into θ˙, i.e., θ˙=J_θ1⋅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

(43)

which is the inverse dynamic model of the RAPM and it contains four equations and six unknowns.

4.2 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 rather than . 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:

(44)

where , , , , and 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

(45)

where is the generalized force vector by the bite force F_B, J_θ2 is the 6 × 6 Jacobian matrix mapping X˙_EE to θ˙, namely, θ˙=J_θ2⋅X˙_EE. Because q˙T_EE is free to vary, it can be rewritten as

(46)

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 , they are directly expressed and computed by . 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

(47)

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-multiplying at the two sides of that of the 6RSS PM, which is used as the core of the dynamic model of the RAPM.

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

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

5.1.1 The end effector

The kinetic energy of the end effector is computed as

(48)

whereM_EEis the inertia dyad that has been given in equation (29), and M_EE1 is its4 × 4 mass matrix that can be expressed as

(49)

Its potential energy has already been derived in equation (31). Regarding these, its Lagrange function is

(50)

One can find that

(51)

 by virtue of which the Lagrangian formulation of the end effector is

(52)

where

See Formula in PDF

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.

5.1.2 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

(53)

where and 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 .

5.1.3 The overall mechanism

For the complete RAPM, the model without constraints can be written as

(54)

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

(55)

 where

and and are the virtual displacement vectors formed by q_r and q_EE, respectively. Because δq_EEis an independent vector, it can be found that

(56)

Substituting equations (17) and (18) into equation (56) and transposing the result, it can be equivalently written as

(57)

where

In equation (57), there are four equations and six unknowns, indicating the mechanism is under actuation redundancy.

5.2 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

(58)

where and 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

(59)

By the D'Alembert formulation, one can obtain

(60)

where

Substituting equation (59) into equation (60) and transposing the result give rise to

(61)

where

and J_θ2 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

(62)

Putting it into equation (61) and “deleting” the free term produce

(63)

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

(64)

is the EOM of the 6RSS PM. In this regard, the model of the RAPM can be formatted by left-multiplying to the two sides of the EOM of the 6RSS PM, which is an identical manipulation as in Section 4.2.

6 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 andγin 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 andβthan those of X, Y, α and γ, and the first and second time-rates of Z andβare 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 α and γ 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 . 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.

Fig. 5 Motions of the end effector in 3D space and time history [22].

Fig. 6 Numerical values of the diagonal elements of JA.

Fig. 7 3D bite force profiles on peanuts [3].

6.1 Solutions of two methods

6.1.1 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 pseudo-inverse of and , respectively, as

(65)

which corresponds to the minimum of actuating torques. From model 1, the input torquesτversus 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.

Fig. 8 Torques from Model 1 and differences between the two models by the generalized momentum approach.

6.1.2 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

(66)

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.

6.2 Discussion

6.2.1 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 γ 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 γ and the four DOFs q_EE is completely alleviated. In other words, Model 2 is free of the complex symbolic computations brought by Z,γ, 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.

6.2.2 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 aforementioned three parts as

(67)

where

and N is the number of sampling points in the time history. Obviously, the total actuating torques τ can also be obtained from Model 2 by the generalized momentum approach, or from the two models by the Lagrange-D'Alembert formulation. 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 γin 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 .

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 andγ, 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 andγare directly used as DOFs which are much more straightforward, and there is no Jacobian matrix M_J.

Fig. 9 Torques required by the end effector, the six cranks, and the six couplers in the RAPM.

7 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 γin 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 table at their appearance order.

{S} Inertia frame: {M} Frame established at the mass centre O_M of the end effector

G_iS_i(i=1,…,6): The ith crank

S_iM_i(i=1,…,6): The ith coupler

{C_i}: Frame established at the fixed point G_i of the ith crank

{N_i}: Frame attached at the mass centre E_i of S_iM_i

T_j (j=L, R): Left and right condyle ball centres

O_SO_M: 3◊1 position vector of O_M in {S}

: Rotation matrix from {S} to {M}

R_X (α),R_Y (β), andR_Z (γ): Rotation matrices about X_M, Y_M, and Z_M axes by Bryant anglesα, β, andγ, respectively

X_EE: 6◊1 vector grouping six motion variables of the end effector

X˙_EE,X¨_EE: First and second time-derivatives of, respectively

p_i(i=1,…,6): Parameters describing the shape of surfaces where condyle ball centres T_L and T_R slide on

X,Y, and Z : Parameters describing translations of the end effector

α,β, and γ: Bryant angles describing rotations of the end effector

q_EE: 4◊1 vector grouping four DOFs of the end effector

q˙_EE,q¨_EE: First and second time-derivatives of q_EE, respectively

ω_EE: Angular velocity of the end effector

M_J: 6◊4 Jacobian matrix between X_EE and q_EE

0₃: 3◊3 zero matrix

E₃: 3◊3 identity matrix

θ: 6◊1 vector containing the displacements of six active joints

β_iandγ_i: Euler angles expressing the rotation of the ith coupler

∥S_iM_i∥: Length of S_iM_i

O_MM_iandMO_MM_i: Vectors ofO_MM_iin {S} and {M}, respectively

O_SG_i: 3◊1 constant position vector of G_i in {S}

: Rotation matrix from {S} to{C_i}at the initial configuration of the mechanism

R_Z (θ_i): Rotation matrix about the axis byθ_i

∥G_iS_i∥: Length ofG_iS_i

: 3◊1 generalized vector specifying the configuration of the ith kinematic chain

M_1i: Jacobian matrix of O_SG_i + G_iS_i + S_iM_i with respect to

M_2i: Jacobian matrix ofO_SO_M + O_MM_i with respect toq_EE

: Angular velocity of S_iM_i

: Translational velocity of mass centre E_i of S_iM_i

: Jacobian matrix ofO_SE_iwith respect to

: Angular velocity of G_iS_i

L (θ, q_EE): 6◊1 constraint vector to analyse singularities

J_A: Jacobian matrix betweenL (θ, q_EE)and θ

J_W: Jacobian matrix betweenL (θ, q_EE)andq_EE

q_EE1: Generalized momentum of the end effector in Model 1 under the generalized momentum approach

M_EE: 6◊6 inertia dyad of the end effector

m_EE and I_EE: Mass and inertia matrix of the end effector with respect to O_M and expressed in {S}, respectively

: Kinetic component of the generalized force acting on the end effector in Model 1 under the generalized momentum approach

P_EE1: Potential energy of the end effector in Model 1 under the generalized momentum approach

g: Gravitational acceleration

: Gravitational component of the generalized force acting on the end effector in Model 1 under the generalized momentum approach

: Generalized momentum of S_iM_i in Model 1 under the generalized momentum approach

: 6◊6 inertia dyad of S_iM_i

and : Mass and inertia matrix of S_iM_i with respect to O_M and expressed in {S}, respectively

: Kinetic component of the generalized force acting on S_iM_i in Model 1 under the generalized momentum approach

: Potential energy of S_iM_i in Model 1 under the generalized momentum approach

: Gravitational component of the generalized force in Model 1 under the generalized momentum approach

and : Mass and rotational inertia of G_iS_i around the axis of frame{C_i}, respectively

: Gravitational component of the generalized force of G_iS_i in Model 1 under the generalized momentum approach

τ: 6◊1 actuating torque vector by the six actuators

: 6◊1 wrench vector by the bite force F_B

 F_B: 3◊1 bite force vector

J_θ1: 6◊4 Jacobian matrix mapping q˙_EE into θ˙

F_EE1: 4◊1 generalized force vector by F_B in Model 1 under the generalized momentum approach

and : Kinetic and gravitational components of the generalized force acting on the end effector, respectively, in Model 2 under the generalized momentum approach

and : Kinetic and gravitational components of the generalized force acting on the ith coupler, respectively, in Model 2 under the generalized momentum approach

and : Kinetic and gravitational components of the generalized force acting on the ith crank, respectively, in Model 2 under the generalized momentum approach

J_θ2: 6◊6 Jacobian matrix mapping X˙_EE into θ˙

T_EE1: Kinetic energy of the end effector in Model 1 under the Lagrange-D'Alembert formulation

L_EE1: Lagrange function of the end effector in Model 1 under the Lagrange-D'Alembert formulation

C_EE1(q_EE,q˙_EE): 4◊4 Coriolis/centrifugal force matrix, gravitational force vector, and generalized force vector of the end effector, respectively, in Model 1 under the Lagrange-D'Alembert formulation

: 3◊3 mass matrix, 3◊3 Coriolis/centrifugal force matrix, 3◊1gravitational force vector, and 3◊1 generalized force vector of the ith chain, respectively, in Model 1 under the Lagrange-D'Alembert formulation

τ_i: Actuating torque provided by the ith actuator

δq_r, δq_EE, δX_EE: Virtual displacement vectors of q_r, q_EE, and X_EE, respectively

F, F_EE, , : Indices about the numerical values of total actuating torques, and torques required by the end effector, six couplers, and six cranks, respectively

HKP: Higher kinematic pair

RAPM: Redundantly actuated parallel mechanism

3D: Three-dimensional

RSS: Revolute-spherical-spherical

TMJ: Temporomandibular joint

EOM: Equations of motion

DOF: Degree of freedom

PUS: Prismatic-universal-spherical

Conflict of interest

The authors declare that they have no conflict of interest.

Data availability

The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Data availability

This work was supported by the Natural Science Foundation of Shaanxi Province (Grant No. 2019JQ-467).

Acknowledgements

The authors also sincerely appreciate the valuable time, suggestions and comments from the editors and the anonymous reviewers for improving the quality of the paper.

References

All Tables

All Figures

	Fig. 1 A schematic view of the RAPM constrained by two point contact HKPs, where ① right condyle ball, ② left condyle ball, ③ articular surface of right TMJ, ④articular surface of left TMJ [21].
In the text

	Fig. 2 The RAPM: (a) CAD model, (b) magnification of the right HKP, (c) prototype of the HKP [2], (d) end effector.
In the text

	Fig. 3 Prototype of the mechanism [2].
In the text

	Fig. 4 The structure of dynamic modelling by the two methods in Sections 4 and 5.
In the text

	Fig. 5 Motions of the end effector in 3D space and time history [22].
In the text

	Fig. 6 Numerical values of the diagonal elements of JA.
In the text

	Fig. 7 3D bite force profiles on peanuts [3].
In the text

	Fig. 8 Torques from Model 1 and differences between the two models by the generalized momentum approach.
In the text

	Fig. 9 Torques required by the end effector, the six cranks, and the six couplers in the RAPM.
In the text