On the Robust Dynamics Identification of Parallel Manipulators: Methodology and Experiments

1. Introduction

The proposed chapter presents a self-contained approach for the dynamics identification of parallel manipulators. Major feature is the consequent consideration of structural properties of such machines in order to provide an experimentally adequate identification method. Thereby, we aim to achieve accurate model parameterization for control, simulation or analysis purposes. Despite the big progress made on identification of serial manipulators, it is interesting to state the missing of systematic identification methodologies for closed-loop and parallel kinematic manipulators (PKM’s). This is due to many factors that are discussed and treated systematically in this chapter.

First, the issue of modelling the dynamics of PKM in a linear form with respect to the parameters to be identified is addressed. As it is already established in the field of classic serial robotics, such step is necessary to ensure model identifiability and to apply computationally efficient linear estimation (Swevers et al., 1997, Khalil & Dombre, 2002, Abdellatif & Heimann, 2007). The case of parallel manipulators is more complicated, since a multitude of coupled and closed kinematic chains has to be considered (Khalil & Guegan, 2004, Abdellatif et al., 2005a). Beside the rigid-body dynamics, friction plays a central role in modelling, since its accurate compensation yields important improvement of control accuracy. If friction in the passive joints is regarded, the dimension of the parameter vector grows and affects the estimation in a negative way. To cope with such problem, a method for the reduction of friction parameter number is proposed, which is based on the identifiability analysis for a given manipulator structure and by considering the actual measurement noise. The calculation procedure of a dynamics model in minimal parametrized form is given in section 2.

Another important issue of PKM’s is the appropriate design of the identification experiment, in order to obtain reliable estimation results is. Two aspects are here crucial: The choice or the definition of the experiment framework at the one hand and its related experiment optimization at the other hand. Regarding the first aspect, the harmonic excitation approach proposed a couple of years ago for serial manipulators is chosen (Swevers et al., 1997). The method provides bounded motion that can be fitted in the usually highly restricted and small workspace of parallel robots. Thus, we propose an appropriate adaptation for PKM’s. The experiment optimization is carried out within a statistical frame in order to accounts for the cross correlation of measurement noise and the motion dependency of the coupled actuators (Abdellatif et al., 2005b). The Experiment design is discussed in section 3.

The typically non measurable information of the end-effector postures, velocities and accelerations are necessary to calculate the dynamics model and therefore to obtain the regression equation. Since in general only actuator measurements are available, there is a need for an adequate estimation of the executed end-effector motion during the identification experiment. However, the numerical computation of the direct and the differential kinematics yields a spectral distortion and noise amplification in the calculated data. Therefore, an appropriate and simple frequency-domain data processing method is introduced in section 4. An accurate and noise-poor regression model is then provided, which is crucial for bias-free estimation of the model parameters. Additionally, we provide useful relationships to evaluate the resulting parameter uncertainties. Here, uncertainties of single parameters as well as the uncertainties of entire parameter sets are discussed and validated.

Figure 1.

Case Study: Hexapod PaLiDA; left: Presentation in the Hannover industrial Fair, Right: CAD-Model

Finally and in section 5, an important part of the chapter presents the experimental substantiation of the theoretical methods. The effectiveness of our approach is demonstrated on a six degrees-of-freedom (dof’s) directly actuated parallel manipulator PaLiDA (see Fig. 1). We address the important issue of exploiting the identification results for model-based control. The impact of accurately identified models on the improvement of control accuracy is illustrated by numerous of experimental investigations.

2. Parameterlinear formulation of the dynamics model

The objective of this section is to derive the inverse dynamics model in a linear form with respect to a set of the parameters to be identified. Such formulation allows for using linear techniques to provide the estimation of model parameters from measurement data. This kind of approach is well established for serial robots (Khalil & Dombre, 2002, Abdellatif & Heimann, 2007). Thereby, the model accounts for the rigid-body as well as for friction dynamics. We consider the case of 6-dof’s parallel manipulator, that is constituted of a moving platform (end-effector platform) attached with six serial and non-redundant actuated kinematic chains to the base platform. Fig. 2 shows a general sketch of such robotic manipulator. Let $$x,v,a$$ be the 6-dimensional vectors denoting the posture, the velocity, and the acceleration of the end-effector, respectively. The posture vector is composed of the cartesian coordinates of the end-effector platform $${}_{(0)}r{}_{\text{E}}^0={\left[\begin{array}{ccc}x& y& z\end{array}\right]}^{\text{T}}$$ and the tilting angles $$\left(\phi ,\theta ,\psi \right)$$ according to the Cardan or the Euler formalism. The velocity vector is defined as $$v={\left[\begin{array}{cc}{}_{(0)}v{}_{\text{E}}^{\text{T}}& {}_{(0)}\omega {}_{\text{E}}^{\text{T}}\end{array}\right]}^{\text{T}}$$ that includes the translational and angular velocities with reference to a cartesian frame. It is known, that $$v\ne \dot{x}$$ holds for systems with two or more rotational dof's (Merlet, 2006). The 6-dimensional vector of actuated joints is denoted by $$q_{\text{a}}$$ . The passive joint variables are grouped in $$q_{\text{p}}$$ . The vector $$q$$ contains all joint variables.

Figure 2.

Scheme of a general parallel manipulator,

The major difference between serial and parallel manipulators is the definition of configuration variables or the configuration space. For classic serial manipulators, the actuation variables $$q_{\text{a}}$$ are sufficient to determine exactly the system's configuration. This is not the case for PKM’s, because the solution of the direct kinematics is ambiguous (Khalil & Dombre, 2002,Merlet, 2006). It is established that the motion of the end-effector given by $$x,v$$ and a is used to derive the dynamics of high-mobility parallel robots (Tsai, 2000, Harib & Srinivasan, 2003, Khalil & Guegan, 2004, Abdellatif et al., 2005a). The solution of inverse kinematics is supposed to be already achieved. It means that for a given dynamic motion of the end-effector, all necessary kinematic quantities are available. The latter include: velocities and acceleration of any body $$i$$ with respect to a body-fixed frame

The body-fixed frames can be defined according to the modified Denavit-Hartenberg (MDH) notation (Khalil & Dombre, 2002).

$$v_i,a_i$$ or those of the center of mass $$v_{{\text{S}}_i},a_{{\text{S}}_i}$$ ; the angular velocities and angular accelerations $${\omega }_i,{\dot{\omega }}_i$$ ; the body Jacobians $$J_{{\text{T}}_i}=\partial v_i/\partial v$$ and $$J_{{\text{R}}_i}=\partial {\omega }_i/\partial v$$ and the inverse Jacobian of the manipulator $$J^{-1}=\partial q_{\text{a}}/\partial v$$ (see (Abdellatif et al., 2005a) and references therein for more details).

The aimed dynamics model consists of the following equation:

$$Q_{\text{a}}=A\left(x,v,a\right)p\iff Q_{\text{a}}=A_{\text{a}\text{,rb}}\left(x,v,a\right)p_{\text{rb}}+A_{\text{a}\text{,f}}\left(x,v,a\right)p_{\text{f}}$$ E1

with $$Q_{\text{a}}$$ being the actuator forces and where the indexes rb and f refer to the rigid-body and friction terms, respectively.

2.1. Parameterlinear formulation of the rigid-body dynamics

Generally, it is recommended to use the Jourdain’s principle of virtual power to derive the dynamics in an efficient manner. In analogy to the virtual work, a balance of virtual power can be addressed:

$$\delta v^{\text{T}}\tau =\delta {\dot{q}}_{\text{a}}^{\text{T}}{\text{Q}}_{\text{a}}\iff \tau {\left(\frac{\partial {\dot{q}}_{\text{a}}}{\partial v}\right)}^{\text{T}}{Q}_{\text{a}}$$ E2

where $$\tau$$ is the vector of the generalized forces, defined with respect to the end-effector generalized velocities $$v$$ . Equation (2) means that the virtual power resulting in the space of generalized velocities is equal to the actuation power. The power balance can be applied for rigid-body forces:

$$Q_{\text{a}\text{,rb}}=\tau {\left(\frac{\partial {\dot{q}}_{\text{a}}}{\partial v}\right)}^{\text{T}}{\tau }_{\text{rb}}=J^{\text{T}}{\tau }_{\text{rb}}$$ E3

The generalized rigid-body forces for a manipulator with $$N$$ bodies are obtained by

$${\tau }_{\text{rb}}={\displaystyle \sum _{i=1}^N\left[J_{{\text{T}}_i}^{\text{T}}\left(m_i{}_{(i)}\dot{v}{}_i+{}_{(i)}\dot{\tilde{\omega }}{}_i{s}_i+{}_{(i)}\tilde{\omega }{}_i{}_{(i)}\tilde{\omega }{}_i{s}_i\right)+J_{{\text{R}}_i}^{\text{T}}\left({}_{(i)}I{}_i^{(i)}{}_{(i)}\dot{\omega }{}_i+{}_{(i)}\tilde{\omega }{}_i\left({}_{(i)}I{}_i^{(i)}{}_{(i)}\omega {}_i\right)+{\tilde{s}}_i{}_{(i)}\dot{v}{}_i\right)\right]}$$ E4

with the dynamic parameters of each body $$i$$ , its mass $$m_i$$ , its statical first moment $$s_i\widehat{=}{\left[\begin{array}{ccc}{s}_{i_x}& s_{i_y}& s_{i_z}\end{array}\right]}^{\text{T}}=m_i{}_{(i)}r{}_{S_i}^i$$ $${}_{(i)}r{}_{s_i}^i$$ : Vector from coordinate frame to centre of mass) and its inertia tensor about the corresponding body-fixed coordinate frame $${}_{(i)}I{}_i^{(i)}$$ . New operators $${\left(\right)}^{\ast }$$ and $${\left(\right)}^{\diamond }$$ are defined:

$${\omega }_i^{\ast }{I}_i^{\diamond }\widehat{=}{}_{(i)}I{}_i^{(i)}{}_{(i)}\omega {}_i$$ E5

with

$${\omega }_i^{\ast }=\left[\begin{array}{cccccc}{\omega }_{i_x}& {\omega }_{i_y}& {\omega }_{i_z}& 0& 0& 0\\ 0& {\omega }_{i_x}& 0& {\omega }_{i_y}& {\omega }_{i_z}& 0\\ 0& 0& {\omega }_{i_x}& 0& {\omega }_{i_y}& {\omega }_{i_z}\end{array}\right]$$ and

$$I_i^{\diamond }={\left[\begin{array}{cccccc}{I}_{i_{xx}}& I_{i_{xy}}& I_{i_{xz}}& I_{i_{yy}}& I_{i_{yz}}& I_{i_{zz}}\end{array}\right]}^{\text{T}}$$ E6

which helps the simplification of the generalized rigid-body forces:

$${\tau }_{\text{rb}}={\displaystyle \sum _{i=1}^N\underset{H_i}{\underbrace{\left[\begin{array}{cc}{J}_{{\text{T}}_i}^{\text{T}}& J_{{\text{R}}_i}^{\text{T}}\end{array}\right]{\Omega }_i}}}\underset{p_i}{\underbrace{\left[\begin{array}{c}{I}_i^{\diamond }\\ s_i\\ m_i\end{array}\right]}}=\left[\begin{array}{cccc}{H}_1& H_2& \dots & H_N\end{array}\right]{\left[\begin{array}{cccc}{p}_1^{\text{T}}& p_2^{\text{T}}& \dots & p_N^{\text{T}}\end{array}\right]}^{\text{T}}$$ E7

with

$${\Omega }_i=\left[\begin{array}{ccc}0& {}_{(i)}\dot{\tilde{\omega }}+{}_{(i)}\tilde{\omega }{}_i{}_{(i)}\tilde{\omega }{}_i& {}_{(i)}\dot{v}{}_i\\ {}_{(i)}\dot{\omega }{}_i^{\ast }+{}_{(i)}\tilde{\omega }{}_i{}_{(i)}\omega {}_i^{\ast }& -{}_{(i)}\dot{v}{}_i& 0\end{array}\right]$$ E8

Equation (7) is already linear with respect to the parameter vector $$p_{\text{rb}}^{\ast }=\left[\begin{array}{cccc}{p}_1^{\text{T}}& p_2^{\text{T}}& \dots & p_N^{\text{T}}\end{array}\right]{}^{\text{T}}$$ . The dimension of the latter has now to be reduced for an efficient calculation and to assure the identifiability of the system. The proposed algorithm in the following is based on former works for serial and parallel manipulators (Grotjahn & Heimann, 2000, Grotjahn et al., 2002). The matrices $$H_i$$ in eq. (7-8) can be grouped in single serial kinematic chains such that a recursive calculation:

$$H_i=H_{i-1}{L}_i+K_i$$ E9

can be achieved. The matrices $$L_i$$ and $$K_i$$ are given in (Grotjahn et al., 2002). The first step considers in eliminating all parameters $$p_{\text{rb}\text{,}j}^{\ast }$$ that correspond to a zero column $$h_j$$ of $$H$$ , since they do not contribute to the dynamics. The remaining parameters are then regrouped to eliminate all linear dependencies by investigating $$H$$ . If the contribution of a parameter $$p_{\text{rb}\text{,}j}^{\ast }$$ depends linearly on the contributions of some other parameters $$p_{\text{rb}\text{,1}j}^{\ast },\dots ,p_{\text{rb}\text{,}kj}^{\ast }$$ , the following equation holds:

$$h_j={\displaystyle \sum _{l=1}^k{a}_{lj}}{h}_{lj}$$ E10

Then $$p_{\text{rb}\text{,}j}^{\ast }$$ can be set to zero and the regrouped parameters $$p_{\text{rb}\text{,}lj,\text{new}}^{\ast }$$ can be obtained by

$$p_{\text{rb}\text{,}lj,\text{new}}^{\ast }=p_{\text{rb}\text{,}lj}^{\ast }+a_{lj}{p}_{\text{rb}\text{,}j}^{\ast }$$ E11

The recursive relationship given in (9) can be used for parameter reduction. If one column or a linear combination of columns of $$L_i$$ is constant with respect to the joint variable and the corresponding columns of $$K_i$$ are zero columns, the parameters can be regrouped. This leads to the rules which are formulated in (Khalil & Dombre, 2002) and in (Grotjahn & Heimann, 2000).

The rules can be directly applied to the struts or legs of the manipulator, since they are considered as serial kinematic chains. For revolute joints the 9^th, the 10^th and the sum of the 1^st and 4^th columns of $$L_i$$ and $$K_i$$ comply with the mentioned conditions. Thus, the corresponding parameters $$I_{i_{\text{yy}}}$$ $$s_{i_{\text{z}}}$$ and $$m_i$$ can be grouped with the parameters of the antecedent joint $$i-1$$ . For prismatic joints however, the moments of inertia can be added to the carrying antecedent joint, because the orientation between both links remain constant.

The end-effector platform closes the kinematic loops and further parameter reduction is possible. The velocities of the platform joint points $$B_j$$ (see Fig. 2) and those of the terminal fixed-body frames of the respective legs are the same, yielding dependencies of the respective energy-functions. The masses of terminal bodies can be grouped to the inertial parameter of the platform according to Steiner's laws.

After Applying every possible parameter reduction the generalized rigid-body forces are obtained from (7) with respect to a minimal set of parameters $${\tau }_{\text{rb}}=A_{\text{rb}}{p}_{\text{rb}}$$ . In combination with (3) the desired form for the rigid-body part of the actuation forces is obtained as

$$Q_{\text{a}\text{,rb}}=J^{\text{T}}{A}_{\text{rb}}{p}_{\text{rb}}=A_{\text{a}\text{,rb}}{p}_{\text{a}\text{,rb}}$$ E12

2.2. Parameterlinear formulation of the friction forces

In analogy to the rigid-body dynamics, the Jourdain’s principle can be applied for friction forces. By defining an arbitrary steady-state model at joint-level $$Q_{\text{f}}=f_{\text{r}}(\dot{q})$$ , a new power balance can be derived:

$$Q_{\text{a}\text{,f}}={\left(\frac{\partial \dot{q}}{\partial {\dot{q}}_{\text{a}}}\right)}^{\text{T}}{Q}_{\text{f}}=J^{\text{T}}{\left(\frac{\partial \dot{q}}{\partial v}\right)}^{\text{T}}{Q}_{\text{f}}$$ E13

Equation (13) means that the friction dissipation power in all joints (passive and active) has to be overcome by an equivalent counteracting actuation power. We notice that the case of classic open-chain robots correspond to the special case, when the joint-Jacobian $$\partial \dot{q}/\partial {\dot{q}}_{\text{a}}$$ is equal to the identity matrix. In the more general case of parallel mechanisms, friction in passive joints should not be neglected (Abdellatif & Heimann, 2006).

For identification purpose, friction in robotics is commonly modelled as superposition of Coulomb (or dry) friction and viscous damping depending on joint velocities $${\dot{q}}_i$$ (Abdellatif et al., 2007, Swevers et al., 1997):

$$Q_{{\text{f}}_i}=f_{\text{r}}({\dot{q}}_i)=\left[\text{sign}({\dot{q}}_i){\dot{q}}_i\right]\left[\begin{array}{cc}{r}_{1_i}& r_{2_i}\end{array}\right]{}^{\text{T}}$$ E14

Regrouping friction forces in all $$n$$ joints yields to

$$Q_{\text{f}}=\underset{D_{\text{f}}}{\underbrace{\left[\begin{array}{cc}{D}_1(\dot{q})& D_2(\dot{q})\end{array}\right]}}\underset{p_{\text{f}}}{\underbrace{\left[\begin{array}{cc}{r}_1^{\text{T}}& r_2^{\text{T}}\end{array}\right]{}^{\text{T}}}}$$ E15

with

$$r_k^{\text{T}}=\left[r_{k_1},\dots ,r_{k_{\text{n}}}\right]$$ E16

$$D_1(\dot{q})=\text{diag}\left\{\text{sign}({\dot{q}}_1),\dots ,\text{sign}({\dot{q}}_n)\right\}$$ E17

$$D_2(\dot{q})=\text{diag}\left\{{\dot{q}}_1,\dots ,{\dot{q}}_n\right\}$$ E18

Considering (13) and (15) the linear form of the resulting friction forces in the actuation space is obtained

$$Q_{\text{a}\text{,f}}=\left[J^{\text{T}}{\left(\frac{\partial \dot{q}}{\partial v}\right)}^{\text{T}}{D}_{\text{f}}\right]p_{\text{f}}=A_{\text{a}\text{,f}}(x,v)p_{\text{f}}$$ E19

Unlike the rigid-body dynamics, there is no uniform or standard approach for the reduction of the parameter vector dimension. In a former publication, we proposed a method that is highly adequate for identification purposes. Thereby, the expected correlation of the friction parameter estimates is analyzed for a given and statistically known measurement disturbance. Parameters whose effects are beneath the disturbance level are eliminated. Parameters with high correlation are replaced by a common parameter. The interested reader is here referred to (Abdellatif et al., 2005c) and (Abdellatif et al., 2007) for a deep insight.

3. Identification experiment design for parallel manipulators

Almost all identification methods in robotics are based on the parameterlinear form that is given by (1) in combination with (12) and (19) (Swevers et al., 1997, Khalil & Dombre, 2002, Abdellatif & Heimann, 2007). Given experimentally collected and noise corrupted $$N$$ measurement sets, the estimation problem can be formulated according to (1) as

$$\underset{\Gamma }{\underbrace{\left[\begin{array}{c}{Q}_{{\text{a}}_{\text{1}}}\\ ⋮\\ Q_{{\text{a}}_N}\end{array}\right]}}=\underset{\Psi }{\underbrace{\left[\begin{array}{c}A\left(x_1,v_1,a_1\right)\\ ⋮\\ A\left(x_N,v_N,a_N\right)\end{array}\right]}}p+\underset{\eta }{\underbrace{\left[\begin{array}{c}{e}_1\\ ⋮\\ e_N\end{array}\right]}}$$ E20

with the measurement vector $$\Gamma$$ , the information or regression matrix $$\Psi$$ and the error vector $$\eta$$ that accounts for disturbances. The most classic and simple solution of the overdetermined equation system (20) can be achieved by the Least-Squares (LS) approach. However, such method assumes that the disurbances of the different actuators are not cross correlated. The assumption does not hold for high-coupled systems like the case of parallel manipulators (Abdellatif et al., 2005b). It is recommended to use the Gauss-Markov (GM) approach that presents a more general case

$${\widehat{p}}_{\text{GM}}={\left({\Psi }^{\text{T}}{\Sigma }_{\eta \eta }^{-1}\Psi \right)}^{-1}{\Psi }^{\text{T}}{\Sigma }_{\eta \eta }^{-1}\Gamma$$ E21

The crosscoupling is regarded by the full covariance matrix $${\Sigma }_{\eta \eta }=\text{E}\left(\eta {\eta }^T\right)$$ of the measurement disturbances $$\eta$$ . Neglecting this fact by applying the simple LS-method will lead to biased estimates (Abdellatif et al., 2005b).

3.1. Design of the excitation trajectory

An important step in identification is the choice of the measurement data to be collected. A classic choice consists in the so-called excitation trajectory, which ensures that the effects of all considered parameters are contained in the measurement data. A challenging issue with parallel manipulators is their restricted and highly constrained workspace. Such property reduces the possibility of highly dynamic and variable motion that is necessary for the excitation of all parameters to be identified. The appropriate choice should be a trajectory that is naturally bounded to fit into a small workspace. An attractive approach is the harmonic excitation approach originally proposed by Swevers et al. (Swevers et al., 1997) and adapted in the following for the case of parallel manipulators.

For each posture coordinate corresponding to the $$i^{th}$$ element of $$x$$ a respective trajectory with $$n_{\text{h}}$$ harmonics is defined as

$$x_i=x_0^i+{\displaystyle \sum _{k=1}^{n_{\text{h}}}\left(\frac{{\mu }_{\text{k}}^{\text{i}}}{k{\omega }_f}\sin \left(k{\omega }_{f}t\right)-{\displaystyle \sum _{k=1}^{n_{\text{h}}}\frac{v_{\text{k}}^{\text{i}}}{k{\omega }_f}}\cos \left(k{\omega }_{f}t\right)\right)}$$ E22

providing a proper trajectory parameter vector

$${\Xi }_i={\left[x_0^i,{\mu }_1^i,\dots ,{\mu }_{n_{\text{h}}}^i,v_1^i,\dots ,v_{n_{\text{h}}}^i\right]}^{\text{T}}$$ E23

with $${\omega }_{\text{f}}$$ being the fundamental frequency. The difference to the implementation for serial robots is that the excitation trajectory is now defined with respect to $$x$$ (and therefore $$v$$ and $$a$$ ) rather than to the actuator coordinates $$q_{\text{a}}$$ . Such modification is necessary, since the dynamics can determined only in the configuration space defined by $$x$$ . With the proposed modification, a direct relationship between the dynamics to be excited and the trajectory is available. If the excitation trajectory is defined with respect to the actuated coordinates $$q_{\text{a}}$$ , the closure constraints of the parallel manipulator and the numerical calculation of the direct kinematics have to be performed while the optimization and design of the trajectory. First ensures a feasible trajectory and second provides the resulting dynamics in form of the regression model. Both operations increase the solution cost and introduce additional numerical errors.

3.2. Optimization of the excitation trajectory

The next step consists in determining the values of all trajectory parameters

$$\Xi ={\left[{\Xi }_1^{\text{T}},\dots ,{\Xi }_6^{\text{T}},{\omega }_{\text{f}}\right]}^{\text{T}}$$ E24

to provide a best possible excitation of the dynamics parameters. Such procedure is called optimal input experiment design. The design is performed by using constrained nonlinear optimization (Swevers et al., 1997, Gevers, 2005). The required constraints are expressed with respect to the actuation variables

$$\begin{array}{l}{q}_{\text{a}}^{\min }\le q_{\text{a}}(t,\Xi )\le q_{\text{a}}^{\max }\\ {\dot{q}}_{\text{a}}^{\min }\le {\dot{q}}_{\text{a}}(t,\Xi )\le {\dot{q}}_{\text{a}}^{\max },\forall \Xi andt\in \left[\mathrm{0,}{T}_{\text{f}}\right]\\ {\ddot{q}}_{\text{a}}^{\min }\le {\ddot{q}}_{\text{a}}(t,\Xi )\le {\ddot{q}}_{\text{a}}^{\max }\end{array}$$ E25

to account for actuator limitation and therefore indirectly for workspace constraints and dynamics capabilities of the manipulator. The inverse kinematics has to be performed while the optimization, which does not introduce any significant computational cost due to its simplicity (Khalil & Guegan, 2004, Abdellatif et al., 2005a, Abdellatif & Heimann, 2007). Of course, it is possible to express the constraints ad-hoc with respect to $$x$$ $$v$$ and $$a$$ . It depends on the considered manipulator, whether such approach is preferable or not, since it results in different constraints than (25), which can accelerate the convergence of the optimization process. The optimization or the experiment design criterion should contribute to the reduction of parameter uncertainty (Gevers, 2005). To account appropriately for disturbances in the information matrix $$\Psi$$ it is recommended to opt for the D-optimal design

$$\Xi =\text{arg}\left\{\underset{\Xi }{\text{ min}}\left(-\text{ln det}\left({\Psi }^{\text{T}}{\Sigma }_{\eta \eta }^{-1}\Psi \right)\right)\right\}$$ E26

that aims increasing the volume of the asymptotic confidence ellipsoid for the parameter estimates, which is equivalent to the determinant of the inverse of the asymptotic parameter covariance matrix $$P={\left({\Psi }^{\text{T}}{\Sigma }_{\eta \eta }^{-1}\Psi \right)}^{-1}$$ or the Cramér-Rao bound (Gevers, 2005). Due to the complexity of the nonlinear dynamics contained in the regressor $$\Psi$$ the optimization is mostly a non-convex one and the obtained results will not correspond to the global minimum. This is however not critical since for experimental identification just a sufficiently good excitation trajectory is needed.

4. Identification procedure: data processing, implementation and parameter uncertainties

At this stage, the dynamics of the manipulator is available in linear form (section 2). Additionally, the appropriate choice of an excitation experiment is proposed (section 3.1) with a recommended method for its optimal design (section 3.2). Therefore, the experiment can be executed and the data can be collected to achieve an estimation according to (21). Here, the next challenge for parallel manipulators is evident. The measurements are provided in the actuation space in form of actuation forces and actuator positions, whereas the information matrix $$\Psi$$ is built up by using $$x$$ $$v$$ and $$a$$ that are not directly measured. Thus, a reconstruction of these variables from the corrupted measurement of $$q_{\text{a}}$$ is necessary.

4.1. Data processing

The first step consists in calculating the direct kinematics to provide a first estimate of the posture $$\widehat{x}$$ . The terminal condition of the numerical calculation has to be set less than the resolution of the used sensors (Merlet, 2006). The obtained estimate is of course noisy and has to be filtered. Filtering the measurement in the time-domain (i.e. by using classic low-pass filters) may cause lost of information, since ideal and exact filtering is not possible. More critical is the calculation of $$v$$ and $$a$$ . Numerical differentiation of the posture data is not convenient. Additionally to the measurement noise, possible oscillations of the direct kinematic solution introduces disturbances, such that the resulting data may be not useful at all (Abdellatif et al., 2004).

By taking advantage of the periodic and harmonic nature of the excitation trajectory, exact filtering in the frequency-domain can be achieved. First, it is recommended to calculate the DFT-transform of each component $$i$$ of the pre-computed posture $${\widehat{x}}_i\to {\widehat{X}}_{\text{i}}\text{(j}\omega \text{)}$$ . Afterwards the spectrum is filtered by a frequency-domain lowpass filter. Ideal filtering can be achieved by means of a rectangular window with a desirable cutoff-frequency $$f_{\text{c}}$$ . The latter may be chosen (but is not limited to) to correspond the nominal fundamental frequency $$f_{\text{c}}\widehat{=}{n}_{\text{h}}{\omega }_{\text{f}}/2\pi$$ . The windowed and filtered spectrum $$X_{\text{i}}\text{(j}\omega \text{)}$$ is multiplied twice by $$\text{j}\omega$$

$$\begin{array}{l}{\dot{X}}_i(\text{j}\omega )=\text{j}\omega X_i(\text{j}\omega )\\ {\ddot{X}}_i(\text{j}\omega )=-{\omega }^2{X}_i(\text{j}\omega )\end{array}$$ E27

Transforming back to the time domain yields the filtered signals $${\dot{\widehat{x}}}_i={\text{DFT}}^{-1}\left({\dot{X}}_i(\text{j}\omega )\right)$$ and $${\ddot{\widehat{x}}}_i={\text{DFT}}^{-1}\left({\ddot{X}}_i(\text{j}\omega )\right)$$ . The posture estimate is also updated according to $${\widehat{x}}_i={\text{DFT}}^{-1}\left(X_i(\text{j}\omega )\right)$$ . The procedure of data processing in the frequency domain is depicted in Fig. 3. The filtered estimates of the velocities and accelerations of the end-effector are provided by using classic kinematic transformations (Merlet, 2006).

Figure 3.

Frequency-domain data processing for filtering and differentiation of non-measured signals

4.2. Parameter uncertainties

To validate the results of the identification, statements on the uncertainties of the obtained parameters are necessary. For the given linear model structure (20) and by assuming Gaussian disturbance vector $$\eta$$ see Abdellatif et al., 2005b), the covariance of the parameter estimate resulting from (21) is

$$P={({\Psi }^{\text{T}}{\Sigma }_{\eta \eta }^{-1}\Psi }^{-1}$$ E28

The confidence area of the estimated parameter set $${\widehat{p}}_{\text{GM}}$$ with respect to the unknown true parameter vector $$p$$ can be calculated for a given quantile $$\alpha \in \left[\text{0}\dots \text{1}\right]$$ as a $$100(1-\alpha )\%$$ confidence ellipsoid:

$${\epsilon }_{\alpha }=\left\{p\in {\Re }^{n_{\text{p}}},{\left(p\text{-}{\widehat{p}}_{\text{GM}}\right)}^{\text{T}}{P}^{-1}\left(p\text{-}{\widehat{p}}_{\text{GM}}\right)\le {\chi }_{\alpha }^2\left(n_{\text{p}}\right)\right\}$$ E29

where $${\chi }_{\alpha }^2\left(n_{\text{p}}\right)$$ denotes the value of the $${\chi }^2$$ distribution with $$n_{\text{p}}$$ degrees of freedom at the quantile $$\alpha$$ and $$n_{\text{p}}$$ is the dimension of the parameter vector (Gevers, 2005). Consequently, the estimate of the single parameter $${\widehat{p}}_k$$ is normally distributed $$N\left(p_k,P_{kk}\right)$$ with variance $$P_{kk}^2$$ , where $$p_k$$ is the true parameter value and $$P_{kk}$$ is the k^th diagonal element of $$P$$ . A 95% confidence interval can be determined as

$$C_k^{95\%}=\left[{\widehat{p}}_k-2{\rho }_k,{\widehat{p}}_k+2{\rho }_k\right]=\left[{\widehat{p}}_k-2\sqrt{P_{kk}},{\widehat{p}}_k+2\sqrt{P_{kk}}\right]$$ E30

Equations (29) and (30) are useful to evaluate the confidence of the estimate results for the complete parameter set or for the single parameters, respectively.

5. Experimental results for model-based control

This section is dedicated to the experimental results achieved on the hexapod PaLiDA.

5.1. Description and modelling of the hexapod

The parallel robot PaLiDA (see Fig. 1) was developed by the Institute of Production Engineering and Machine Tools at the University of Hannover as a Stewart–Gough platform. It is designed with electromagnetic linear direct drives used as extensible struts for use in fast handling and light cutting machining like deburring. The actuation principle has several advantages compared to conventional ball screw drives: Fewer mechanical components, no backlash, low inertia with a minimized number of wear parts. Furthermore, higher control bandwidth and extremely high accelerations can be achieved. A commercial electromagnetic linear motor originally designed for fast lifting motions is improved for use in the struts. Each strut of the hexapod is composed of three bodies as depicted in Fig. 4. Thus, the system is modelled with 19 bodies: The movable platform (index E), 6 identical movable cardan rings (index 1), 6 identical stators (index 2) and 6 identical sliders (index 3).

Figure 4.

MDH-frames and parameters of the struts

The dynamics model in parameterlinear form results by applying the rules discussed in section 2. The rigid-body part contains 10 base parameters (see Table 1). According to the friction modelling approach (14) the actuated joints $$q_{\text{a}}$$ correspond to 6 different dry friction and also 6 different viscous damping coefficients. Friction in the passive joints is modelled only as dry friction with a common parameter for all $${\alpha }_j$$ and another one for all $${\beta }_j$$ -joints. The friction model contains therefore 14 different parameters. Its structure was optimized according to the statistical analysis mentioned in section 2.2 and presented in (Abdellatif et al., 2005c).

rigid-body	friction
$$p_1=I_{{\text{zz}}_1}+I_{{\text{yy}}_2}+I_{{\text{zz}}_3}\left[{\text{kgm}}^{\text{2}}\right]$$	$$p_{11}=r_{\alpha }\left[\text{N}\right]$$
$$p_2=I_{{\text{xx}}_2}+I_{{\text{xx}}_3}-I_{{\text{yy}}_2}-I_{{\text{zz}}_3}\left[{\text{kgm}}^{\text{2}}\right]$$	$$p_{12}=r_{\beta }\left[\text{N}\right]$$
$$p_3=I_{{\text{zz}}_2}+I_{{\text{yy}}_3}\left[{\text{kgm}}^{\text{2}}\right]$$	$$p_{13}=r_{1_1}\left[\text{N}\right]$$
$$p_4=s_{{\text{y}}_2}\left[\text{kgm}\right]$$	$$p_{14}=r_{1_2}\left[\text{N}\right]$$
$$p_5=s_{{\text{y}}_3}\left[\text{kgm}\right]$$	$$p_{15}=r_{1_3}\left[\text{N}\right]$$
$$p_6=I_{{\text{xx}}_{\text{E}}}+m_3{\displaystyle {\sum }_{j=1}^6\left(r_{{\text{B}}_{{\text{y}}_j}}^2+r_{{\text{B}}_{{\text{z}}_j}}^2\right)}\left[{\text{kgm}}^{\text{2}}\right]$$	$$p_{16}=r_{1_4}\left[\text{N}\right]$$
$$p_7=I_{{\text{yy}}_{\text{E}}}+m_3{\displaystyle {\sum }_{j=1}^6\left(r_{{\text{B}}_{{\text{x}}_j}}^2+r_{{\text{B}}_{{\text{z}}_j}}^2\right)}\left[{\text{kgm}}^{\text{2}}\right]$$	$$p_{17}=r_{1_5}\left[\text{N}\right]$$
$$p_8=I_{{\text{zz}}_{\text{E}}}+m_3{\displaystyle {\sum }_{j=1}^6\left(r_{{\text{B}}_{{\text{x}}_j}}^2+r_{{\text{B}}_{{\text{y}}_j}}^2\right)}\left[{\text{kgm}}^{\text{2}}\right]$$	$$p_{18}=r_{1_6}\left[\text{N}\right]$$
$$p_9=s_{{\text{z}}_{\text{E}}}+m_3{\displaystyle {\sum }_{j=1}^6{\text{r}}_{{\text{B}}_{{\text{z}}_j}}}\left[\text{kgm}\right]$$	$$p_{19}=r_{2_1}\left[{\text{Nsm}}^{\text{-1}}\right]$$
$$p_{10}=m_{\text{E}}+6m_3\left[\text{kg}\right]$$	$$p_{20}=r_{2_2}\left[{\text{Nsm}}^{\text{-1}}\right]$$
	$$p_{21}=r_{2_3}\left[{\text{Nsm}}^{\text{-1}}\right]$$
	$$p_{22}=r_{2_4}\left[{\text{Nsm}}^{\text{-1}}\right]$$
	$$p_{23}=r_{2_5}\left[{\text{Nsm}}^{\text{-1}}\right]$$
	$$p_{24}=r_{2_6}\left[{\text{Nsm}}^{\text{-1}}\right]$$

Table 1.

Rigid-body and friction model parameters for the parallel robot PaLiDA

5.2. Experiment design and data processing

The experiment design has been carried out according to the method given in section 3. An example of a resulting excitation trajectory with the order $$n_{\text{h}}=5$$ is depicted by Fig. 5. The obtained measurements of the actuator lengths are transformed numerically by the direct kinematics. The resulting estimation of posture elements are then filtered and differentiated in the frequency domain as proposed in section 4.1. Fig. 6 illustrates exemplarily such procedure for the reconstruction of the second translational degree of freedom $$y$$ corresponding to the excitation trajectory, shown in Fig. 5.

The left side of Fig. 5 depicts the frequency-discrete spectral amplitudes of the signals along with the used selection window that corresponds to an ideal lowpass filter. The respective signals in the time-domain are given on the right side of the picture. The effectiveness of the proposed filter is obvious, since the calculated signals exhibit almost no noise or disturbance corruption. Such property is a central requirement for a robust and reliable identification of parallel manipulators, because the necessary but non-measurable information has to be extracted from corrupted and limited measurements of the actuator displacements.

Figure 5.

Example of a periodic excitation trajectory $$n_{\text{h}}=5$$ ; top left: Translational coordinates, bottom left: Rotational coordinates, right: 3-D presentation

Figure 6.

Reconstruction of the end-effector displacement, velocity, and acceleration, with respect to the inertial $$y$$ -axis by using frequency-domain filtering

In the following three models are compared, that all result from the identification using the same trajectory but after implementing three different data-processing techniques. The first one results directly from rough data without any filtering. For the second, the measurements of the actuator displacements were filtered in the time domain. The third model has been identified according to the proposed frequency domain method. The validation of the models on a circular bench-mark trajectory, that was not used for identification, is depicted in Fig. 7. The frequency-domain processing yields the best prediction quality corresponding to the smallest error variance $${\sigma }^2$$ . Time-domain filtering is not accurate enough to extract all information at the relevant frequencies.

Figure 7.

Prediction Accuracy of three different models for an arbitrarily chosen actuator; left: By using rough data, middle: By using time-domain filtering, right: By using frequency-domain filtering.

5.3. Estimation results and parameter uncertainties

The filtered data resulting from the investigated trajectory (Fig. 5) are used to compute the regressor matrix $$\Psi$$ . The corresponding actuation forces can be obtained from the measurement of the motor currents. The case of PaLiDA reveals high noisy and cross-correlated force measurements (Abdellatif et al., 2005b). Therefore, the Gauss-Markov estimate has been proposed earlier (see (21)) that yields the parameter set given in Table 2. It is important to notice, that the provided a priori values do not present the true parameters, since they were calculated by using uncertain CAD-Data. The quality of the results is in general very high, despite that the parameters with small values exhibit higher uncertainties. This is however a known and general problem of experimental estimation in practice. We refer to former publications for detailed discussions on the different aspects of the estimation results (Abdellatif et al., 2005b, Abdellatif et al., 2005c).

The validation of the parameter estimation robustness can be provided, e.g. after repeating the identification experiment 100 times. The resulted parameter sets are compared to the 95% confidence intervals (see eq. (30)). Such investigation is depicted for some exemplarily chosen parameters in Fig. 8. The history of the weighted parameter estimate $${\widehat{p}}_i/{\overline{p}}_i$$ are illustrated over the measurement trials, where $${\overline{p}}_k$$ is the mean value of all estimates. The corresponding weighted upper and lower bounds $$C_{\text{M}}$$ and $$C_{\text{m}}$$ of the confidence intervals are additionally shown. The robustness of the identification is proven, since the estimates remain mostly within the confidence intervals. Some exceptions are observed though, such the very small first rigid-body parameter and the first few measurement trials. The latter is

$$p_k$$	$${\widehat{p}}_k$$	$${\rho }_k=\sqrt{P_{\text{kk}}}$$	$$C_k^{95\%}$$	a priori
$$p_1\left[{\text{kgm}}^{\text{2}}\right]$$	-0.0447	0.0039	[-0.0526 -0.0369]	0.0074
$$p_2\left[{\text{kgm}}^{\text{2}}\right]$$	1.0892	0.0070	[1.0753 1.1032]	0.9439
$$p_3\left[{\text{kgm}}^{\text{2}}\right]$$	1.0077	0.0045	[0.9988 1.0166]	0.9458
$$p_4\left[\text{kgm}\right]$$	0.5995	0.0036	[0.5922 0.6068]	0.6201
$$p_5\left[\text{kgm}\right]$$	-1.2885	0.0056	[-1.2998 -1.2772]	1.2295
$$p_6\left[{\text{kgm}}^{\text{2}}\right]$$	0.3078	0.0061	[0.3049 0.3106]	0.2878
$$p_7\left[{\text{kgm}}^{\text{2}}\right]$$	0.3021	0.0014	[0.2996 0.3045]	0.2878
$$p_8\left[{\text{kgm}}^{\text{2}}\right]$$	0.1176	0.0012	[0.1152 0.1201]	0.1217
$$p_9\left[\text{kgm}\right]$$	1.8896	0.0012	[1.8774 1.9017]	1.9012
$$p_{10}\left[\text{kg}\right]$$	16.3081	0.0460	[16.2161 16.4002]	16.1920
$$r_{\alpha }\left[\text{Nm}\right]$$	0.5756	0.0158	[0.5440 0.6072]	-
$$r_{\beta }\left[\text{Nm}\right]$$	0.9195	0.0179	[0.8837 0.9552]	-
$$r_{1_1}\left[\text{N}\right]$$	11.9772	0.2485	[11.4803 12.4742]	-
$$r_{1_2}\left[\text{N}\right]$$	4.8071	0.1861	[4.4350 5.1793]	-
$$r_{1_3}\left[\text{N}\right]$$	20.1528	0.3226	[19.5075 20.7980]	-
$$r_{1_4}\left[\text{N}\right]$$	5.1518	0.1817	[4.7884 5.5151]	-
$$r_{1_5}\left[\text{N}\right]$$	1.5857	0.2618	[1.0620 2.1094]	-
$$r_{1_6}\left[\text{N}\right]$$	5.0057	0.3519	[4.3018 5.7096]	-
$$r_{2_1}\left[{\text{Nsm}}^{\text{-1}}\right]$$	16.8771	0.5268	[15.8235 17.9307]	-
$$r_{2_2}\left[{\text{Nsm}}^{\text{-1}}\right]$$	16.7406	0.3712	[15.9981 17.4830]	-
$$r_{2_3}\left[{\text{Nsm}}^{\text{-1}}\right]$$	6.3408	0.5720	[5.1968 7.4848]	-
$$r_{2_4}\left[{\text{Nsm}}^{\text{-1}}\right]$$	23.1662	0.3799	[22.4065 23.9259]	-
$$r_{2_5}\left[{\text{Nsm}}^{\text{-1}}\right]$$	26.4675	0.4461	[25.5754 27.3596]	-
$$r_{2_6}\left[{\text{Nsm}}^{\text{-1}}\right]$$	22.8053	0.5539	[21.6974 23.9131]	-

Table 2.

Estimated dynamics parameters of the hexapod parameters $$\widehat{p}$$ with corresponding standard deviations, confidence intervals and a priori values for the rigid-body model parameters

due to the variation of friction at the beginning of the measurement process until a nearly stationary state is reached. Additionally to the single parameters, the confidence of the entire parameter set can be validated. The outer bound of the 95% confidence ellipsoid $${\epsilon }_{5\%}$$ is given by $${\chi }_{5\%}^2\left(\dim \left(p\right)=24\right)=36.42$$ . Its comparison with distribution $${\chi }_{}^2\left({\widehat{p}}_i\right)$$ of the vector estimates $${\widehat{p}}_i$$ over the measurement trials is given by Fig. 9. Excepting the first trial, the set of all parameters lays clearly within the confidence ellipsoid, which demonstrates the effectiveness and robustness of the estimation.

Figure 8.

History of parameter estimates over different measurement trials: For the sake of uniform illustration, the parameters are given as weighted terms with respect to their respective mean values.

Figure 9.

Comparison of the $${\chi }_{}^2$$ distribution of the estimated parameter sets with the radius $${\chi }_{5\%}^2\left(24\right)$$ of the 95% confidence ellipsoid

Figure 10.

Control errors for both test trajectories at increasing end-effector velocity and by implementing different control strategies

5.4. Identification and model-based control

The impact of identification on the control and tracking accuracy of the hexapod PaLiDA is studied in the following. Hereby three control strategies are investigated. The first variation passes on any model knowledge, i.e. by implementing only linear controller for the single actuators. The second uses the inverse dynamics model to compensate for the nonlinear dynamics by considering only nominal parameter values. The third variation uses the identified model for the feedforward compensating control. All approaches are substantiated experimentally on two different trajectories: The first trajectory is a circular one and allows reaching high actuation forces, whereas the second is quadratic and is characterized with high actuator velocities.

Both trajectories were executed at different velocities $$v_{\text{E}}$$ of the end-effector and for the three mentioned control variations. The evolution of the rooted mean squares errors $$e_{\text{RMS}}$$ of all actuator deviations is depicted in Fig. 10 with respect to $$v_{\text{E}}$$ .

Figure 11.

Tracking accuracy of actuator 5 for the two studied trajectories at maximal velocity; comparison between the compensation of nominal model (thin line) and identified model (thick line)

As expected, the use of standard linear control (variation 1) exhibits a significant decreasing accuracy with increasing speeds, since the impact of nonlinear and coupled dynamics increases with higher velocities and accelerations. Using model-knowledge (variation 2 and 3) improves always the tracking performance. Furthermore, the compensation of identified model (variation 3) outperforms clearly variation 2 that just uses the nominal parameter values. The latter statement can be proven at the level of actuator tracking accuracy like depicted in Fig. 11. For the same arbitrarily chosen actuator, the tracking accuracy is higher if the identified model is implemented. The same results are noticeable for the cartesian tracking accuracy $$\Delta x$$ , like depicted in Fig. 12. It may be concluded that only accurately identified model allows keeping good tracking performance over a wide range of the robot dynamics.

Figure 12.

Calculated cartesian tracking accuracy $$\Delta x$$ for the two studied trajectories at maximal velocity; comparison between the compensation of nominal model (thin line) and identified model (thick line)

6. Conclusions

The present chapter discussed most significant aspects to achieve accurate and robust dynamics identification for parallel manipulators with 6 dof's. Hereby, the adequate consideration of structural properties of such systems has been stressed out. First, an efficient methodology to determine the inverse dynamics in a parameterlinear form has been presented, which enables the use of linear estimation techniques. Periodic excitation has been proved to be a powerful method for parallel robots, since it allows for appropriate consideration of hard workspace constraints. Due to measurement noise and cross coupling between the actuators, the achievement of the identification in a statistical framework is recommended. This includes the consideration of disturbance covariances in the experiment design, the use of Gauss-Markov estimation approach as well as the frequency-domain filtering to extract non measurable information from rough data. The robustness of the results has been substantiated on a direct driven hexapod. The obtained estimates have presented high confidence in terms of single parameters, as well as in terms of the whole parameter set. Additionally, the benefits of accurate identification on the enhancement of control performance have been clearly and experimentally demonstrated.