Kinematic and Stiffness Modeling of a Novel 3-DOF RPU+UPU+SPU Parallel Manipulator

A novel 3-degrees of freedom (DOF) RPU+UPU+SPU parallel manipulator (PM) is proposed in this study, and its complete kinematics and stiffness are studied systematically. First, the architecture description is discussed, and the inverse and the forward positional posture analysis are studied based on the constraints in the PM. Second, the Jacobian matrix, the velocity model, the Hessian matrix and the acceleration model are derived in explicit and compact forms. Third, based on the virtual work principle, the static model and the deformation decompose method, the stiffness model is built. Meanwhile, the stiffness matrix and the compliance matrix are obtained. Finally, the correctness of the models built in this study are verified by simulative PMs. This study is expected to provide new ideas for the design of PM machine tools.


I. INTRODUCTION
As an essential branch of limited-DOF PMs, 3-DOF PMs have attracted much attention due to their excellent merits, such as simplicity in structure, low cost of manufacturing and easy control [1]. According to these characteristics, 3-DOF PMs are widely used in the field of machine tool. In this field, one of the famous inventions is a 3PRS PM [2], which is the primary mechanism of Sprint Z3. In addition, Neumann proposed a 3UPS+UP PM (Tricept) [3], which owns high dynamics, stiffness and ample workspace. Then, the inventor presented the concept of a 2UPR+SPR PM (Exechon) [4], which can be applied to fill in the gap between the traditional machine tools and the industrial robots.
In recent years, the design and analysis of 3-DOF PMs are still a hot topic in academia. By integrating one active UPU leg into a passive one, Huang et al. [5] designed a modified Tricept robot-TriVariant. Li et al. [6] presented an A3 head applied to manufacturing large structural components. Wang et al. [7] proposed a 3PUU PM possessing good motion/force transmissibility and orientation capability, and studied its optimal design [8]. Sun synthesized some 3-DOF PMs [9], and presented a simple and highly visual approach The associate editor coordinating the review of this manuscript and approving it for publication was Guilin Yang . for type synthesis of 3-DOF over-constraint PMs [10]. Lu et al. [11] innovated some 3-DOF PMs with planar sub-chains using the revised digital topological graphs and arrays. Gallardo and Rodriguez [12] proposed a 3-DOF 3RPRRC+RRPRU robot. Jing et al. [13] designed a redundant collaborative manipulator containing a class of 3-DOF PM heads with high rotational capability. Hu et al. [14] studied the kinematic characteristics of a 3-DOF 3UPU+UP PM with coupling parallel platforms. Yang et al. [15] discussed the kinematics of a 3-DOF 3PPS PM. Here and throughout, R, P, U and S denotes a revolute joint, a prismatic joint, a universal joint and a sphere joint, respectively. However, existing 3-DOF PMs are not enough to satisfy the continuous innovation demands of the machine tools. Especially in recent years, the demands for irregular shape components in the aerospace field are increased, the stiffness requirements of each direction are inconsistent during processing. Therefore, the demands for innovative design of asymmetric 3-DOF PMs are increased significantly. In addition, existing asymmetric 3-DOF PMs virtually all contain over-constraint wrenches, which have high requirements for manufacture and assembly. Once the manufacture and assembly are slightly careless, it will seriously affect the accuracy of the machine tools. Compared with over-constraint PMs, non-over-constraint PMs have a simple structure and lower requirements in manufacture and assembly, which can better adapt to different working requirements and realize specific requirements. However, so far there are few reports for asymmetric non-over-constraint 3-DOF PMs, which provide motivations for this study.
Stiffness refer to the ability to resist elastic deformations, which is one of the PM's most important performance parameters and a core factor that must be considered in the machine tool design. Many scientists commit to studying the stiffness of PM. Gosselin [16] established the stiffness model of PMs considering the active factors. Huang et al. [17] studied the stiffness of a tripod-based PM containing the rigidity of the machine frame. Liu et al. [18] conducted an optimum design of a 3-DOF spherical PM concerning the conditioning and stiffness indices. Zhang and Gosselin [19] and Zhang [20] built the kinetostatic and stiffness model for some 3, 4 and 5-DOF PMs. Han et al. [21] analyzed the stiffness of a 4-DOF PM basing on the screw theory. Merlet [22], a famous French mechanism scholar, pointed out that: for lack of considering the role of constraints, many mechanism analysis problems may exist faultiness. For example, constrained forces/torques are bound to cause deformations of the leg and then significantly affect the stiffness of PMs. Therefore, they must be considered when stiffness modeling. Unfortunately, the previously mentioned researches have ignored.
In recent years, the influences of constraints on the force and stiffness of limited-DOF PMs is gratifying caused extensive concerns. Wojtyra [23] used the singular value decomposition and QR decomposition methods to solve the joint constraint reaction forces of over-constraint PMs. Utilizing the instantaneous screw theory, Lian et al. [24] formulated the stiffness model of a 5-DOF PM considering the gravitational effects. Sun et al. [25] established the semianalytic stiffness model of a hybrid manipulator as a friction stir welding robot composed of a 3-DOF PM module and a 2-DOF rotating head. Hu and Huang [26] proposed the constraints and deformations decomposition matching method for solving the force and stiffness problems of limited-DOF PMs. Liu et al. [27] established the static and stiffness models of over-constraint PMs by combining the weighted Moore-Penrose inverse. Li et al. [28] studied the analytic solution of the elastic stiffness model for limited-DOF PMs using the geometric algebra and strain energy methods. Cao and Ding [29] solved the joint reaction forces of an over-constraint PM with flexible joints. Considering both constrained wrenches and active wrenches, Li and Xu [30] presented the stiffness characteristics of a 3-DOF 3-PUU translational PM; Chen et al. [31] systematically studied the stiffness of a 6-DOF 3CPS PM; Shan and Hen [32] investigated the stiffness of a 2(3PUS+S) PM with two moving platforms. Each of the above method has its own merits, which lays solid foundations for this study.
For these reasons, this study proposes a novel asymmetric non-over-constraint 3-DOF RPU+UPU+SPU PM, the complete kinematics and stiffness analysis of this novel PM are carried out, and is expected to provide new ideas for the design of PM machine tools. Since this novel PM has three different types of legs, and each leg contains different constraints, the researches of the complete kinematic analysis and stiffness analysis are still challenging works.
The remainder of this study is organized as below. In section II, the architecture description of the proposed PM is given, then the constraints analysis is performed. In section III, the complete positional posture analysis is discussed, then the velocity and acceleration kinematics analysis are studied. In section IV, the stiffness model is established. In section V, numerical examples are given to verify the correctness of the analytic models established in this study. Finally, conclusions are drawn.

II. DESCRIPTION OF THE RPU+UPU+SPU PM A. MECHANISM ARCHITECTURE
The proposed PM includes a moving platform m, a base platform B and three legs r i (i = 1, 2, 3), as shown in FIGURE 1. The architectures of m and B are equilateral triangles. Each r i (i = 1, 2, 3) has one active P joint. Each U joint is composed of two perpendicularly intersecting R joints.
In r 1 , the bottom of the P joint is attached to B by a R joint whose axis is in the plane of B and perpendicular to one side of B. The other end of the P joint is fixed to m with a U joint. In this U joint, one R joint is parallel with the R joint in the plane of B, the remaining one R joint is perpendicular to m.
In r 2 , the bottom of the P joint is attached to B by a U joint. In this U joint, one R joint is perpendicular to B, the remaining one R joint is vertical to r 2 . The other end of the P joint is fixed to m with a U joint. In this U joint, one R joint is parallel with the R joint that is perpendicular to r 2 , the remaining one R joint is in the plane of m and perpendicular to one side of m.
In r 3 , the bottom of the P joint is attached to B by a S joint. The other end of the P joint is fixed to m with a U joint.  The mechanism under study is a 3-DOF PM, as can be calculated by the revised Kutzbach-Grübler equation [1]: where M is the DOF number of the PM, n is the number of links in the PM, g is the number of joints, m i is the DOF number of the i-th joint and m 0 is the passive DOF. For the proposed PM, there are one S joint, four U joints, one R joint and three P joints, n = 8 and g = 9. Application of (1), it leads to: O is the center of m. X is parallel with A 1 A 3 . Y also lies in the plane of m while Z is normal to m and points upward, thereby forming a right-handed orthogonal frame. As shown in FIGURE 2. Let R ij be the j-th R joint from B to m in r i (i = 1, 2). Based on the mechanism architecture, following relationships can be written: where and ⊥ denotes parallel and perpendicular constraints, respectively. Because of the above relationships, there are constrained wrenches in the proposed PM, which can be determined by the geometrical rules [33], [34]: (a) In each leg, the constrained forces should be perpendicular to all P joints and coplanar with all R joints.
(b) In each leg, the constrained torques should be perpendicular to all R joints.
Utilizing the rules (a) and (b), the constrained forces/ torques in this PM can be determined. In r 1 , there exists one constrained force F p1 which is parallel with R 12 and R 11 ,and one constrained torque T 1 which is perpendicular to R 12 and R 13 as well as passes through A 1 . In r 2 , there exists one constrained force F p2 which is parallel with R 23 and R 22 as well as passes through a point C. Where C is an intersection point of R 21 and R 24 . As shown in FIGURE 2.
To simplify expression, denote δ i be the unit vector of r i , e i be the vector from O to A i , F pi /T 1 be the value of F pi / T 1 , f i / τ 1 be the unit vector of F pi /T 1 , d i be the vector from O to an arbitrary point on F pi . As shown in FIGURE 2.

A. COMPLETE POSITIONAL POSTURE MODEL
In limited-DOF PMs, there exists coupling relationships between the 6-dimensional positional posture of the moving platform, which need to be analyzed first before the kinematic modeling. Denote R ij , X, Y, Z, X , Y and Z be the unit vector of R ij , X , Y , Z , X , Y and Z in{O-XYZ}, respectively. Based on (3), it leads to: Let A i and B i (i = 1, 2, 3) be the vector of A i and B i in{O-XYZ}. And denote m A i (i = 1, 2, 3) be the vector of A i in {O X Y Z }. Based on the mechanism architecture, B i and m A i can be expressed as below: From (4), (5) and (6), the positional posture coupling relationships of m for the proposed PM can be obtained.
Furthermore, set m B R be formed by YXZ-type Euler rotations with α, β and λ as three Euler angles. Combined with (7), m B R can be specific expressed as below: where s θ denotes sin (θ ) and c θ denotes cos (θ ) with θ is α or λ, then (7) can be simplified as below: Equation (9) shows the explicit coupling relationships between the 6-dimensional positional posture of m. According to (9), α, λ and Z o can be considered as the independent kinematic parameters for this PM.
Complete positional posture analysis includes the forward and inverse analysis. When the actuators are set, the positional posture of the moving platform can be determined by the forward model. By contrast, the inverse analysis determines the required actuators variables from a given positional posture of the moving platform.
For the proposed PM, r i (i = 1, 2, 3) can be derived by the following formulas: Bring (5), (6), (8) and (9) into (10), it leads to: Equation (11) is the explicit inverse positional posture model of the proposed PM. According to (11), the inverse position solutions can be directly solved by the independent kinematic parameters (α, λ and Z o ) of m.
Eliminate Z o in (10), it leads to: Since (12) contains many trigonometric functions, to simplify expression, let t 1 = tan (α/2) and t 2 = tan (γ /2), then . Combined with (8) and (9), (12) can be simplified as below: where a i (i = 1, 2) and b i (i = 1, 2, 3, 4) only contain t 2 , which can be easily refined by MATLAB. Multiplying both sides of (13 a) by t 2 1 and t 4 1 , respectively, it leads to: Since (13a), (13b) and (14) form a system of four linearly independent equations with four variables t 1 , t 2 1 , t 4 1 and t 6 1 , it can be expressed in a matrix form as below: where Q is a coefficient matrix of the above linear equations about t 1 . Since t 1 must exist, there must be a nontrivial solution corresponding to (15). Based on the theory of linear algebra, it leads to: Equation (16) is a nonlinear equation with only regard to t 2 , which can be easily solved by MATLAB. After t 2 is VOLUME 10, 2022 solved from (16), t 1 can be solved from (13a). Then, the independent kinematic parameters α and γ corresponding to t 1 and t 2 can be obtained. Subsequently, the last one independent kinematic parameter Z o can be solved from (11).
Of course, the above forward positional posture model will get multiple solutions. However, combined with the computer-aided design variation geometry method [35], the unique solution can be determined.

B. INVERSE VELOCITY MODEL
Denote v r be the velocity vector of actuators and V be the velocity vector of m, respectively. For general n-DOF(n < 6) non-over-constraint PM with linear active leg, v r can be expressed as below [34]: where J is the Jacobian matrix which contains two submatrices, one is the traditional n×6 Jacobian matrix of limited-DOF PMs (named J α ), and the other one is the velocity constraints Jacobian matrix (named J v ). For this PM, since F pi (i = 1, 2) and T 1 do no work to m, it leads to: where C is the coordinate of C in {O-XYZ}. Since C is the intersection point of R 21 and R 24 , C should meet R 21 and R 24 line equations simultaneously. From (3), the line equation of R 24 can be expressed as below: Since the point on R 21 satisfies: x = 0 and y = E, substituting x = 0 and y = E into (19), C can be derived.
Combined (4) with (18), it leads to:  Then, combined (17) with (21), the inverse velocity model of the proposed PM is built: Since coupling relationships between the 6-dimensional positional posture of the moving platform are existed in limited-DOF PMs, there must exists velocity coupling relationships of the moving platform, which should be added to the inverse velocity model.
The relations between V and the velocity of independent kinematic parameters θ i (i = 1, . . . , n) for n-DOF(n < 6) PMs can be expressed as below [34]: where J 01 and J 02 are the linear and angular velocity decoupling Jacobian matrices, respectively. For this PM, from (8) and (9), it leads to: Likewise, from (8) and (9), it leads to: where: Combination of (22), (24) and (25) are the explicit inverse velocity model of the proposed PM. Whenα,λ andŻ 0 are given, V can be obtained according to (24) and (25). Subsequently, v r can be obtained according to (22).

C. INVERSE ACCELERATION MODEL
Denote a r be the acceleration vector of actuators and A be the acceleration vector of m. Deriving time from both sides of (17), a r can be expressed as below: where H is the Hessian matrix which contains two sub-matrices, one is the traditional Hessian matrix of limited-DOF PMs (named H α ), and the other one is the constraints Hessian matrix (named H v ). In order to solve H v , the derivatives of some vectors in (21) should be derived. For r 1 : Based on the vector algorithm, H v corresponding to r 1 can be expressed as below: For r 2 :ḟ In order to get an explicit expression of (29),θ 21 andḋ 2 need to be solved. Considering ω is a composite of all R joints in r 2 , it leads to: To eliminate unnecessary parameters, dot multiply both side of (30) by R 21 and R 24 , respectively, it leads to: Combined with (3) and (4), (31) can be reduced to the following form. Subsequently,θ 21 can be obtained: In addition, considering B 2 A 2 , B 2 C and CA 2 are a closedloop, it leads to: Simultaneously derived time from both sides of (33), and let v B 2 C and v CA 2 be the velocity value of B 2 C and CA 2 , respectively. It leads to: Meanwhile,ḋ 2 can be derived. Combined with (34), it leads to:ḋ Since (35) still contains uncertainties parameters v CA 2 , dot multiply both sides of (34) by X, it leads to: From (32), (35) and (36), (29) can be expressed. Then, based on the vector algorithm, H v corresponding to r 2 can be expressed as below (37), as shown at the bottom of the next page.
Combined with (26), (28) and (37), the inverse acceleration model of the proposed PM is established.
It is the same reason that there must exists acceleration coupling relationships of the moving platform, which should be added to the inverse acceleration model. Deriving time from both sides of (23), it leads to: where: For this PM, from (9) and (24), it leads to: Likewise, from (25), it leads to: where: Whenα,λ,Ż 0 ,α,λ andZ 0 are given, a and ε can be obtained according to (39) and (40). Subsequently, a r can be obtained according to (26), (28) and (37).

IV. STIFFNESS ANALYSIS OF THE RPU+UPU+SPU PM A. STATICS ANALYSIS AND CONSTRAINTS DECOMPOSITION
Utilizing the principle of static equilibrium, the relations between the active forces (F ai ) (i = 1, 2, 3), the constrained forces/torque (F pi /T 1 ) (i = 1, 2) and the loading force/torque (F/T ) can be expressed as below: To find the wrenches directly causing deformations, F p2 should be decomposed [26]. Based on the principle of force equivalence, F p2 can be equivalent to a force F s which is parallel with F p2 as well as passes through A 2 and a torque T 2 which is perpendicular with A 2 C and F p2 . F s and T 2 can be expressed as below: To find the deformations due to T i , decompose T i into two elements (i = 1, 2) [26], one is along r i (named T ip ), and the other one is perpendicular to r i (named T iq ), as shown in FIGURE 3. Let τ 2 , τ ip and τ iq be the unit vector of T 2 , T ip and T iq , respectively. Based on (4), it leads to: 23 (43) According to (43), the value of T ip and T iq (i = 1, 2) can be expressed as below: where s ip /s iq (i = 1, 2) is a coefficient representing the relationship between the constrained forces/torques and their components.
Though the above analysis, it leads to: where W is a coefficient matrix representing the relationship between the combination of the constrained wrenches and the active forces and the wrenches directly causing deformations.

B. DEFORMATION ANALYSIS
Suppose that three elastic r i (i = 1, 2, 3) elastically suspend m and all joints are a rigid body. Corresponding to different wrenches, deformations of r i can be analyzed based on the material mechanics. a) The forces directly causing deformations: F ai produces the longitudinal deformation δ ri along r i (i =  1, 2, 3). Let k ri be a coefficient mapping the relationships of δ ri and F ai , it leads to: where E is the modular of elasticity and S i is the area of r i . F p1 produces the flexibility deformations δ dp1 in r 1 and F s produces the flexibility deformations δ dp2 in r 2 . Let k p1 be a coefficient mapping the relationship of δ dp1 and F p1 , k p2 be a coefficient mapping the relationship of δ dp2 and F s , it leads to: where I is the moment inertia. b) The torques directly causing deformations: T ip produces the torsional deformation δ θ ip in r i (i = 1, 2). Let k ip be a coefficient mapping the relationships of δ θ ip and T ip , it leads to: where G is the shear modulus and I p is the polar moment of inertia. T iq produces the bending deformation δ θ iq in r i (i = 1, 2). Let k iq be a coefficient mapping the relationships of δ θ iq and T iq , it leads to: Though the above analysis, it leads to:  where K r is a coefficient matrix mapping the relationships between the deformations and the wrenches directly causing deformations.

C. STIFFNESS MODEL
Let δp = [δxδyδz] T , δ = [δ x δ y δ z ] T be the linear and angular deformations of m. According to the principle of virtue work, it leads to: Substituting (45) into (51), and combining with (41), it leads to: According to (52), the deformation relations can be derived as below: Furthermore, from (41), (45) and (53), it leads to: Substituting (54) into (53), it leads to: From (55), the stiffness model of the proposed PM is built: where K 6×6 and C 6×6 is the stiffness matrix and compliance matrix of this PM, respectively.

V. NUMERICAL EXAMPLES
The correctness of the previously established models is verified in this section.

A. FORWARD POSITIONAL POSTURE NUMERICAL EXAMPLE
For this example PM, the dimension parameters and the active parameters of each leg are shown in TABLE 1. According to TABLE 1,(15) can be written in a specific expression and Q can be obtained. Then t 2 can be solved form (16), the detailed values are listed in TABLE 2. In order to determine the acceptable analytic solution from TABLE 2, the simulative PM is created using the computer-aided design variation geometry method [35], as shown in FIGURE 4.  In the CAD software when given the identical settings of TABLE 1 to the simulative PM, the value of λ can be measured, which is in excellent agreement with the 5 th solution in TABLE 2. Subsequently, bring the 5 th solution t 2 = 0.1612 (corresponding to λ = 18.3146) into (9), (11) and (13a), other positional posture parameters can be solved. The detailed analytic values are shown in TABLE 3. Meanwhile, the simulative value of the positional posture parameters can be measured form the simulative PM, the   detailed simulative values are shown in TABLE 3. From  TABLE 3, it can be seen that the analytic values and the VOLUME 10, 2022   From an analytic point of view, when motions of the independent motion parameters are given, the initial positional posture of m can be calculated based on (9). The velocity v and ω of m can be calculated based on (24) and (25), FIGURE 5(a) and 6(a) shows the corresponding analytic values curves plotted by MATLAB.
The acceleration a and ε of m can be calculated based on (39) and (40), FIGURE 7(a) and 8(a) shows the corresponding analytic values curves plotted by MATLAB. Furthermore, according to the above expected motions of m, the initial length of actuators can be calculated based on (11). v ri (i = 1, 2, 3) can be calculated based on (22). And a ri (i = 1, 2, 3) can be calculated based on (26), (28) and (37).
From a simulative point of view, a simulative PM by MATLAB/SIMULINK is set, as shown in FIGURE 9. In the simulative PM, three simulation-driven modules are set on r i (i = 1, 2, 3), a velocity/acceleration sensor is set on m and the sensor values can display by scope. Appling identical v ri , a ri calculated above to the corresponding simulation-driven module, the scope displays the corresponding simulative values curves, as shown in FIGURE 5,6,7,and 8 (b).
The quantitative comparisons between the analytic values and the simulative values at t= 0s, t= 1.25s, t= 2.5s are shown in TABLE 5.
By comparing the curves in FIGURE 5,6,7 and 8 and the values in TABLE 5, it can be seen that the analytic values and the simulative values are in excellent agreement, which   verifies the correctness and precision of the kinematic model established in this study.

C. STIFFNESS NUMERICAL EXAMPLE
For this example PM, the dimension parameters, the initial independent motion parameters, the load force/torque and the mechanical parameters of materials are shown in TABLE 6.   Since the results of the FE simulated model depend on multiple key factors (e. g. dimensions and types, boundary constraints, solver, connection constraints), it is well known that the results of the FE model are approximate numerical results. Therefore, the analytic deformation values and the simulative deformation values in TABLE 6 are basically coincident, which is acceptable for stiffness analysis, and verifies the correctness of the stiffness model established in this study.
In future work, a RPU+UPU+SPU machine tools prototype will be designed and fabricated, and the experimental study will be performed to further validate the correctness of the kinematic and stiffness models built in this study.
The complete kinematic models for the RPU+UPU+SPU PM are built. Simulative PMs are set to validate the accuracy of the conducted kinematic models.
A 6 × 6 form Jacobian matrix is derived. And a 6 × 6 × 6 form Hessian matrix is derived. The velocity/ acceleration coupling relationships of the moving platform are derived to supplement the inverse kinematic models.
Both considering the active forces and the constraint wrenches, the stiffness model of the RPU+UPU+SPU PM are established.
The stiffness matrix and compliance matrix are derived. A FE simulative PM is built to validate the accuracy of the conducted stiffness model.

Symbol
Meaning R, P, U and S The revolute joint, the prismatic joint, the universal joint, and the sphere joint, M The DOF number of the mechanism n The number of links in the mechanism g The number of joints m i The DOF number of the i-th joint m 0 The passive DOF B 1 , B 2  The j-th R joint from the base platform to the moving platform in the i-th leg, the unit vector of R ij r i , δ i The length of the i-th leg, the unit vector of r i F p1 and T p The constrained force and torque in the 1 st leg F p2 The constrained force in the 2 nd leg f i and τ The unit vectors of F pi (i = 1,2) and T P C, C The The output velocity/acceleration vector J 01 and J 02 The linear and angular velocity decoupling Jacobian matrices H 1 and H 2 The linear and angular acceleration decoupling Hessian matrix of n-DOF PM.

F ai
The active force F/T Loading force/torque I p The polar moment of inertia F p1 , F s , T ip and T iq The force/torque directly casing deformations τ ip and τ iq The unit vector of T ip and T iq s A coefficient that represents the relationship between the constrained force/torque and its component δ ri The longitudinal deformation along r i k ri A coefficient for mapping the relationships of δ ri and F ai E The modular of elasticity S i The area of r i . δ dpi The flexibility deformation in r i k pi A coefficient for mapping the relationships of δ dpi and force I The moment inertia δ θip The torsional deformation about r i , k ip A coefficient for mapping the relationships of δ θ ip and torque G The shear modulus k iq A coefficient for mapping the relationships of δ θ iq and torque δp = [δxδyδz] T , δ = [δ x δ y δ z ] T The linear and angular deformations of the moving platform W A coefficient matrix that represents the relationship between the constrained force/torque and the force/torque directly casing deformations K r A coefficient matrix for mapping the relationships of deformation and force/torque K 6×6 and C 6×6 The stiffness matrix and compliance matrix