Solving the dynamic equations of a 3-PRS Parallel Manipulator for efficient model-based designs

Introduction of parallel manipulator systems for different applications areas has influenced many researchers to develop techniques for obtaining accurate and computational efficient inverse dynamic models. Some subject areas make use of these models, such as, optimal design, parameter identification, model based control and even actuation redundancy approaches. In this context, by revisiting some of the current computationally-efficient solutions for obtaining the inverse dynamic model of parallel manipulators, this paper compares three different methods for inverse dynamic modelling of a general, lower mobility, 3-PRS parallel manipulator. The first method obtains the inverse dynamic model by describing the manipulator as three open kinematic chains. Then, vector-loop closure constraints are introduced for obtaining the relationship between the dynamics of the open kinematic chains (such as a serial robot) and the closed chains (such as a parallel robot). The second method exploits certain characteristics of parallel manipulators such that the platform and the links are considered as independent subsystems. The proposed third method is similar to the second method but it uses a different Jacobian matrix formulation in order to reduce computational complexity. Analysis of these numerical formulations will provide fundamental software support for efficient model-based designs. In addition, computational cost reduction presented in this paper can also be an effective guideline for optimal design of this type of manipulator and for real-time embedded control.


Introduction
Seminal research in Parallel Manipulators (PMs) described architectures of 6 degrees of freedom (DOF) which are mainly used to perform industrial tasks.Nevertheless, not all applications (e.g., commercial, space exploration, entertainment or even industrial) require full 6-DOF capabilities, thus, cost-effective PMs with less than 6 DOF (i.e., lowermobility) have been developed.One such architecture is the 3-PRS manipulator which has a platform and a fixed base connected through three identical sets of links and joints (i.e., legs).Each leg has a slider attached to the base by an actuated prismatic joint (P), a coupler connected to the slider by a passive rotational joint (R) and to the platform by a passive spherical joint (S).The 3-PRS manipulator was first de-scribed in Carretero et al. (2000b) for a telescope application, and then it was proposed as a machining centre in Fan et al. (2003) and for a medical application in Merlet (2001).The kinematics and workspace analyses of the manipulator have been extensively studied in Carretero et al. (2000a), Carretero et al. (2000b), Tsai et al. (2002), Li and Xu (2007), to name a few.On the other hand, despite the fact that inverse dynamic modelling is essential for optimal design, parameter identification (Mata et al., 2008), model-based control (Díaz-Rodríguez et al., 2013), and internal redundancy (Parsa et al., 2013) amongst others, few papers have focused on revisiting and comparing computationally-inexpensive methods in order to obtain the dynamic models that are cost-effective for real-time applications.
Published by Copernicus Publications.
Lagrangian formulations allowed to develop the inverse dynamics model of the 3-PRS manipulator (Li and Xu, 2004).The formulation uses the Lagrange multiplier to include the constraints forces that lead to a modelling approach not only intricate but also computationally complex.Li and Xu (2004) applied the Principle of Virtual Work (PVW), but they simplify the dynamics of the coupler link by dividing its mass into two portions located at its extremes.Tsai and Yuan (2010) solved the inverse dynamic model along with the reaction forces through a special decomposition of the reaction forces at the joints that connect the leg with the platform.A similar approach was used in Yuan and Tsai (2014) for solving direct dynamics including friction effect.However, the later approach considers the calculation of reaction forces which are may be needed for structural design of a manipulator but its computation increases computational complexity which is unnecessary for parameter identification or model-based control.Staicu (2012) analyses and compares the power consumption of the 3-PRS vs. the 3-PRS configuration using the PVW with recursive modelling.The method obtains the Jacobian by differentiating the vector loop equation.Carbonari et al. (2013) solved the inverse dynamics of a 3-DOF parallel manipulator via screw theory and the PVW.On the other hand, the 3-RPS manipulator presents similar characteristics to the 3-PRS manipulator, in this respect, Mata et al. (2008) implement recursive velocity equations used in serial manipulator analysis to find the Jacobian of the manipulator for the inverse dynamic modelling, and Ibrahim and Khalil (2007) exploit architectural characteristics of the 3-RPS to give a closed form solution for the inverse and direct dynamics modelling.
The inherent complexity of the dynamic models lies on the way the system is modelled and how the Jacobian matrix is put forward.In this context, this paper compares the computational number of operation of three formulations for inverse dynamic modelling of a 3-PRS.The first formulation applies the general solution of PMs dynamic modelling proposed in Khalil and Ibrahim (2007).The second method considers the manipulator as a set of open kinematic chains and finds the Jacobian in joint space coordinates by taking into account the vector loop constraints at the split joints (Mata et al., 2008).The third method relies on the modelling approach originally presented in Li and Xu (2004).
The ultimate goal of contrasting these numerical formulations is focused on supporting the implementation of emerging model-based designs which not only depends of the discrete inverse dynamic model but also the numerical finite realization in a given computational platform (Williamson, 1991).In fact, software architecture, for model based control (Díaz-Rodríguez et al., 2013), relies on minimizing computational task timing commonly constrained by fast sampling periods (Goodwin et al., 1992).Similarly, in other applications, such parameter identification (Mata et al., 2008), and internal redundancy (Parsa et al., 2013) few papers have focused on revisiting and comparing computationally-inexpensive methods in order to obtain the dynamic model of the 3-PRS configuration for cost-effective real-time applications.
To this end, this paper is organized as follows: Sect. 2 presents in general terms how the dynamic model for a parallel manipulator is developed for the three approaches investigated in this paper.Section 3 presents the implementation of the approaches for solving the dynamics problem of a 3-PRS spatial parallel manipulator.Section 4 summarizes and discusses the complexity and the computational load of these three formulations.Finally, the conclusions are drawn.

Development of the dynamic models
The dynamic model of a closed chain mechanical system such as a parallel manipulator can be obtained by virtually cutting or splitting the manipulator at one or more of its joints until the complete dynamic model of a tree-like system with several open chains is obtained.Newton-Euler formulations are then used for solving the dynamics of each serial chain.Finally, constraint equations obtained by means of the Lagrange Multipliers are incorporated to include the necessary forces at the splitting points as to ensure the kinematic chains remain closed.On the other hand, the Lagrangian approach can be applied by using the Lagrange equation with respect to a minimum set of generalized coordinates.Yiu et al. (2001) showed that either application of Lagrange Equation or treelike system analysis are equivalent to one another and lead to the same set of equations when applied to a parallel manipulator.Moreover, similar results were obtained by Murray and Lovell (1989) using the D'Alembert's principle and the principle of virtual work.
Regardless, of the dynamics equation (Newton-Euler, Lagrange or Principle of Virtual Work) used for developing the model, the dynamics of closed chain system can be written as: This equation establishes that the generalized forces corresponding to an open chain system (h) are related to the actuated forces of a closed chain system (τ ) by a linear transformation G T .This linear transformation is based on the Jacobian matrix.Thus, the dynamics equation of a parallel manipulator essentially relies on finding the Jacobian matrix that relate the passive generalized coordinates to the actuated ones.
The following subsections revisit three approaches for solving the dynamics equations of a 3-PRS parallel manipulator.The differences in the approaches are based on two facts: (1) how the parallel manipulator is split into open chains and (2) the type of generalized coordinates used for modelling the manipulator.Each approach leads to establishing different formulations of the Jacobian matrix.Since each Jacobian matrix holds different formulation there are differences in the number of algebraic operations to compute each.The objective of this work is mostly revisiting these approaches to find the cost-effectiveness of each when solving the inverse dynamics through an example.This is of particular interest as inexpensive computer models are essential when adaptive control algorithms are used to control manipulator at fast update rates.For instance, in the fields of parameter identification and actuation redundancy, a fast estimation of the inverse model implies cost-effectiveness for real-time applications.

Dynamics considering the legs and platform as subsystems
This approach is based on the general formulation for modelling parallel manipulators presented by Khalil and Ibrahim (2007), which is based on the following aspects: (1) the manipulator is split open at the spherical joints so that the moving platform is separated from the legs and (2) the local joint coordinates systems q can be used to develop the dynamic equations of each leg h i while the Cartesian coordinates x are used to obtain the dynamic equations of the platform h p .Then, the dynamic equations are combined and projected onto the active joint space as follows: where J p is the Jacobian projecting the task space coordinates (6 in the general case) to the n active joint coordinates, while m is the number of joints for each leg.Likewise, where m p is the mass of the platform, I p denotes the inertial matrix of the moving platform about its centre of mass, a p stands for the acceleration of the end effector, and ω p , ωp respectively denote the angular velocity and angular acceleration of the platform.

Method I
In order to develop the model in actuated joint space one has to project the passive joint variables to the active ones.That is: where indices a and p stand for the active and passive joints, respectively, while G I is a l ×n matrix projecting the dynamics from the passive to the active joints.Here, l represents the number of passive joints while n is the number of active joints on each leg.Equation ( 4) can be written in the form of Eq. ( 1).That is: where I is the identity matrix with dimension equal to the degree of freedom of the manipulator.In Eq. ( 4), matrix G I can be obtained by considering the fact that the distance among spherical joints at the platform is constant due to the rigid body assumption.This distance is calculated based on the norm of the vector obtained by subtracting the position vector identifying the location of the spherical joints.The partial derivatives of each equation with respect to the joints coordinates yields matrix G I .Figure 1 shows how to establish the closed chain equations.
Note that in Eq. ( 4), h p is a 6 × 1 vector and, in the particular case of a 3-PRS parallel manipulator, τ is a 3 × 1 vector.Therefore, matrix J p is not square and cannot be obtained using previous methods for solving the inverse kinematics of this kind of manipulator.For instance, when developing the inverse dynamics model using the Principle of Virtual Work, Li and Xu (2004) found a 3 × 3 square matrix J p .They did so by only considering three of the components of the platform inertial forces; they considered those associated with the desired degrees of freedom of the end effector.In Sect. 3 a method obtaining J p considering the 6 components of h p for the particular case of the 3-PRS manipulator is shown.

Method II
Another approach to formulate the dynamic model is to find the Jacobian matrices according to the approach presented by Khalil and Ibrahim (2007).The method is based on projecting the dynamics equation of the passive joint onto the task www.mech-sci.net/7/9/2016/Mech.Sci., 7, 9-17, 2016 space, and then, project them back to the active joint space so that: where G II is a l × 6 matrix that holds new definition that can be written as follow: where J v i and J −T q i can be obtained respectively from the direct Jacobian J x and the inverse Jacobian J q of the manipulator.
Equation ( 6) can be written in the form of Eq. ( 1).That is:

Dynamics considering the manipulator as open kinematic chains
A parallel manipulator can be split open into m − 1 joints yielding m open chain systems.In this approach the platform is attached to one of the legs, see Fig. 2. Algorithms for obtaining the dynamics model of serial manipulators may now be applied to obtain the dynamics of each leg.
The cut joints introduce constraint forces, which can be included into the model by means of the Lagrange multipliers: where A is the Jacobian that can be found by analysing the constraint equations.Mata et al. (2008) presented a method for the Jacobian matrix by considering the linear velocity at the split joints.The velocities can be computed through the Jacobian analysis of each leg.Then, the velocity obtained at the split joint following each leg are the same.In this approach, recursive modelling of the velocity analysis from conventional serial manipulator methods can be applied for each leg.
Once the Jacobian matrix is found, the Lagrange multipliers in Eq. ( 9) are eliminated by multiplying matrix C, so that, CA = 0.One way to find C is by obtaining the natural orthogonal complement of C. On the other hand, the matrix can be found by separating matrix A = 0 into the passive and the active joints.That is: where G III has dimensions l × n and can be computed as follows: where subscripts p and a respectively refer to the passive and active variable terms in Jacobian A. Equation ( 10) can be written in the form of Eq. ( 1).That is: As Eqs. ( 5), ( 8) and ( 12) show, the three methods establish a linear relation between the cut-open model (h) to the original system.The different lies in how matrix G T and vector h are solved in each model.This is illustrated in the next section for the particular case of a general 3-PRS parallel manipulator.

Inverse dynamics of the 3-PRS PM
A schematic representation of the 3-PRS PM is shown in Fig. 3 where the 7 moving rigid bodies are shown.There, links 0, 1, and 2 can be seen as 2-DOF serial manipulator with PR joints, also it applies to links 0, 4, and 5 as well as links 0, 6, and 7.The platform is indicated by the number 3. The manipulator is a lower mobility (i.e., less that 6-DOF) spacial PM with 3-DOF.This manipulator holds the characteristic of zero torsion at its platform because the three spherical joints move in vertical planes intersecting at a common line (Liu and Bonev, 2008).In addition, the topology of its legs provides 2-DOF of angular rotation (2R) in two axes (A/B axis, rolling and pitching) and 1-DOF translation (1T) motion (heave) at the end effector.
As Fig. 3 shows, reference frame {O} defines the global coordinate system (x, y, z) while {p} defines the local coordinate system attached to the moving platform.In the typical symmetrical configuration, the centre of the spherical joints at the platform form an equilateral triangle circumscribed in a circle with centre p and radius r p .The line of action of the prismatic joints intersects the base Oxy plane at A i also forming an equilateral triangle.The distance from O to A i is the base platform radius r b .
In order to take advantage of the dynamic equation algorithms developed for serial manipulators, when modelling the manipulator with leg and platform as subsystem, the joint coordinates q can be used to develop the dynamics equations of each leg (Mata et al., 2002).To this end, a local coordinate system {O i,0 } is defined at the bottom of leg i where the leg meets the based plane.Table 1 lists the Denavit-Hartenberg (D-H) parameters, according to Craig's notation (Craig, 2005), of each leg i from {O i,0 } up to the location of the axis {O i,2 } at the revolute joint.
In addition, the roll-pitch-yaw (α, β and φ) Euler angles represent the orientation of the moving frame {p} with respect to the global coordinate system {O}.The rotation matrix O R p is defined as: where c * = cos ( * ) and s * = sin ( * ).
Table 1.D-H parameters for a 3-PRS PM when modelling the manipulator with leg and platform as subsystem.
The task space (x) and joint space (q i ) coordinates are given by: ) where s i represents the displacement along the axis of the prismatic joint i, and θ i the angle of the link 2 and the axis of the prismatic joint i in the plane of movement of the leg i.
The components of b i with respect to the local coordinate system {p} are given by: while the components of a i with respect to the global coordinate system {O} are given by: , and The position problem for the considered PM is not included in this paper since it can be found in Carretero et al. (2000b), Tsai et al. (2002), Li and Xu (2004), Mata et al. (2008).The following subsections focus on the computation of the Jacobian matrix for the aforementioned approaches.

Jacobian matrix for Model I
Matrices J p and G for computing Eq. ( 4) are found following the approach presented in Li and Xu (2004).The vector loop equation of the ith leg can be written as: where u 1 i and u 2 i are unit vectors, l i is the distance between the rotational joint and the spherical joint, b i is the position vector from the origin O to the centre of the spherical joint.Vector a i is the position vector between O and local frame of www.mech-sci.net/7/9/2016/Mech.Sci., 7, 9-17, 2016 each leg.Differentiating Eq. ( 20) with respect to time and after some algebraic manipulation, the following equation can be obtained: where: and Due to the constraints imposed by the fact that each legs moves on a plane, the following set of equations holds: From these equations, a 6 × 3 Jacobian matrix mapping the dependent task space coordinates to the independent ones can be found such that: It is important to note that ẋ = ẋp ẏp żp ω p x ω y p ω z p T .In order to apply Eq. ( 25), one has to find the angular velocity of the platform through the rate of change of the generalised coordinates α β φ T .That is: After some algebraic manipulation the following equation can be obtained: where J c = J −1 q J x J r .The 3 × 3 matrix G I is found by considering the fact that the distance among spherical joints at the platform is constant due to the rigid body assumption.That is: with i = 1, 2, 3 and when i = 3, i + 1 = 1.
Thus, by obtaining the partial derivatives of the above set of equation matrix G I can be written as:

Jacobian matrix for Model II
Another approach to compute matrices J p and G is to consider that the spherical joints in each leg is constrained to move on a plane normal to the revolute joint.The motion of each leg at point B i can be found in terms of the joint coordinates.Moreover, it can also be expressed with respect to O i,2 .Due to the constraints provided by the P-R pair, the third row of the Jacobian matrix consist of zero entries.That is: where The linear velocity of the end effector can be related to the linear velocity of points B i as follows: stands for the skew symmetric matrix substituting the cross product b i ×, and The first two rows of Eqs. ( 31) and (32) relate the task space to the joint spaces coordinates.That is: where Table 2. D-H Parameters for the first leg of the 3-PRS PM when modelling as a three open chains. and, In the above equations x = 1 0 0 T , y = 0 1 0 T , and 0 is a 2 × 2 zero matrix.
From the above equation, matrix G II can be obtained as follows: where A (:, 2) denotes the 2nd column of matrix A.
The Jacobian matrix, relating the task space coordinates J r , is found by considering the third column of velocity equation following each leg and the platform.The Jacobian matrix is obtained by following method 2 which is graphically represented in Fig. 4.

Jacobian matrix for Model III
The inverse dynamic is computed as a function of three open chains which are obtained after disassembling two of the three spherical joints.The platform is attached to one of the legs, and the spherical joint is modelled as three intersecting revolute joints with the three axes mutually perpendicular to one another.Therefore, the chain with the end effector platform is modelled by using the set of D-H parameters presented in Table 2.The remaining legs have only sets of two variables which are the same as those presented in Table 1.
One of the advantages of considering the manipulator as tree-like serial chains is that the Jacobian for each leg can be computed by using well-known recursive modelling from serial manipulator.In this respect, the velocities at the cut joints can be computed through recursive formulation (Angeles, 2002).That is: where i = 1, 2, 3 and when i = 3, i+1 = 1, A i is the Jacobian matrix for the ith leg, and V B i,i+1 is the velocity at the cut joint connecting leg i and i + 1.From Eq. ( 39) the following equation can be obtained: This equation provides a set of three linear systems relating the joint coordinates following each loop.If the set of linear equations is appended together the following equation holds: The relationship between the active and passive coordinates can be obtained from Eq. ( 41) as follows:

Results and discussion
In order to solve the inverse dynamics of the 3-PRS parallel manipulator, each term of Eqs. ( 4), (6), and (10) are found in closed form by using a Computer Algebra Symbolic (CAS) program, such as Maple.One of the advantages of developing the model in a CAS program relies on the fact that the mathematical operations can be performed symbolically and simplified.In the present case, built-in functions such as simplify and combine (with option=trig) of Maple programming environment were used to reduce the number of operations for solving each model.
A second advantage of obtaining the model in closed form is that the code can be written automatically for Matlab by using the code generation capabilities of the software.The Matlab procedure with optimize=tryhard option of the package CodeGeneration were used in this case to develop Matlab code.The number of algebraic operations (i.e., additions/subtractions and divisions/multiplications) necessary for solving the dynamic problem was obtained through the cost function of the codegen package.Without any loss of generality, the number of operations for matrix inversion (i.e., when obtaining J p ) were computed by considering the number of operations conventional LU decomposition takes to solve a linear system: about +, −n 3 /3− n 2 /2 + 5n/6 and × n 3 /3 + n 2 − n/3 (Chapra and Canale, 2006).The objective of this paper is to evaluate the complexity and the computational load of these three formulations which is presented in Table 3.
www.mech-sci.net/7/9/2016/Mech.Sci., 7, 9-17, 2016  As Table 3 shows, the models in Eqs. ( 4) and ( 6), which are based on splitting the platform from the legs, hold fewer number of operation than those of the model considering the platform attached to one of the legs.This fact is due to the inversion of a 6 × 6 matrix when finding the matrix G III .On the other hand, by having the platform attached to one leg, the projection of the platform generalised forces onto the actuated joints is cumbersome.The result shows that either using Eqs.(4) or (6) a reduction of about 30 % in the number of multiplication and about 25 % in additions is obtained when comparing to the model given in Eq. 10).These results indicate that considering the platform and legs as subsystem can improve speed in the calculation of dynamics for applications such as model-based control, see for example Díaz-Rodríguez et al. (2013).

Conclusions
Three strategies for dynamic modelling of parallel manipulators were revisited and applied for developing the inverse dynamic model of a general 3-PRS parallel manipulator.The first method exploited parallel manipulators particularities such that the platform and the links are considered as subsystems, namely: the legs and the platform, thus, joint coordinates were used for modelling the legs dynamics and task space coordinates for the platform.The core concept of the second method is based on solving the dynamic model through this approach reduces to how the Jacobian matrix relating the task and the joint space coordinates is computed and it was explained in this paper.A third method was presented which is similar to the second method but uses a different formulation of the Jacobian matrix.In the third approach the dynamic model was obtained by describing the manipulator as three open kinematic chains.The vector-loop closure constraints introduced the relationship between the dynamics of the open kinematic chains and the original closed chains.A Computer Algebraic Software allowed to find each term of the dynamic in symbolic form and to compute the computational burden of each model.
The results showed that the approaches based on splitting the manipulator in two sub-systems (platform and legs) require about 30 % in the number of multiplication and about 25 % in additions are obtained when comparing to the model given by Eq. ( 10).These results indicate that in problems when the model is needed to be computed on-line at high rates of speed, method 1 and 2 can be useful.This work has provided numerical guidelines for implementing computationally efficient models for use in numerically intensive optimal mechanical synthesis problems or in resourceconstraint embedded computers, particularly for control and model identification.The software support presented for solving the inverse dynamic problem efficiently also provides some insight on some of the advantages and/or disadvantages on these revisited methods.

Figure 1 .
Figure 1.Closed chain equation for finding G I .

Figure 2 .
Figure 2. Several open chains obtained after cutting open the parallel manipulator.

Figure 4 .
Figure 4. Formulation of the Jacobian matrix using Method II.

Table 3 .
Computational effort for solving the inverse dynamic problem.