A Redundantly Actuated Chewing Robot Based on Human Musculoskeletal Biomechanics: Differential Kinematics, Stiffness Analysis, Driving Force Optimization and Experiment

: Human masticatory system exhibits optimal stiffness, energy efﬁciency and chewing forces needed for the food breakdown due to its unique musculoskeletal actuation redundancy. We have proposed a 6PUS-2HKP (6 prismatic-universal-spherical chains, 2 higher kinematic pairs) redundantly actuated parallel robot (RAPR) based on its musculoskeletal biomechanics. This paper studies the stiffness and optimization of driving force of the bio-inspired redundantly actuated chewing robot. To understand the effect of the point-contact HKP acting on the RAPR performance, the stiffness of the RAPR is estimated based on the derived dimensionally homogeneous Jacobian matrix. In analyzing the inﬂuence of the HKP on robot dynamics, the driving forces of six prismatic joints are optimized by adopting the pseudo-inverse optimization method. Numerical results show that the 6PUS-2HKP RAPR has better stiffness performance and more homogenous driving power than its non-redundant 6-PUS counterpart, verifying the beneﬁts that the point-contact HKP brings to the RAPR. Experiments are carried out to measure the temporomandibular joint (TMJ) force and the occlusal force that the robot can generate. The relationship between these two forces in a typical chewing movement is studied. The simulation and experimental results reveal that the existence of TMJs in human masticatory system can provide more homogenous and more efﬁcient chewing force transmission.


Introduction
A chewing robot in an oral context has been found in many applications regarding the emotion and expression of social robotics, the analysis of the chewing process, dental and joint prosthetic material tests and food science [1][2][3][4][5]. In order to design a high anthropomorphic chewing robot, the working principle of human masticatory system should be studied. Simultaneously, in turn, the chewing robot can provide a useful tool for studying biomechanical behavior of human masticatory system [6]. The masticatory system itself is actuation redundant as the mandible is driven by greater groups of muscles than required and is constrained by two temporomandibular joints (TMJs). The human mandible is attached to the skull by muscles and guided by passive structures such as posterior teeth and TMJ at each side of the jaw [7]. There have been a variety of chewing machines and devices available for replicating the human masticatory movements and occlusal forces. However, the functions of the TMJ and redundant actuation features of the masticatory system have not been given sufficient consideration during the design of chewing robots. However, the TMJ not only affects the mandibular movement, but also affects the force of the mandible. Ignoring the TMJ will result in poor bio-imitability of chewing robot.
In the area of human-machine interaction, most of the researches [8][9][10] used one single servo motor to raise and lower the lower jaw. The problem with these robots is that they lack sufficient DOFs for reproducing complex human jaw gestures. Only a few chewing robots designed for robotic head have multiple DOFs, like Flores's jaw [1]. However, the compound motions of this chewing robot were realized by simple combination of drive units. In the area of rehabilitation, Waseda-Yamanashi (WY) series of robots can imitate the movement of doctor's hands to train the patients who suffer disorders in opening/closing the mouth [3]. The WY robot showed good performance in dental training, but the driving structure is not completely consistent with human masticatory muscles. The Waseda Jaw (WJ) robot was developed to cooperate with a WY robot for dental training and was also be used to explore the physical phenomena and resistance force during jaw training of patients [11]. However, the WJ robot used eleven artificial muscle actuators (AMA) which can only generate contraction forces. Aimed at quantifying food texture changes during chewing, a number of chewing robots have been developed [12][13][14]. However, either the TMJ was not considered or the relationship between occlusal force and TMJ force were ignored in the design of these robots. An in vitro wear simulator based on Stewart platform was developed to replicate dental wear formations on dental components, such as crowns, bridges or a full set of teeth [15]. Typical movements and occlusal forces can be realized on this simulator, but the influence of TMJ force to occlusal force was ignored. The wear behavior of condyle/fossa pairs in TMJ was tested in a mandibular movement simulator for an equivalent of two years of clinical use [16]. Nevertheless, the TMJ articulation force affected by occlusal force during human chewing movement cannot be replicated on this simulator. Moreover, as far as we know, there is no such mechanical device, which can reproduce the musculoskeletal actuation redundancy of human masticatory system and no literature has yet studied the relationship between occlusal force and TMJ force on a robotic simulator.
Basically, actuation redundancy can be achieved in three ways: Actuating the passive joints, adding additional kinematic chains, or hybrid of the first two types [17,18]. That is, adding more actuators than necessary to control the robot for a given task [19]. Redundant actuation has been proved to be an effective method for eliminating singularities, improving manipulability, stiffness and torque distribution in parallel robots (PRs) [20][21][22]. Liu et al. [23] investigated a two degrees of freedom (DOFs) planar mechanism and analyzed the distribution of singularity of this mechanism. Zhu and Dou studied the kinematic model and Jacobian matrix of a 2-DOF 3-RRR planar PRs [24]. Marquet et al. [25] presented a 3-DOF and four actuators redundant parallel mechanism, ARCHI, which has three PUS chains and a PRR chain, improving the mechanical strength, stiffness and singularities. Wu et al. [26,27] compared the performance of a 3-DOF planar PR with its corresponding non-redundant one, showing the redundantly actuated one has better dexterity, littler singular configurations and higher stiffness. It should be noted that most of studies are about planar type redundantly actuated parallel robots (RAPRs). So far, only a few studies are available about spatial ones. Kim et al. [28] added two redundant actuators to the Eclipse robot eliminating the singularities within the workspace. Zhao et al. [29] discussed the dynamic performance comparison between the 8-PSS redundant parallel manipulator and its counterpart 6-PSS parallel manipulator. It can be seen that these redundantly actuated robots are all achieved by the above three ways, but the redundant actuation of the robot studied in this paper is achieved by reducing the DOFs of the end-effector by inclusion of additional constraints.
Bionic robot has been a hot research area for many years. To some extent, human body is a very successful evolved cyborg. Literatures about robots which are bio-inspired by human musculoskeletal structure have attracted much attention from the academic community, including mechanical design, stiffness analysis, dynamics modeling, and energy-efficiency analysis, etc. [30][31][32]. For human masticatory system, the maximum bite force that human being is able to produce with one's maximum strength can be up to 50-60 kgf [33], which is almost one's own weight. While, the size of lower jaw and TMJ is very small to burden such a heavy weight compared with the leg and knee joint. Although the human masticatory system is relatively small, it has good performance in stiffness, force output and energy saving. It seems the musculoskeletal actuation redundancy characteristics and the existence of TMJs play significant roles in masticatory system. In the fields, such as man-machine interaction, implantology, pathology and masticatory biomechanics research [34,35], it is very necessary to understand how this nonlinear redundantly actuated system works and what role of TMJs play in human masticatory system. In addition, the relation between occlusal force and TMJ force during chewing is also a point of interest in academia. The hypothesis that the TMJ force increases with the rise of occlusal force has been widely accepted in recent years [36]. In terms of the principle of operation of the whole redundantly actuated masticatory system, there is still much work to be done. Computational and mathematical modeling, such as static model [37,38], dynamic model [39], and finite-element model [40,41], were utilized in studying the masticatory system these days. Various objective functions for solving this actuation redundant system were evaluated: Minimization of joint loads, minimization of muscle effort [42] and minimization the sum of squared muscle activations [43]. These computational models are beneficial to the design of chewing robot. Meanwhile, the chewing robot can trace their foundation back to the physiological structure and computational model of the masticatory system and in turn, modeling, analysis and experiments of chewing robot may also verify the hypothesis of computational models.
In order to study the good performance of optimal stiffness and energy efficiency in redundantly actuated human masticatory system from a robotics perspective. Section 2 introduces the bio-inspired 6PUS-2HKP (6 prismatic-universal-spherical chains, 2 higher kinematic pairs) redundantly actuated chewing robot which was proposed based on masticatory biomechanics. This robot can not only mimic the musculoskeletal actuation redundancy of human masticatory system, but also reproduce occlusal force and TMJ force more realistically by introducing two HKPs. In the modeling and analysis of this nonlinear system, the differential kinematics model of this RAPR is derived in Section 3. For this RAPR with mixed DOFs, dimensionally homogeneous non-square Jacobian matrix formulation method based on three end-effector points are proposed. In order to study the significance of TMJs in redundantly actuated human masticatory system, the stiffness characteristic based on Jacobian matrix is analyzed and the optimization of driving forces is investigated in Section 4. The simulation results of the 6PUS-2HKP RAPR and its counterpart 6-PUS PR are compared to prove the RAPR with point contact constraints has better performance in stiffness and energy saving. Section 5 presents the chewing experiment of the RAPR. The chewing forces of the robot are verified and the relation between occlusal force and the corresponding TMJ force are studied via chewing simulated foods.

The Masticatory System
The masticatory system consists of an upper and a lower jaw. The lower jaw is attached to the skull pivoted at the condyle via two TMJs. There are a large number of muscles of various shapes and sizes driving the lower jaw to perform chewing movement, among which the main muscles include masseter, temporalis and lateral pterygoid muscles [7]. Each of them is made up of multiple groups of muscle fibers, i.e., the masseter muscle is composed of MAS-S and MAS-P, as shown in Figure 1a. Different muscle with its line-ofaction contributes to masticatory movement in a unique manner. The muscle insertion coordinates and angle of line-of-action of each muscle group are estimated by referring to [14,44]. The resultant force line-of-action of each muscle group (MAS, TEM and LPT) is indicated by the purple arrow in Figure 1a. The TMJ is a joint between the temporal bone and the condyle of the lower jaw. As the joint is treated as a rigid body, the disk and ligament are ignored. The condyle moves anterior and inferior along the mandibular fossa, enabling the mandible to accomplish the translation movement, as shown in Figure 1b. The condyle also revolves on the fossa surface to complete the rotation movement. In general, mastication movement involves both translation and rotation movements at the same time. Different with other joints of human body, like hip or knee joint, these two TMJs are correlative, which are guided by two articular surfaces linked by the rigid lower jaw (see Figure 1). referring to [14,44]. The resultant force line-of-action of each muscle group (MAS, TEM and LPT) is indicated by the purple arrow in Figure 1a. The TMJ is a joint between the temporal bone and the condyle of the lower jaw. As the joint is treated as a rigid body, the disk and ligament are ignored. The condyle moves anterior and inferior along the mandibular fossa, enabling the mandible to accomplish the translation movement, as shown in Figure 1b. The condyle also revolves on the fossa surface to complete the rotation movement. In general, mastication movement involves both translation and rotation movements at the same time. Different with other joints of human body, like hip or knee joint, these two TMJs are correlative, which are guided by two articular surfaces linked by the rigid lower jaw (see Figure 1).  As there are more muscles driving the lower jaw than required, and two TMJs constraining the motion of the condyle, the human masticatory system itself is mechanically actuation redundant [45]. There is a very short distance in the direction which the condyle moves perpendicular to the articular surface. The articular surface of the TMJ is assumed to be rigid and the condyle keeps contact with the surface all the time. The DOFs of the lower jaw is four.

The Redundantly Actuated Chewing Robot
Inspired by human masticatory system, a 6PUS-2HKP RAPR has been proposed [46]. The major challenge in designing the chewing robot deals with the TMJs and the masticatory muscles. In human TMJ, the two condyles are constrained in moving along the fossa surfaces, as shown in Figure 2a. Although a ball-socket joint has three rotational DOFs, it does not allow for translations. While, the condyle actually moves almost freely along the articular surface of the temporal bone in its articular capsule. Therefore, two point-contact HKPs are designed to model the articular contact of human TMJs. The HKP consists of a condylar rod and a fossa structure, as shown in Figure 2b,c. The condylar rods keep contact with the fossa articular surfaces.
Six PUS chains are arranged in parallel corresponding to the resultant line-of-action of six muscles (MAS, TEM and LPT in Figure 1a) to mimic the main masticatory muscles. The determined angles between each linkage and the vertical plane are shown in Figure  1d. Taking advantage of the bi-directional driving of the engineering solution, six PUS linkage realizes both mouth-closing and mouth-opening functions. To avoid singularity (dotted lines in Figure 2d) when the lateral pterygoid PUS chain is near perpendicular to As there are more muscles driving the lower jaw than required, and two TMJs constraining the motion of the condyle, the human masticatory system itself is mechanically actuation redundant [45]. There is a very short distance in the direction which the condyle moves perpendicular to the articular surface. The articular surface of the TMJ is assumed to be rigid and the condyle keeps contact with the surface all the time. The DOFs of the lower jaw is four.

The Redundantly Actuated Chewing Robot
Inspired by human masticatory system, a 6PUS-2HKP RAPR has been proposed [46]. The major challenge in designing the chewing robot deals with the TMJs and the masticatory muscles. In human TMJ, the two condyles are constrained in moving along the fossa surfaces, as shown in Figure 2a. Although a ball-socket joint has three rotational DOFs, it does not allow for translations. While, the condyle actually moves almost freely along the articular surface of the temporal bone in its articular capsule. Therefore, two point-contact HKPs are designed to model the articular contact of human TMJs. The HKP consists of a condylar rod and a fossa structure, as shown in Figure 2b,c. The condylar rods keep contact with the fossa articular surfaces.
Six PUS chains are arranged in parallel corresponding to the resultant line-of-action of six muscles (MAS, TEM and LPT in Figure 1a) to mimic the main masticatory muscles. The determined angles between each linkage and the vertical plane are shown in Figure 1d. Taking advantage of the bi-directional driving of the engineering solution, six PUS linkage realizes both mouth-closing and mouth-opening functions. To avoid singularity (dotted lines in Figure 2d) when the lateral pterygoid PUS chain is near perpendicular to the sliding rail of the prismatic joint, the tilt angle of lateral pterygoid chain is adjusted to slightly larger than its angle of line-of-action. Figure 2e illustrates the kinematical diagram of the 6PUS-2HKP mechanism with one PUS linkage, as the rest of them are the same. Figure 2f shows the composition of one PUS linkage.
The 3D model of the RAPR and its corresponding prototype are shown in Figure 3. This RAPR is bio-inspired by human evolved masticatory system. In particular, the redundant actuation is achieved by adding passive kinematic constraints, the HKPs. The DOFs of the robot reduce while the number of actuators remain the same. The passive point-contact HKP in the RAPR brings new challenges in stiffness analysis and driving force optimization, but it also helps to reproduce the chewing force more realistically and provide an in vitro experimental platform to study the function of TMJs in human masticatory system.    The left and right articular surface Ls and Rs (indicated in Figure 4) are the surfaces where the centers of the left and right condyles can move along. Each curved surface is made up of two sections, which can be described by two quadratic polynomial curves in XM-ZM plane as follows with respect to frame ΣL0 or ΣMa (see Section 2.3 Figure 4). The 3D model of the RAPR and its corresponding prototype are shown in Figure 3. This RAPR is bio-inspired by human evolved masticatory system. In particular, the redundant actuation is achieved by adding passive kinematic constraints, the HKPs. The DOFs of the robot reduce while the number of actuators remain the same. The passive point-contact HKP in the RAPR brings new challenges in stiffness analysis and driving force optimization, but it also helps to reproduce the chewing force more realistically and provide an in vitro experimental platform to study the function of TMJs in human masticatory system.
The left and right articular surface L s and R s (indicated in Figure 4) are the surfaces where the centers of the left and right condyles can move along. Each curved surface is made up of two sections, which can be described by two quadratic polynomial curves in X M -Z M plane as follows with respect to frame Σ L0 or Σ Ma (see Section 2.3 Figure 4).
where x could be x L or x R and z could be z L or z R . They are coordinates of the center of two TMJs. The subscripts L and R represent the left and right TMJ center, respectively. The left and right articular surface Ls and Rs (indicated in Figure 4) are the surfaces where the centers of the left and right condyles can move along. Each curved surface is made up of two sections, which can be described by two quadratic polynomial curves in XM-ZM plane as follows with respect to frame ΣL0 or ΣMa (see Section 2.3 Figure 4).

Reference Frames
Three reference frames are set up in order to describe the motion of the robot, as shown in Figure 4 where x could be xL or xR and z could be zL or zR. They are coordinates of th TMJs. The subscripts L and R represent the left and right TMJ center, respe

Reference Frames
Three reference frames are set up in order to describe the motion o shown in Figure Plane P 2 The pose of end-effector with respect to the base can be described by non-collinear points. The three points may be selected anywhere on th Generally, the platform joints of a PR are in a same plane. While, the redund chewing robot is a 3D spatial one and all spherical joints Mi and universal robot are not in the same plane. The expression of coordinates of platform The pose of end-effector with respect to the base can be described by three arbitrary non-collinear points. The three points may be selected anywhere on the end-effector. Generally, the platform joints of a PR are in a same plane. While, the redundantly actuated chewing robot is a 3D spatial one and all spherical joints M i and universal joints C i of this robot are not in the same plane. The expression of coordinates of platform joints in global frame in terms of three end-effector points will be more complicated than the robot whose platform joints are in same plane.

Inverse Kinematics
The three end-effector points are selected as L, R and I, as shown in Figure 4. The coordinates of platform spherical joints in mandible frame in terms of three end-effector points can be expressed as: where LD = LR × LI. Vector LD is perpendicular to both LR and LI. Then, for each i, the three unknowns (k 1,i , k 2,i and k 3,i ) in Equation (2) can be obtained in terms of initial position coordinates of points L, R, I and M i . In consideration of the HKPs constrains, four generalized variables are able to describe the position and pose of the robot. Assuming G is the generalized variables, four independent parameters x I , y I , z I and x L are chosen to describe the motion of the robot. The rest of three end-effector points' coordinates y L , z L , x R, y R and z R can be calculated based on the geometric constraints. These nine coordinate parameters of three end-effector points are defined as E.
Left HKP keeps contact with constrained surface, we have: where, f indicates the constraint Equation (1). z L can be calculated. The length of LI is constant, then we have: From Equation (4), y L can be obtained. Similarly, as RL and RI all have fixed length and right HKP keeps contact with right constrained surface, the following equations can be written: where x R, y R and z R can also be obtained by using Equation (5). Then the coordinates of platform joints M i can be calculated from Equation (2). Based on the geometric constraints of parallel kinematic chains, the joint variables Q i (i=1, 2, . . . , 6) of the 6 PUS kinematic chains can be obtained [46].

Velocity Analysis and Jacobian Matrix
Jacobian matrix is the generalized transmission ratio of the velocity or force transmission from the input of the joint variables to the output of the generalized variables [47]. This RAPR has mixed translational and rotational DOFs. For this type of RAPR, this section presents a dimensionally homogeneous non-square Jacobian matrix formulation method based on three end-effector points.
For the robot studied in this paper, the length of leg C i M i is constant and leads to: Differentiating (6) yields: Machines 2021, 9, 171 8 of 22 where V M i is the velocity of spherical joint, and V C i the velocity of prismatic joint, indicated by . Q i . The position of M i in global frame can be expressed as: Differentiating (8) gives: where x I , Combining (7) and (9), yields: where z L ] · LI (11) in which, LI = [x LI , y LI , z LI ] T . The instantaneous velocity vector of a condyle point is in a plane that goes through the condyle center and is tangent to the curved surface as shown in Figure 4, Plane P 1 and P 2 .
The slope of P 1 in X m -Z m plane can be calculated by the equation of left curved surface L s . So does the slope of P 2 : Since the resultant velocity vector V L must be in the plane P 1 , the slope m P 1 can also be calculated by: z R can also be obtained. Writing in matrix form, we have: where P is a matrix that can be obtained by velocity projection theory aforementioned, and Equation (10) can be rewritten as: Combining (14) and (15), gives: .
The inverse Jacobian matrix of the RAPR can be written as:

Stiffness Analysis
Stiffness is a very important performance index in parallel robots [48]. Redundant actuation has been proved to be an effective way to improve its stiffness [26]. The redundant actuation of this 6PUS-2HKP RAPR is achieved by adding two passive point-contact constraints HKPs to the 6-PUS parallel mechanism. In order to investigate the effect of pointcontact HKP constraint on the redundantly actuated chewing robot, the stiffness analysis of redundantly actuated 6PUS-2HKP and non-redundantly actuated 6-PUS mechanism are compared.
The relation between Cartesian generalized velocity and joint velocity can be obtained by combining Equations (16) and (17): The infinitesimal displacements in the Cartesian space and joint space can be written as: Based on virtual work principle, the relation between external force F and prismatic joint force τ can be written as: The joint force τ can also be computed by: where K C is a diagonal matrix which consists of stiffness of the driving system. Combining Equations (19)-(21), the following equation can be obtained.
Therefore, the stiffness matrix of the RAPR can be defined as: The generalized infinitesimal displacements of the end-effector can be expressed by: In this derivation, the basic assumptions in this paper are: (1) The end-effector and the base are assumed as rigid bodies, as the rigidities of them are much larger than that of the six driving rods; (2) the weights of the legs are neglected; (3) the joints are assumed rigid and frictionless; (4) the stiffness of the motors and coupler is equal and is assumed as rigid. All of the six motors actuate the six prismatic joints through ball screw and they have the same rigidity. In one PUS linkage, only the ball screw and the driving legs are taken into account. As the prismatic joints are the direct actuated joints, the model can be simplified as follows. The stiffness matrix of one linkage is approximately calculated by the following equation: where K C is the scalar stiffness parameter. K bs and K dr stand for the stiffness parameters of the ball screw and driving rod, respectively. K bs and K dr can be computed as: where E is the Young's modulus of the material, A the sectional area of the link, L the length of the link, and j = bs, dr. The scalar stiffness parameters of six linkages can computed as follows. K C1 = K C2 = 2.11 × 10 7 N/m, K C3 = K C4 = 1.88 × 10 7 N/m, K C5 = K C6 = 2.57 × 10 7 N/m. In the numerical simulation, a mouth-opening trajectory was implemented to estimate the deformation displacement occurring to 6PUS-2HKP RAPR and 6-PUS PR. The motion of the RAPR in terms of the incisor point (point I) with respect to frame Σ B is shown in Figure 5. The motion lasted 1 s and a resultant force (80 N in Z B -direction and 20 N in X B -direction) was continuously applied to point I. During this process, the deformation displacements ∆G at different configuration of the robot were calculated, as shown in Figure 6. The largest deformation displacement is 4.047 × 10 −6 m in Z I direction and the smallest one is 2.573 × 10 −7 m in -Y I direction. It can be seen that the RAPR at each configuration has small infinitesimal displacements, which means that it owns a high stiffness. In the numerical simulation, a mouth-opening trajectory was implemented to estimate the deformation displacement occurring to 6PUS-2HKP RAPR and 6-PUS PR. The motion of the RAPR in terms of the incisor point (point I) with respect to frame ΣB is shown in Figure 5. The motion lasted 1 s and a resultant force (80 N in ZB-direction and 20 N in XB-direction) was continuously applied to point I. During this process, the deformation displacements G at different configuration of the robot were calculated, as shown in Figure 6. The largest deformation displacement is 4.047 × 10 −6 m in ZI direction and the smallest one is 2.573 × 10 −7 m in -YI direction. It can be seen that the RAPR at each configuration has small infinitesimal displacements, which means that it owns a high stiffness.
To highlight the role of redundant actuation in improving the stiffness, the deformation displacement of the non-redundantly actuated 6-PUS PR was also computed. The displacement of the 6-PUS PR under the same mouth-opening trajectory is shown in Figure 7. Having the same geometric parameters with the 6PUS-2HKP RAPR, the smallest deformation displacement of the non-redundantly actuated 6-PUS PR is 1.451 × 10 −5 in XI direction, which is larger than any deformation displacements of the RAPR at each configuration. By comparison of these two robots, one can see that the redundantly actuated 6PUS-2HKP RAPR experiences smaller deformation displacements under the same force. It can be concluded from the numerical simulation that the 6PUS-2HKP RAPR with two HKPs mimicking human TMJs owns better stiffness than the non-redundantly actuated 6-PUS PR counterpart.   In the numerical simulation, a mouth-opening trajectory was implemented to estimate the deformation displacement occurring to 6PUS-2HKP RAPR and 6-PUS PR. The motion of the RAPR in terms of the incisor point (point I) with respect to frame ΣB is shown in Figure 5. The motion lasted 1 s and a resultant force (80 N in ZB-direction and 20 N in XB-direction) was continuously applied to point I. During this process, the deformation displacements G at different configuration of the robot were calculated, as shown in Figure 6. The largest deformation displacement is 4.047 × 10 −6 m in ZI direction and the smallest one is 2.573 × 10 −7 m in -YI direction. It can be seen that the RAPR at each configuration has small infinitesimal displacements, which means that it owns a high stiffness.
To highlight the role of redundant actuation in improving the stiffness, the deformation displacement of the non-redundantly actuated 6-PUS PR was also computed. The displacement of the 6-PUS PR under the same mouth-opening trajectory is shown in Figure 7. Having the same geometric parameters with the 6PUS-2HKP RAPR, the smallest deformation displacement of the non-redundantly actuated 6-PUS PR is 1.451 × 10 −5 in XI direction, which is larger than any deformation displacements of the RAPR at each configuration. By comparison of these two robots, one can see that the redundantly actuated 6PUS-2HKP RAPR experiences smaller deformation displacements under the same force. It can be concluded from the numerical simulation that the 6PUS-2HKP RAPR with two HKPs mimicking human TMJs owns better stiffness than the non-redundantly actuated 6-PUS PR counterpart.   To highlight the role of redundant actuation in improving the stiffness, the deformation displacement of the non-redundantly actuated 6-PUS PR was also computed. The displacement of the 6-PUS PR under the same mouth-opening trajectory is shown in Figure 7. Having the same geometric parameters with the 6PUS-2HKP RAPR, the smallest deformation displacement of the non-redundantly actuated 6-PUS PR is 1.451 × 10 −5 in X I direction, which is larger than any deformation displacements of the RAPR at each configuration. By comparison of these two robots, one can see that the redundantly actuated 6PUS-2HKP RAPR experiences smaller deformation displacements under the same force. It can be concluded from the numerical simulation that the 6PUS-2HKP RAPR with two HKPs mimicking human TMJs owns better stiffness than the non-redundantly actuated 6-PUS PR counterpart.

Optimization of Driving Force
The inclusion of HKP into 6-PUS mechanism brings new features to the c robot. On the one hand, redundant actuation of the RAPR causes the inverse dyna robot lacks a unique solution. On the other hand, it also allows redistribution of force by dynamics. This section will investigate the effect of point-contact HKP con in this redundant actuation system from the perspective of dynamics. The driving of the RAPR are optimized by finding the minimum Euclidean norm solutio optimization flow chart is shown in Figure 8. The inverse kinematics has been dis in Section 3.1. From this foundation, the velocity, kinetic energy, potential energy six PUS linkages and end-effector can be calculated, the details of which are not pr due to the limited space.

Constraints of two HKPs
Inverse kinematics

Optimization of Driving Force
The inclusion of HKP into 6-PUS mechanism brings new features to the chewing robot. On the one hand, redundant actuation of the RAPR causes the inverse dynamics of robot lacks a unique solution. On the other hand, it also allows redistribution of driving force by dynamics. This section will investigate the effect of point-contact HKP constraint in this redundant actuation system from the perspective of dynamics. The driving forces of the RAPR are optimized by finding the minimum Euclidean norm solution. The optimization flow chart is shown in Figure 8. The inverse kinematics has been discussed in Section 3.1. From this foundation, the velocity, kinetic energy, potential energy of the six PUS linkages and end-effector can be calculated, the details of which are not provided due to the limited space.
Lagrange's equations are used to derive the dynamics of the RAPR. The Lagrange's equations of the first type can be written in two sets. The first equation set can be written as: where L is the Lagrange function; Q gj represents the generalized force on the end-effector and k is the number of generalized coordinates, k = 4; i indicates the number of constraint equations, i = 6; n is the sum of i and k, n = 10; Γ i is the ith constraint equation; q j indicates the four generalized position parameters X L , Y L , β L and γ L in Cartesian space. βL and γ L are the Euler angles. Let G C be [x L , y L , β L , γ L ] T . As the robot has four DOFs, the generalized forces Q gj = [F gx , F gy , M gβ , M gγ ] T are defined as forces in X L and Y L direction and torques around Y L and Z L . Considering the geometric constraints of each kinematic chain, link C i M i has a fixed length L i and leads to: The second equation set, relating the driving torque/force (driving force of slider in this paper), can be expressed by: where τ i (Q aj ) is the driving force. According to the optimization algorithm of pseudo-inverse method, the minimum 2-norm solution τ 0 can be obtained: . τ 0 denotes the six minimum 2-norm driving forces.
The trajectory of left lateral movement was analyzed in this simulation, shown in Figure 9. The position, velocity, and acceleration of the six prismatic joints for the given pose of mandible can be determined by inverse kinematics.

Optimization of Driving Force
The inclusion of HKP into 6-PUS mechanism brings new features to the c robot. On the one hand, redundant actuation of the RAPR causes the inverse dyna robot lacks a unique solution. On the other hand, it also allows redistribution of force by dynamics. This section will investigate the effect of point-contact HKP co in this redundant actuation system from the perspective of dynamics. The drivin of the RAPR are optimized by finding the minimum Euclidean norm soluti optimization flow chart is shown in Figure 8. The inverse kinematics has been di in Section 3.1. From this foundation, the velocity, kinetic energy, potential energ six PUS linkages and end-effector can be calculated, the details of which are not p due to the limited space.

Constraints of two HKPs
Inverse kinematics Lagrange's equations are used to derive the dynamics of the RAPR. The Lag equations of the first type can be written in two sets. The first equation set can be as: is the right pseudo-inverse matrix of . τ 0 denotes the six minimum 2-norm driving forces. The trajectory of left lateral movement was analyzed in this simulation, shown in Figure 9. The position, velocity, and acceleration of the six prismatic joints for the given pose of mandible can be determined by inverse kinematics. Following the same lateral movement, the driving force and driving power of 6PUS-2HKP RAPR and its counterpart non-redundantly actuated 6-PUS parallel robot were solved by dynamics equations for comparison. The driving force and power of 6-PUS robot without optimization algorithm are illustrated in Figure 10 and Table 1. The driving forces from the most significant to the least significant in sequence are masseter, lateral pterygoid, and temporalis. The maximum force is that of LM which is 10.78 N and the minimum force is that of RT which is 6.27 N. As the velocity of RLP is the largest, the RLP has a maximum power of 0.487 W.
The driving forces of the 6PUS-2HKP RAPR adopting pseudo-inverse optimization method are shown in Figure 11 and Table 1. Compared with the 6-PUS robot, the driving forces of the six prismatic joints reduces dramatically, especially the lateral pterygoid and the masseter. The maximum force is that of RM which is 10.26 N and the minimum force occurs at the RLP which has a value of 3.36 N. As we can see from the standard deviation of the power in Table 1, the driving power also becomes more balanced and the power consumption has improved. The RT has a maximum power 0.224 W. None of any kinematic chain of 6PUS-2HKP RAPR has a very sharp power curve compared with the 6-PUS robot.  Following the same lateral movement, the driving force and driving power of 6PUS-2HKP RAPR and its counterpart non-redundantly actuated 6-PUS parallel robot were solved by dynamics equations for comparison. The driving force and power of 6-PUS robot without optimization algorithm are illustrated in Figure 10 and Table 1. The driving forces from the most significant to the least significant in sequence are masseter, lateral pterygoid, and temporalis. The maximum force is that of LM which is 10.78 N and the minimum force is that of RT which is 6.27 N. As the velocity of RLP is the largest, the RLP has a maximum power of 0.487 W. Following the same lateral movement, the driving force and driving power of 6PUS-2HKP RAPR and its counterpart non-redundantly actuated 6-PUS parallel robot were solved by dynamics equations for comparison. The driving force and power of 6-PUS robot without optimization algorithm are illustrated in Figure 10 and Table 1. The driving forces from the most significant to the least significant in sequence are masseter, lateral pterygoid, and temporalis. The maximum force is that of LM which is 10.78 N and the minimum force is that of RT which is 6.27 N. As the velocity of RLP is the largest, the RLP has a maximum power of 0.487 W.
The driving forces of the 6PUS-2HKP RAPR adopting pseudo-inverse optimization method are shown in Figure 11 and Table 1. Compared with the 6-PUS robot, the driving forces of the six prismatic joints reduces dramatically, especially the lateral pterygoid and the masseter. The maximum force is that of RM which is 10.26 N and the minimum force occurs at the RLP which has a value of 3.36 N. As we can see from the standard deviation of the power in Table 1, the driving power also becomes more balanced and the power consumption has improved. The RT has a maximum power 0.224 W. None of any kinematic chain of 6PUS-2HKP RAPR has a very sharp power curve compared with the 6-PUS robot.   Table 1. Comparison of driving force and power between 6-PUS and 6PUS-2HKP mechanism (F min , F max , F mean -The minimum, maximum and mean value of the driving force, respectively; P max -The maximum value of the power; P SD -The standard deviation of the power). The driving forces of the 6PUS-2HKP RAPR adopting pseudo-inverse optimization method are shown in Figure 11 and Table 1. Compared with the 6-PUS robot, the driving forces of the six prismatic joints reduces dramatically, especially the lateral pterygoid and the masseter. The maximum force is that of RM which is 10.26 N and the minimum force occurs at the RLP which has a value of 3.36 N. As we can see from the standard deviation of the power in Table 1, the driving power also becomes more balanced and the power consumption has improved. The RT has a maximum power 0.224 W. None of any kinematic chain of 6PUS-2HKP RAPR has a very sharp power curve compared with the 6-PUS robot.    The highlight of the role of the point-contact HKP constraints in the 6PUS-2HKP RAPR is homogenizing driving force and driving power. The driving force decreases and the driving power becomes more homogenous by applying the pseudo-inverse optimization method.

Measurement Methods for Chewing Simulated Food
In this paper, the reaction forces on the occlusal surface of incisor and on the mandibular condyle are regarded as the occlusal forces and the TMJ (HKP) forces, respectively. For the purpose of studying the relation between occlusal force and TMJ force of the chewing robot, a real-time force measurement system was built. Due to the size of the point-contact HKP constraints, the strain gauge is selected as the force sensor. Two strain gauges are glued perpendicular to each other on the upper jaw, as shown in Figure 12a. The vertical one is used for measuring the occlusal force and the lateral one is for lateral force on teeth. The other two strain gauges are glued on left condylar rod and right condylar rod, respectively, as shown in Figure 12b. These two strain gauges sense the HKP force in real time. As the chewing robot was intended to perform a variety of trajectories betw lateral and vertical chewing, the trajectories may vary. A desired chewing trajectory profiled according to characteristics of human chewing trajectory. The given trajecto of the mandible are determined by the four independent position parameters. The mo controller computes the kinematics in a time interval of 2 ms to get the parameters of six servo motors. The chewing trajectory of the chewing robot is shown in Figure 13 w respect to frame ΣB.
Simulated food is adopted to conduct chewing test on the chewing robot in orde obtain the occlusal force and TMJ force. Two kinds of materials silicone (SI) polyethylene (PE) are chosen as the simulated food, as shown in Figure 14. As diffe thicknesses of simulated food can produce different value of force, silicone with 4 mm mm and 12 mm thick and polyethylene with 1 mm, 2 mm thick are adopted during chewing cycle in this paper.  As the chewing robot was intended to perform a variety of trajectories between lateral and vertical chewing, the trajectories may vary. A desired chewing trajectory was profiled according to characteristics of human chewing trajectory. The given trajectories of the mandible are determined by the four independent position parameters. The motion controller computes the kinematics in a time interval of 2 ms to get the parameters of the six servo motors. The chewing trajectory of the chewing robot is shown in Figure 13 with respect to frame Σ B .  As the chewing robot was intended to perform a variety of trajectories between lateral and vertical chewing, the trajectories may vary. A desired chewing trajectory was profiled according to characteristics of human chewing trajectory. The given trajectories of the mandible are determined by the four independent position parameters. The motion controller computes the kinematics in a time interval of 2 ms to get the parameters of the six servo motors. The chewing trajectory of the chewing robot is shown in Figure 13 with respect to frame ΣB.
Simulated food is adopted to conduct chewing test on the chewing robot in order to obtain the occlusal force and TMJ force. Two kinds of materials silicone (SI) and polyethylene (PE) are chosen as the simulated food, as shown in Figure 14. As different thicknesses of simulated food can produce different value of force, silicone with 4 mm, 8 mm and 12 mm thick and polyethylene with 1 mm, 2 mm thick are adopted during the chewing cycle in this paper.  Simulated food is adopted to conduct chewing test on the chewing robot in order to obtain the occlusal force and TMJ force. Two kinds of materials silicone (SI) and polyethylene (PE) are chosen as the simulated food, as shown in Figure 14. As different thicknesses of simulated food can produce different value of force, silicone with 4 mm, 8 mm and 12 mm thick and polyethylene with 1 mm, 2 mm thick are adopted during the chewing cycle in this paper. In the experiment, the forces in one chewing cycle was rec prototype. At the opening phase, the simulated food was put the lower jaw (incisor position). When the mouth closed simulated food was chewed by the robot. Figure 15

Experimental Results of Occlusal Force and TMJ Force
Once the chewing robot started chewing the simulated structure and the condyle rod were stressed. The occlusal for by strains were collected by the strain data acquisition system In the experiment, the forces in one chewing cycle was recorded on the chewing robot prototype. At the opening phase, the simulated food was put between the upper jaw and the lower jaw (incisor position). When the mouth closed to the occlusal phase, the simulated food was chewed by the robot. Figure 15 shows a silicone simulated food is being chewed by the chewing robot. Different simulated food was tested one by one. In the experiment, the forces in one chewing cycle was recorded prototype. At the opening phase, the simulated food was put betw the lower jaw (incisor position). When the mouth closed to th simulated food was chewed by the robot. Figure 15 shows a silic being chewed by the chewing robot. Different simulated food was

Experimental Results of Occlusal Force and TMJ Force
Once the chewing robot started chewing the simulated food structure and the condyle rod were stressed. The occlusal force an by strains were collected by the strain data acquisition system durin The strains of the four sensors with 1 mm thick polyethylene in o shown in Figures 16 and 17.

Experimental Results of Occlusal Force and TMJ Force
Once the chewing robot started chewing the simulated food, both the upper jaw structure and the condyle rod were stressed. The occlusal force and TMJ force expressed by strains were collected by the strain data acquisition system during the chewing course. The strains of the four sensors with 1 mm thick polyethylene in one chewing cycle are shown in Figures 16 and 17.

Experimental Results of Occlusal Force and TMJ Force
Once the chewing robot started chewing the simulated food, both the upper jaw structure and the condyle rod were stressed. The occlusal force and TMJ force expressed by strains were collected by the strain data acquisition system during the chewing course. The strains of the four sensors with 1 mm thick polyethylene in one chewing cycle are shown in Figures 16 and 17.  The force can be calculated with the following methods. The material of the TMJs and upper jaw is aluminum alloy and the modulus of elasticity of them is 0.7 × 10 5 MPa. Based on the definitions of stress and strain and Hooke's Law, the stress can be obtained: 11 6 5 E 0.7 10 12.5 10 8.75 10 where σ is the stress of the upper jaw. E is the Young's modulus of aluminum alloy and ε is the strain shown in Figure 16b. The occlusal force can be calculated as follows: 5 4 8.75 10 1.2 10 105 where F is the occlusal force and A is the cross-sectional area of the middle position of strain gauge on upper jaw (Figure 12a). Similarly, the cross-sectional area of the condylar block is 49 × 10 −6 m 2 and the TMJ force can also be calculated. Multi-group experiments were carried out and the occlusal forces and TMJ forces were also calculated. The results are shown in Figures 18 and 19. As can be seen from the bar chart, the occlusal force increases with thicknesses of two kinds of simulated food increases directly. The increase of the occlusal force causes the TMJ forces become larger. The results verify that the loading of the TMJ increases along with occlusal force.
The maximum occlusal force occurred when the robot chews the 2-mm thick polyethylene, which is 303.2 N and the corresponding left and right TMJ force is 70.0 N and 72.7 N, respectively. The minimum occlusal force is 104.2 N and corresponding left The force can be calculated with the following methods. The material of the TMJs and upper jaw is aluminum alloy and the modulus of elasticity of them is 0.7 × 10 5 MPa. Based on the definitions of stress and strain and Hooke's Law, the stress can be obtained: where σ is the stress of the upper jaw. E is the Young's modulus of aluminum alloy and ε is the strain shown in Figure 16b. The occlusal force can be calculated as follows: where F is the occlusal force and A is the cross-sectional area of the middle position of strain gauge on upper jaw (Figure 12a). Similarly, the cross-sectional area of the condylar block is 49 × 10 −6 m 2 and the TMJ force can also be calculated. Multi-group experiments were carried out and the occlusal forces and TMJ forces were also calculated.
The results are shown in Figures 18 and 19. As can be seen from the bar chart, the occlusal force increases with thicknesses of two kinds of simulated food increases directly. The increase of the occlusal force causes the TMJ forces become larger. The results verify that the loading of the TMJ increases along with occlusal force.  The maximum bite force that human can apply to the molars which has been measured about 700 N or higher in healthy people [49,50], approximately equal to human's own weight. While the size of the condyle is only about 20 × 10 × 10 mm [33], which is too small to burden such a large loading compared with the knee joint of human body. From the perspective of theory of mechanism, redundant actuation allows to partially control the internal forces and the actuating force and power consumption of actuators can be redistributed by optimization methods.
Redundant actuation in masticatory system may help us to reveal the physiological phenomenon. In this experiment, it shows the increment of occlusal force has priority over the TMJ force. It means the two point-contact HKPs in the redundant actuated system have a positive effect on generating process of chewing force caused by actuators. This suggests the forces of these two TMJs play a role in optimizing the internal force of the redundant actuated system, enabling the muscle group to output more efficient occlusal

Force (N)
Occlusal force Right HKP force  The maximum bite force that human can apply to the molars which has been measured about 700 N or higher in healthy people [49,50], approximately equal to human's own weight. While the size of the condyle is only about 20 × 10 × 10 mm [33], which is too small to burden such a large loading compared with the knee joint of human body. From the perspective of theory of mechanism, redundant actuation allows to partially control the internal forces and the actuating force and power consumption of actuators can be redistributed by optimization methods.
Redundant actuation in masticatory system may help us to reveal the physiological phenomenon. In this experiment, it shows the increment of occlusal force has priority over the TMJ force. It means the two point-contact HKPs in the redundant actuated system have a positive effect on generating process of chewing force caused by actuators. This suggests the forces of these two TMJs play a role in optimizing the internal force of the redundant actuated system, enabling the muscle group to output more efficient occlusal force during chewing process.

Force (N)
Occlusal force Right HKP force The maximum occlusal force occurred when the robot chews the 2-mm thick polyethylene, which is 303.2 N and the corresponding left and right TMJ force is 70.0 N and 72.7 N, respectively. The minimum occlusal force is 104.2 N and corresponding left and right TMJ force is 32.9 N and 38.4 N when the robot bites the 4-mm thick silicone. The occlusal force increases by 190.98% from silicone with 4 mm thick to polyethylene with 2 mm thick. Nevertheless, the corresponding left TMJ force only increases by 120.97% and the right TMJ force only increases by 82.29%. It can be seen from the changing process of the two forces that the growth rate of occlusal force is larger than that of TMJ force.
The maximum bite force that human can apply to the molars which has been measured about 700 N or higher in healthy people [49,50], approximately equal to human's own weight. While the size of the condyle is only about 20 × 10 × 10 mm [33], which is too small to burden such a large loading compared with the knee joint of human body. From the perspective of theory of mechanism, redundant actuation allows to partially control the internal forces and the actuating force and power consumption of actuators can be redistributed by optimization methods.
Redundant actuation in masticatory system may help us to reveal the physiological phenomenon. In this experiment, it shows the increment of occlusal force has priority over the TMJ force. It means the two point-contact HKPs in the redundant actuated system have a positive effect on generating process of chewing force caused by actuators. This suggests the forces of these two TMJs play a role in optimizing the internal force of the redundant actuated system, enabling the muscle group to output more efficient occlusal force during chewing process.
Besides, the average value of occlusal force on incisor of a healthy young adult is about 146.17 N [51]. As we can see the force that the chewing robot can produce in this experiment is much larger than the average force. The robot is able to output enough occlusal force which can be potentially used in new dental materials test or food properties test. The internal TMJ forces during chewing provides the possibility to study the performance test of TMJ prosthesis.

Conclusions
The physical redundantly actuated chewing robot was specified in terms of six PUS linkage and two HKPs. The introduction of HKP allows the chewing robot to be able to reproduce the TMJ movement and influence of the TMJ force on occlusal force. Accordingly, it can improve the bio-imitability of the chewing robot. Unlike the traditional ways of achieving redundant actuation, this RAPR was realized through the inclusion of additional constraints to reduce DOFs, resulting in the number of actuators is larger than the number of DOFs.
In order to investigate the advantage of musculoskeletal actuation redundancy in human masticatory system, such as optimal stiffness, energy efficiency and chewing forces needed for the food breakdown, this paper studies the stiffness and optimization of driving force and presents the chewing experiment of the RAPR. In stiffness analysis, the three endeffector points method was proposed for use in differential kinematics of the RAPR with mixed DOFs, mapping the reduced number of the end-effector's generalized variables to the joint variables. Based on the dimensionally homogeneous Jacobian matrix, the stiffness simulation results confirmed the 6PUS-2HKP RAPR with two point-contact HKPs has better stiffness than the non-redundantly actuated counterpart 6-PUS PR. In optimization of driving force, Lagrange's equations were used to derive the dynamics of the RAPR. The driving force of the RAPR was optimized by finding the minimum Euclidean norm solution. By comparison, one can get the same result that the 6PUS-2HKP RAPR with two HKPs performs better than the 6-PUS PR. These simulation results of chewing robot show the musculoskeletal structure of human masticatory system itself has the advantage of good stiffness and energy saving. They also reveal the significance of TMJs in the redundantly actuated human masticatory system.
As the TMJs are modeled as two point-contact HKPs in the proposed robot, the occlusal force generated has more realistic implication. The experiments of chewing simulated food with different thicknesses are carried out to evaluate the occlusal force and the TMJ force. The occlusal force increases with the thicknesses of simulated food, and the corresponding TMJ force also increases, but at a slower rate. From a biomechanical point of view, the experiment also illustrates that the existence of TMJ in human redundantly actuated masticatory system ensures more efficient transmission and transformation from muscle force to occlusal force.
Ongoing research on the robot prototype involves a wear-resistant test of the alloy ceramic crown, which evaluates the wear condition of the occlusal surfaces of the crowns. Future work may focus on developing a new generation of chewing robot. As in the current version, rigid linkages and rigid contacts are used to mimic the muscles and TMJs, respectively, the chewing force and TMJ force generated by the chewing robot are still not completely coincident with human masticatory system. Emphasis in future studies will involve adopting parallel flexible drive system considering the muscle elasticity, and designing soft TMJ structure considering the disc and ligament in the chewing robot.