Elastodynamic modeling and parameter sensitivity analysis of a parallel manipulator with articulated traveling plate

This paper deals with the elastodynamic modeling and parameter sensitivity analysis of a parallel manipulator with articulated traveling plate (PM-ATP) for assembling large components in aviation and aerospace. In the elastodynamic modeling, the PM-ATP is divided into four levels, i.e., element, part, substructure, and the whole mechanism. Herein, three substructures, including translation, bar, and ATP, are categorized according to the composition of the PM-ATP. Based on the kineto-elastodynamic (KED) method, differential motion equations of lower levels are formulated and assembled to build the elastodynamic model of the upper level. Degrees of freedom (DoFs) at connecting nodes of parts and deformation compatibility conditions of substructures are considered in the assembling. The proposed layer-by-layer method makes the modeling process more explicit, especially for the ATP having complex structures and multiple joints. Simulations by finite element software and experiments by dynamic testing system are carried out to verify the natural frequencies of the PM-ATP, which show consistency with the results from the analytical model. In the parameter sensitivity analysis, response surface method (RSM) is applied to formulate the surrogate model between the elastic dynamic performances and parameters. On this basis, differentiation of performance reliability to the parameter mean value and standard variance are adopted as the sensitivity indices, from which the main parameters that greatly affect the elastic dynamic performances can be selected as the design variables. The present works are necessary preparations for future optimal design. They can also provide reference for the analysis and evaluation of other PM-ATPs.


Introduction
Parallel manipulator with articulated traveling plate (PM-ATP) is one of the most well-recognized mechanisms in the research community of parallel manipulators [1,2]. The common parallel manipulator is composed of one fixed base, one moving platform, and several kinematic chains linking to them. The moving platform is usually a single plate. Alternatively, the articulated traveling plate (ATP) is a special type of moving platform consisting of two or more in-parts and one out-part [3]. Besides the mobility provided by the kinematic chains, PM-ATP gains extra motions from the ATP. Therefore, PM-ATPs are more flexible in terms of motion capability. One typical example of PM-ATPs is the parallel manipulator with Schönflies motion (three translations and one rotation around the vertical axis, i.e., 3T1R) whose rotation is provided by the relative translation of the two in-parts. The well-known Par4 [4], I4 [5], and the four degree of freedom (DoF) parallel robot [6] belong to this group of PM-ATPs. The extra rotation from the in-parts is up to 720 degree, making the PM-ATP attractive to posture changing of disorder products. In practice, these PM-ATPs have been successfully applied for the high-speed pick-andplace in food packaging, medicine, and semiconductor manufacturing.
Inspired by the successful applications of parallel manipulators with Schönflies motion, more and more attention has been drawn to the investigation of PM-ATPs. Referring to the industrialization of conventional parallel manipulators in the sequence of topology innovation [7,8], optimal design [9][10][11][12][13], calibration, and control [14][15][16][17][18], the developments of PM-ATPs start from the topology synthesis. The aim is to invent new PM-ATPs for wider industrial scenarios [19,20].
In this regard, Sun [3] discussed the kinematic constraints within the ATP and proposed a group of PM-ATPs that can be applied as tracking mechanism, docking equipment, or machine tools. The parameterized topological models were further analyzed [21][22][23]. By filling in the gap between finite and instantaneous screw theory, Sun [24,25] succeeded in connecting topology analysis to the following performance analysis and even optimal design, which is a milestone in the topology synthesis of parallel manipulators.
A PM-ATP (details shown in Section 2) is then selected from the topology synthesis and used as a pose-adjusting mechanism for assembling large components in aviation and aerospace. The next problem for the development of the PM-ATP is the optimal design that builds an actual mechanism from certain topological structure. The concerned performances would be optimized by adjusting the structural parameters which are regarded as design variables. Therefore, the two essential elements for the optimal design are the performance and the structural parameters.
The commonly concerned performance indices of parallel mechanisms include workspace [26], singularity avoidance [27,28], stiffness [29,30], and dynamic [31][32][33]. Since the studied PM-ATP is targeted for assembling large components, elastodynamic performance [34,35] catering on large loadcarrying, lightweight structure and small deflections is of importance. Performance indices such as natural frequency or elastic deformations can be adopted as objectives in the optimal design, which require for the mapping model between the elastodynamic performances and the structural parameters.
The existing elastodynamic modeling methods are mainly for the high-speed pick-and-place PM-ATPs whose links are made of light weight material. The elastic deformations of the links are coupled with the pick-and-place motions, resulting in an ongoing trend of applying flexible multibody dynamic methods to model the elastodynamic performances [36][37][38]. Although the obtained elastodynamic models are with high accuracy, the modeling process is computationally expensive due to the nonlinear couplings between link deformations and rigid motions. For the rigid structures moving in a low speed, however, elastic deformations of parts are much smaller than the rigid body motions. The deformations are assumed not to affect the mechanism motions and the kineto-elastodynamic (KED) method can be applied [39,40]. In the KED framework, the motions of PM-ATP are firstly analyzed by the rigid body kinematics, and then the elastic deformations computed by the structural dynamics at each instantaneous moment are added. The modeling procedure is greatly simplified while the effectiveness in describing the elastodynamic performance of the whole mechanism can be kept.
The KED method is applicable under the assumption that the parts are relatively rigid and the mechanism moves slowly. Since the pose-adjusting motions of the studied PM-ATP are relatively slow and the parts are designed towards high stiffness, the KED method can be adopted. Two challenges need to be addressed in the elastodynamic modeling of the studied PM-ATP by the KED method. (1) The parts are usually with irregular shapes, increasing difficulty in analytically computing the dynamic behavior. (2) The ATP contains complex structures and multiple joints, complicating the whole system.
In order to address the modeling difficulties, the studied PM-ATP is divided into four levels, i.e., element, part, substructure, and the whole mechanism. By applying the element as basic unit, the elastic deformation of irregular parts can be captured. The elastodynamic model of each level is built by the KED method and the model of upper level is assembled by the model of lower level. The obtained elastodynamic model of the PM-ATP will be verified by simulations in finite element software and elastic dynamic experiments in the following sections. Herein, the relation of structural parameters and elastodynamic performance are the major concern in the modeling. Noises [41,42] that impose influence on the dynamic performance in practical application are not included. They are regarded as system disturbances and solved by in the controller development after the optimal design of the PM-ATP.
However, the optimal design is still challenging if the elastodynamic model is directly applied. Large amounts of structural parameters are involved because of the irregular parts, the complex structures, and the compositions of the PM-ATP. Parameter sensitivity analysis is usually implemented to exclude the trivial parameters and simplify the model [43]. By analyzing the change of the output performances when varying the input parameters, parameter sensitivity identifies the effects of parameters to the performances. Parameters with high sensitivity impose more influence to the performance and should be chosen as design variables while the parameters with low sensitivity can be eliminated in the optimal design. Current parameter sensitivity methods mainly compute performance changes by randomly changing the values of one parameter at a time, in which the rest of the parameters remain unchanged [44,45]. The coupling effects of multiple parameters are ignored in these methods. In order to efficiently and effectively analyze parameter sensitivity, more comprehensive parameter sensitivity indices are required. They should (1) consider the possible coupling effects among parameters, and (2) evaluate the change of each parameter via statistic technique, for instance mean value and standard variance.
In summary, this paper focuses on the preparations for the optimal design of a PM-ATP, i.e., the elastodynamic modeling and the parameter sensitivity analysis. The difficulties of this work are resulted from the complicated composition and the substantial parameters. To illustrate the adopted methods, this paper is organized as follows. Section 2 briefly introduces the PM-ATP and carries out the inverse kinematics analysis. Based on the KED method, elastodynamic modeling of the PM-ATP is implemented in Section 3. Translation, bar, and ATP substructures are assigned to assemble the dynamic model of the whole mechanism. In Section 4, simulation and experiment are conducted to verify the elastodynamic model. Section 5 is dedicated to the parameter sensitivity analysis, from which the main parameters are identified and selected as design variables in future optimal design. Conclusions are drawn in Section 6.

Mechanism description and inverse kinematics
The studied PM-ATP is named as PaQuad PM (Fig. 1a). The PaQuad PM is composed of a fixed base, an ATP, and four identical PRS limbs. Herein, P, R, and S denote actuated prismatic, revolute, and spherical joint. The PRS limbs connect to the fixed base and the ATP by P joint and S joint, respectively. The ATP consists of in-part 1, in-part 2, and out-part. The inpart 1 links to the out-part through helical (H) joint, whereas the in-part 2 joins to out-part by R joint. The axes of the H joint and R joint are collinear. According to the mobility analysis, the ATP has one translational and two rotational capabilities provided by the PRS limbs. Additionally, the relative translation of in-part 1 and in-part 2 results in the rotation of H joint, which adds an extra rotation to the out-part. Hence, the PaQuad PM has one translation and three rotations.
In order to formulate the kinematic model of the PaQuad PM, some denotations and coordinate frames are defined as shown in Fig. 1b. Point A i is assigned to the connecting point of the ith (i = 1, 2, ⋯, 4) PRS limb and the fixed base. The fixed base is defined by a circle whose center is point O and radius is a. Points B i , C i , and D i denote centers of P joint, R joint, and S joint, respectively. The lengths of in-part 1 (D 1 D 3 ) and in-part 2 (D 2 D 4 ) are both 2b and the vertical distance between them is e. The traveling distance of P joint and the length of bar are represented by q i and l. A fixed reference frame O − xyz is assigned to point O. The x-axis is collinear with OB 2 and the z-axis is vertical to the fixed base. A moving reference frame O ′ − uvw is attached to the point O ′ on the outpart. Its w-axis points to the same direction as H joint and the u-axis is parallel to D 4 D 2 at home position.
The rotation matrix of frame O ′ − uvw with respect to frame O − xyz can be computed by where α, β, and γ are the three Euler angles. c and s denote cosine and sine, respectively. Considering the projection of point D i on the O ′ − uv plane, the position vector of point D i in frame O ′ − uvw can be expressed as where γ 1 denotes additional rotation angle produced by Euler angles α and β. Point D i can also be described in frame O − xyz as where The point D i moves within the plane spanned by s 1, i and s 3, i due to the limitation of the R joint, which can be mathematically described by where (1) and Eq. (2) into Eq. (3) yields where e ¼ P h Âγ 2π þ e 0 . P h denotes screw pitch, e = e 0 when PaQuad PM is at home position. d 0 represents distance between point O ′ and D 2 D 4 . The closed-loop equation can be formulated as For the inverse kinematic of the PaQuad PM, the z value of point O ′ and the three Euler angles (α, β, γ) are the known parameters. By solving Eq. (1) to Eq. (5), the traveling distance of the P joint is obtained as follows where

Elastodynamic modeling
As has been mentioned in Section 1, the KED method is adopted to formulate the elastodynamic model of the PaQuad PM. Based on the inverse kinematics, any configuration within workspace can be computed, with which the elastodynamic modeling is implemented. The PaQuad PM is divided into four levels, i.e., element, part, substructure, and the whole mechanism. Euler-Bernoulli beam is applied to be the basic element, whose dynamic performance can be analytically formulated.
The elastic dynamics of parts are denoted by the elements and then assembled to form the elastodynamic model of substructures. Finally, the model of the whole mechanism is established by the models of substructures. During the assembling process, displacements at connecting points of parts and deformation compatibility conditions of the substructures are concerned. According to the KED method, some basic assumptions are made as follows. (1) Figure 2 shows a spatial beam element whose two nodes are E 1 and E 2 . An element coordinate frame E 1 −xyz is assigned to point E 1 . Elastic deformation of any point on the beam can be expressed as

Differential motion equation of beam element
U(t), V(t), and W(t) are the linear deformations along x, y, and z -axis while θ x (t), θ y (t), and θ z (t) are the angular deformations. The details are referred to [39]. N is a type function matrix [39], and u ¼ u 1 u 2 ⋯ u 12 ½ T is the elastic deformations of node E 1 and E 2 .
The beam element generates tension, compression, bending, and torsional deformations if an arbitrary force is applied. Elastic energy of the beam is thus computed by where E, A, and G are Young's modulus, area of cross section, and shear modulus. I x , I y , and I z denote the polar moment of inertia about each axis. U ′ and θ 0 x are the first-order differentiation of U and θ x , and V ″ and W ″ are the second-order differentiation of V and W.
In addition, gravity potential energy is given by where ρ, g, and r p2 are density, gravity acceleration, and vector Moreover, kinetic energy can be expressed as where u ai ¼ N u r þ u ð Þ¼Nu a . Herein, u ai is the absolute velocity of point E i . u r and u denote velocities of rigid body motion and elastic deformation.
Hence, Lagrange equation is formulated as where E p = E p1 + E p2 , f is the vector of external forces. Differential motion equation of the beam can be formulated by substituting Eq. (8) to Eq. (10) into Eq. (11) as follows.
where U e is the generalized coordinates of nodes. M e , K e , F e , and Q e are the mass matrix, stiffness matrix, and vector of external and internal forces in frame O − xyz.
herein T is the transformation matrix of frame E 1 −xyz with respect to frame O − xyz.

Differential motion equation of substructures
The differential motion equation of the beam is applied to assemble the elastic dynamic of the parts by considering DoFs of the connecting points. Then, the parts will be used to construct the substructure. The PaQuad PM is divided into three substructures, i.e., translation, bar, and ATP. The first two substructures are from the PRS limb, where the former is the main body of the P joint and the latter is the connecting structure of the R and S joints. The effects of P, R, and S joints will be taken into account in the mechanism level. Unlike the other two substructures, ATP contains internal joints (H and R joints) whose influences need to be addressed in the substructure level.

Translation substructure
The P joint is composed of screw pair and guide slider, as is shown in Fig. 3. Ball screw and slider are the major components to deform; thus, they are assumed to be elastic and represented by two and three beam elements, respectively. The elements are named from iE1 to iE5 and the seven nodes are denoted by iN1 to iN7. Element coordinate frames are firstly established for establishing the differential motion equation of elements. Frame iN1 − x i y i z i is assigned to node iN1. The x i -axis is collinear with the rotational axis of ball screw, the y i -axis is parallel to the x-axis in frame O − xyz, and z i -axis satisfies right hand rule. Frames of element iE2, iE3, and iE4 are parallel to frame iN1 − x i y i z i . Frame iN7 − xi y i z i is established at node iN7, whose x i -axis points to the axis of element iE5 and z i -axis is in the same direction as z-axis of O − xyz. Based on Eq. (12), the differential motion equation of ball screw can be computed as where Similarly, the differential motion equation of slider is where Next, the boundary conditions and relations of connecting points are analyzed. The elastic deformations of node iN1 are restricted because one end of the ball screw is fixed to the servo motor. The other end is linked to the fixed base by bearings. Thus, the five elastic deformations of node iN3 are zeros except for the rotation about the screw axis. Node iN2 of the ball screw and node iN5 of the slider are connected by H joint, hence where p h is the pitch of ball screw. Finally, the differential motion equation of translation substructure is obtained by Eq. (13) to Eq. (16) as where U ip is generalized coordinates of translation substructure. U ip , M ip , K ip , F ip , and Q ip are shown under Fig. 3, in which B ipj is referred to Appendix.

Bar substructure
According to the features of the bar, ten elements with nine nodes are assigned to the bar, as is shown in Fig. 4. The element coordinate frames iNj − x i y i z i (i = 5, ⋯, 8; j = 1, ⋯, 10) are defined, where the x i -axis is along the length of the elements. Differential motion equation of the bar is expressed as where and Q ip are the mass matrix, stiffness matrix, and vector of external and internal forces, herein J iEj (j = 1, 2, ⋯, 10) is shown in Appendix.

ATP substructure
The ATP is shown in Fig. 5. Concerning the effects of the internal joints, the procedure for formulating differential motion equation of ATP is summarized as follows: (1) apply beam element to describe the major features of each component and establish element coordinate frames, (2) formulate differential equation of each component separately, (3) analyze the relations of connecting modes according to the assembling conditions, and (4) assemble to get the differential equation of ATP by step (2) and (3). First of all, elements are assigned to the parts. The in-parts are represented by four beam elements (R1E1, R2E1, R3E1, and R4E1). The two cylindrical structures are designed to enhance translational capability of the H joint and assessed by concentrated masses at their centers of mass (R5 and R6). The screw of the H joint is denoted by two elements (R7E1 and R7E2). The out-part is regarded as concentrated mass (R8). Then element coordinate frames are defined. For the elements R1E1, R2E1, R3E1, and R4E1, frame D i − x Ri y Ri z Ri . The x Ri -axis is collinear with direction of element length and z Ri -axis is perpendicular to the plane of in-parts. For the elements R5, R6, and R8, frame O ′ − x Ri y Ri z Ri (i = 5, 6)is defined at point O ′ . The x Ri -axis is along the axis of screw, and z R5 -axis and z R6 -axis are parallel to x R4 -axis and x R1 -axis.
Next, differential motion equation of each component can be formulated by the elements as where U Ri is the generalized coordinates of components; M Ri , K Ri , F Ri , and Q Ri are the mass matrix, stiffness matrix, and vector of external and internal forces. Differential motion equation of R7 and R8 is formulated The relations at connecting nodes of each part are then assessed by the assembling conditions as shown in Fig. 5. Node R8 and node R7N3 have the same generalized coordinates since the out-part and the screw is fixed. Node R8 and node R5 are linked by R joint; thus, the five coordinates are the same except for the angular deformation about rotational axis. Cylindrical joint is formed between node R5 and node R6, and the coordinates about/along the joint axis are different. Node R6 and node R7 is connected by H joint. Node R6 is attached to node R1 and R3 rigidly, while node R5 is fixed to node R2 and R4. These assembling conditions can be mathematically expressed as where R 56 is the rotation matrix of frame O ′ − x R5 y R5 z R5 with respect to frame O ′ − x R6 y R6 z R6 . R 61 and R 63 are rotation matrices of frame O ′ − x R6 y R6 z R6 with respect to frame D 1 − x R1 y R1 z R1 and frame D 3 − x R3 y R3 z R3 . R 52 and R 54 are the rotation matrices of frame O ′ − x R5 y R5 z R5 with respect to frame D 2 − x R2 y R2 z R2 and frame D 4 − x R4 y R4 z R4 . S(a 1 ) and S(a 3 ) are the skew matrices of node R6 to node R1N2 and node R3N2 in the frame O ′ − x R6 y R6 z R6 . S(a 2 ) and S(a 4 ) are skew matrices of node R5 to node R2N2 and node R4N2 in the frame O ′ − x R5 y R5 z R5 . Finally, the differential motion equation of ATP can be formulated from the equations of parts and relations of connecting nodes as follows.
where U 9p is a 40 × 1 vector representing generalized coordinates.

Dynamic model of PaQuad PM
With the differential motion equations of substructures available at hand, the elastodynamic model of the whole mechanism is assembled by the deformation compatibility conditions among substructures. Node iN7 of translation substructure connects to node iN8 and iN9 of bar substructure by R joints. Hence, the other five generalized coordinates of these nodes are the same except for the rotational deformations about the axis of R joint.
where k = 0, 2. Node iN2 (i = 5, ⋯8) links to the in-parts through S joints. Considering the stiffness of the S joints, the generalized coordinates of connecting nodes between bar structure and ATP are expressed as Therefore, the elastodynamic model of the PaQuad PM is obtained by assembling differential motion equations of substructures as follows where M, K, and F are the mass matrix, stiffness matrix, and external forces of PaQuad PM.
4 Natural frequencies

Case study
The virtual and physical prototypes of the PaQuad PM were built, which can be used to verify the elastodynamic model formulated in Section 3. The natural frequency of the PaQuad PM can be computed by the elastodynamic model shown in Eq. (34) as where ω i (i = 1, 2, ⋯, n) is the natural frequencies.
Dimensional parameters of the PaQuad PM are shown in Table 1. The major material for most components is 45# steel, whose density is 7.8 × 10 3 kg/m 3 , Young's modulus is 2.2 × 10 11 Pa and shear modulus is 7.938 × 10 10 Pa. The mass of inpart 1 is 13.13 kg. Its moment of inertia about x-, y-, and z-axis is 0.0598, 0.0598, and 0.0572 kg ⋅ m 2 , respectively. The mass of in-part 2 is 20.46 kg, whose moment of inertia is 0.146, 0.146, and 0.278 kg ⋅ m 2 . The mass of out-part is 34.27 kg, the moment of inertia is 0.348, 0.348, and 0.681 kg ⋅ m 2 .
By applying the parameter value to the elastodynamic model, the distribution of natural frequencies within workspace is shown in Fig. 6. The first frequency is symmetrical to the plane β = 0 ∘ and the maximum value is at β = 0 ∘ , then it decreases with the increasing of β. With the increment of α, the first frequency drops when β is fixed. The second frequency is symmetrical to the axis α = β = 0 ∘ . It decreases with the changes of α and β. The peak value is at the symmetrical axis.
The third frequency monotonously climbs up as α increases and obtains the maximum value at α = 10 ∘ when β keeps the same. However, the third frequency decreases if α increases to α = 30 ∘ , then it increases again. The fourth frequency is plane symmetrical. Maximum and minimum values are mainly on the boundary of the workspace. The former are at α = 0 ∘ , β = ± 40 ∘ , and α = ± 40 ∘ , β = 0 ∘ , while the latter are at α = ± 40 ∘ and β = ± 40 ∘ . The fifth frequency is axial symmetrical to α = β = 0 ∘ , where it gets the minimum value. Distribution of the sixth frequency is similar to the first frequency but the change is sharper.
In summary, the natural frequencies change versus configurations and they show plane-symmetrical features due to the plane-symmetrical structure of the PaQuad PM.

Simulation and experiment
Simulations on the virtual prototype by FEA software are applied to verify the elastodynamic model of PaQuad PM. SAMCEF from SEMTECH Inc. is chosen to analyze the natural frequencies of eight typical configurations within workspace. The simulation is implemented as follows. Result. (6) Choose the second configuration and proceed to step (1) to step (5). Repeat until PaQuad PM under all eight configurations is simulated.
The simulation results are shown in Table 2. The changing tendencies of the first to sixth frequencies are similar for both theoretical analysis and simulations. Generally, the simulated frequencies are smaller because non-standard features, such as  Elastic dynamic experiments are also carried out as shown in Fig. 7. Measuring points are firstly assigned to different parts of PaQuad PM, on which acceleration sensors are attached. Hammer is then applied to excite the PaQuad PM at the end reference point. Through collecting exciting forces and response signals from the sensors, the modes are fitted and analyzed by LMS dynamic testing system (including SCSASA III data collecting hardware and LMS Test. Lab software from SIEMENS). Finally, the frequencies of PaQuad PM are obtained. The experimental procedure can be summarized as follows.
(1) Geometrical modeling. In the Geometry interface, the overall and local coordinate frames (fixed and element frames) are defined according to Section 3. Measuring points are assigned in the local frames based on the structures. By connecting the measuring points in the overall frame, the geometrical model of PaQuad PM is obtained. There are 68 measuring points in total. (2) Sensor setting. In the Channel Setup interface, channel 1 is assigned to measure force. The actual sensitivity is 0.2838 mv/N. Channels 2 to 4 are chosen to collect signals from acceleration sensors in three directions. Their sensitivity are 9.822, 9.99, and 9.864 mv/(m/s 2 ), respectively.  (3) Excitation setting. In the Impact Scope interface, the excitation of hammer is set as free form mode with 200 Hz bandwidth. The pre-excitation time and signal setups are also defined. The experimental results are shown in Table 2 and Fig. 8. Comparing with the analytical models and simulations, experimental frequencies are the smallest due to (1) detail features such as chamfer and groove are included in the physical prototype while they are ignored in the other two cases; (2) contacts between parts, especially within joints are not ideal in experiments; and (3) there are external influences like measurement noise, non-rigid fixed base, or low bearing preload.
The modals of the PaQuad PM are further analyzed to check the consistency of the analytical, simulation, and experiment results. Taking pose I in Table 2 as an example, the first and second modals are the relative vibrations of opposite PRS limbs along xand y-axis. The third and fourth modals are the torsions of opposite PRS limbs and the connected in-parts. The fifth and sixth modals are the vibrations along z-axis. The first to sixth modals have the same rules.
In summary, frequencies from the analytical model, simulation, and experiment are close, and the changing tendencies are similar. In addition, the modal of the eight typical poses are the same. The elastodynamic model is validated. The parametric elastodynamic model can be applied to the re-design of the PaQuad PM under specific requirements from different application scenarios.  The complicated composition and structures of the PaQuad PM lead to large amounts of parameters. In the re-design or the optimization process, however, it is not necessary to optimize all the parameters since some parameters impose little effects to the elastodynamic performance. Besides, the substantial parameters increase the difficulties of the optimization. Based on the engineering experience, the parameters are scaled down as shown in Table 3. Parameter sensitivity analysis is then implemented to categorize the main or subordinate parameters according to their effects to the natural frequencies. In our previous work, we proposed a parameter sensitivity analysis method by response surface method (RSM) based model, parameter mean value, and variancebased indices [44]. The RSM method is to build the surrogate model between parameters and the performances. Based on the explicit mapping model, the parameter mean value and variance-based indices are applied as the criteria for evaluating the parameter sensitivity. The merits of this method are twofold.
(1) Coupling contributions among different parameters are included in the analytical expression of RSM model.
(2) Comprehensive evaluation is achieved by considering the statistic features of parameters. According to this method, parameter sensitivity analysis of PaQuad PM can be divided into two parts: formulation and assessment of RSM model, computation, and analysis of sensitivity indices.

RSM model
Based on design of experiment (DoE), RSM employs a set of experiments to establish the polynomial surface function for mapping the relations between parameters and targeted performances [46]. Herein, experiments are performed by the elastodynamic model built in Section 3. The experimental setup, including different combinations of parameters, number of experimental sets, is decided by the DoE strategy. Due to the amount of parameters and uncertain order of polynomial functions, the Latin hypercube design (LHD) strategy [47] is chosen since it is capable of dealing computation-intensive problems with limited number of experiments. After setting up by LHD strategy and implementing the experiments by the elastodynamic model, the linear, quadratic, cubic, and quartic RSM functions can be formulated as Þ T is the vector of parameters. a 0 , b i , c i , d ij , e i , and f i are the estimated coefficients obtained from the least square regression. x i x j is the interaction of any two parameters. x 2 i , x 3 i , and x 4 i are the second, third, and fourth nonlinearity of x.
Accuracy assessment is then carried out to verify the RSM models. Additional parameter sets are randomly generated. The errors between the actual responses and the results from RSM models are compared via four metrics [48], i.e., R square (RS), relative average absolute error (RAAE), relative maximum absolute error (RMAE), and root mean square error (RMSE), as follows where y i is the natural frequency obtained from analytical model,ŷ i is the value computed by RSM model, and y is the mean value of y i . In these accuracy metrics, RS, RAAE, and RMSE evaluate the overall accuracy of RSM models within the parameter ranges while RMAE shows the maximum error. Through the simultaneous consideration of global and worst accuracy, RSM model that has smaller RS, RAAE, and RMSE and larger RS are selected as the surrogate model for parameter sensitivity analysis.
For the PaQuad PM, the parameter values shown in Table 1 are regarded as the baseline and the range of parameters is set as ±10%. Since the number of involving parameters is 20, the  Table 4). The errors of quadratic and cubic RSM models are smaller than the linear and quartic models. RS values of quadratic and cubic RSM models are the same. The RAAE and RMSE values of quadratic RSM model are lower than cubic model, but the RMAE value is higher. From the comparison, the quadratic RSM model is finally chosen.

Parameter sensitivity indices
Based on the explicit RSM model, parameter sensitivity indices are defined by the differentiation of performance reliability to the parameter mean value and variance. Performance reliability describes the probability of the studied PM achieving target performance with given range of parameters. It is expressed as where P{g(x) > 0} denotes the probability of g(x) > 0. And g(x) is the subtraction of RSM model and the allowable values. Y = [g(x) − μ r ]/σ r is a random variable determined by β r = μ r /σ r . Herein, Y ∈ [0, 1] is subjected to normal distribution. μ r and σ r are the mean value and variance of g(x) > 0. Taking differentiation of Eq. (43) to the mean value and covariance of parameters yields where μ j and σ j (j = 1, 2, ⋯, 20) are mean value and standard variance of jth parameter (see Table 3).
Considering both parameter mean value and variance, a global sensitivity index is defined as where κ j ¼ . κ max and σ max are the maximum sensitivity of mean value and variance among all the parameters. Based on Eq. (37) to Eq. (39), the parameter sensitivity of PaQuad PM is computed and shown in Table 5 and Fig. 9. In summary, ATP imposes great effects to the elastodynamic Fig. 9 Global sensitivity and parameter sensitivity to mean value and variance Table 5 Parameter sensitivity of first-order natural frequency performance, especially central screw and in-part 2 that directly link to the end-effector. Comparatively, PRS limbs are more rigid and hence have little influence to the natural frequency. For the sensitivity to mean values, diameter of central screw in H joint (d hj ), parameters of in-part 2 (l ip2 , h ip2 , and b ip2 ), diameter of S joint (d sj ), and width of fixed bar (b b ) have larger values than the rest of the parameters. Herein, h sd2 is the maximum, indicating it has the most significant influence to the natural frequency.
For the sensitivity to parameter variance, the top parameters are h ip1 , h ip2 , d sj , l sj , d hj , and b b . Considering both parameter mean value and variance, global indices show that b b , d sj , l ip2 , h ip2 , b ip2 , and d hj are the major parameters that would greatly affect the natural frequencies. By increasing the values of these parameters, the resulted natural frequency is expected to increase. Therefore, they can be chosen as the design parameters in the future optimization.

Conclusions
Aiming at the optimal design of a PM-ATP named as PaQuad PM, this paper carries out the elastodynamic modeling by KED method and parameter sensitivity analysis by RSM method and reliability sensitivity indices. Conclusions are drawn as follows.
(1) The PaQuad PM is divided into four levels, i.e., beam element, parts, substructure, and the whole mechanism. Differential motion equation of lower level is assembled to formulate the model of upper level. Through this layerby-layer modeling method, multiple parts and joints within ATP can be considered explicitly. The proposed elastodynamic modeling method can also be applied to other types of PMs. (2) Natural frequencies are regarded as the dynamic performance indices. Simulation by finite element software and experiments by LMS dynamic test system are implemented to verify the elastodynamic model of PaQuad PM. Although the results from the three methods are slightly different, the changing tendency and modal are the same. The analytical model is confirmed.
(3) RSM models are investigated to get explicit polynomial functions between natural frequency and parameters of PaQuad PM. The coupling effects among parameters can be evaluated. Parameter sensitivity indices are defined based on the differentiation of performance reliability to the parameter mean value and standard variance. The statistic features of parameters are considered; hence, the obtained parameter sensitivity is comprehensive.
Elastodynamic modeling is the basis of formulating objective functions while the parameter sensitivity analysis determines the design variables. They are both necessary preparations for the optimal design which is the future work for the development of the PaQuad PM.
The connecting relations of nodes within ATP can be calculated as ;  ;