Model Calibration for a Rigid Hexapod-Based End-Effector with Integrated Force Sensors

Integrated single-axis force sensors allow an accurate and cost-efficient force measurement in 6 degrees of freedom for hexapod structures and kinematics. Depending on the sensor placement, the measurement is affected by internal forces that need to be compensated for by a measurement model. Since the parameters of the model can change during machine usage, a fast and easy calibration procedure is requested. This work studies parameter identification procedures for force measurement models on the example of a rigid hexapod-based end-effector. First, measurement and identification models are presented and parameter sensitivities are analysed. Next, two excitation strategies are applied and discussed: identification from quasi-static poses and identification from accelerated continuous trajectories. Both poses and trajectories are optimized by different criteria and evaluated in comparison. Finally, the procedures are validated by experimental studies with reference payloads. In conclusion, both strategies allow accurate parameter identification within a fast procedure in an operational machine state.


Introduction
In production, many manufacturing applications, such as process diagnosis and quality assurance, or adaptive process control for milling, grinding or thermo-smoothing, require in-process force measurement, where, in particular, the measurement of spatial forces and moments in 6 degrees of freedom (DoF) is requested [1]. As one approach, Refs. [1][2][3] present a new measuring system using a hexapod-based end-effector with structure-integrated sensors. Six 1 DoF force sensors are integrated into a rigid bar framework and the measured forces are transformed to Cartesian forces and moments at the tool centre point (TCP) by a control-integrated model.
Compared to commercial 6 DoF force/torque transducers (FT sensors), which are mounted close to the TCP, the new measuring system does not reduce usable workspace, causes no restrictions to spindle mounting and is also less expensive. Compared to a drive current based force measurement in spindle drives, it is far more accurate (factor 100 smaller uncertainty [2]) and requires less modelling effort, since no non-linear and stochastic mechanical influences of drive components interfere such as friction, stick-slip effects or elastic deformation. Altogether, structure-integrated force measurement has many benefits and can complement or replace the measurement by the use of conventional FT sensors as long as it provides a 6 DoF measurement, fulfils accuracy and stiffness requirements, and is simple to parametrise [1]. In addition, integrated sensors can be supplied by the machine manufacturer and make force measurement and control an easyto-use machine-integrated feature in the future. However, to reach this point, further studies and evidence regarding the parametrisation are required, as hardly any suitable literature on parameter identification for force measurement models is available.
Machine-integrated force measurement solutions have been studied since the beginning of force measurement and are outlined in a detailed state-of-the-art overview in [1]. In states. The larger the distance between TCP and sensors within the machine structure, the more parameters need to be taken into account, which, as a result, lead to higher modelling and parametrisation efforts, and higher measurement uncertainties. A proof of concept for all setups evaluating quasi-static measurements is presented in [2], which shows models, classification, parameter identification, sensitivity analysis, and experimental results. A more detailed view of the end-effector, including the dynamic model as well as accuracy and precision analyses during quasi-static, dynamic, and milling process force measurements, is discussed in [1]. Kinematic-integrated approaches (setups 2x) and force control using these measuring systems are topics of future work.   Figure 1. Force sensor integration into hexapod structures and kinematics: approach and sensor placement; 1a: Rigid not moving hexapod structure (e.g., clamping table), 1b: Rigid moving hexapod structure (e.g., end-effector platform), 2x: Hexapod kinematics with sensors in the struts at various positions [1][2][3]. However, every model needs a parametrisation and to make the promising results of [1] available for practical purposes, this paper applies and evaluates parameter identification procedures on the example of the rigid hexapod-based end-effector (setup 1b), Figures 1 and 2. Nevertheless, the introduced methods are also applicable for the remaining structure-integrated setups and even for the use of FT sensors.
For all approaches, the parameters contain sensor positions and orientations (36 parameters, constant for rigid frameworks) and sensor zero offsets (6 parameters). Additionally, for the focused moving rigid bar structures (setup 1b), mass, centre of gravity and inertia of the framework (10 parameters), as well as its actual pose, velocity and acceleration in real-time, are included. Moreover, additional parameters need to be introduced to include influences concerning the sensibility of force sensors towards lateral forces, torques, elastic structure deformation or temperature effects. Table 1 gives an overview of the presented steps and variations according to the methods and the structure of the approach: Based on the measurement model (Section 3), at first, different identification models are derived, which include rigid body, sensor and error compensation parameters (Section 4). As the description and optimisation of the parameter set are crucial for identification, next, multiple parametrisations are analysed and compared using measurement simulation, Pearson correlation and QR decomposition (Section 5). Hereby, two excitation approaches are presented: a quasi-static measurement at independent poses (Section 6), and a dynamic excitation using continuous trajectories (Section 7). The two excitations are created both, in a manual way and in an algorithmic way by the use of different optimisation criteria and methods. Finally, both methods are compared in terms of the application and identification of reference payloads (Section 8). Figure 2. Left: Hexapod machine tool with rigid hexapod-based end-effector platform (setup 1b) including six structureintegrated 1 DoF force sensors of the type Althen ALF256 5 kN. Right: Corresponding values for measurement model. The integration of the sensors, which are screwed in after cutting, drilling, and thread-cutting the bars of the existing end-effector, is fast, cheap, and does not change mass or size of the end-effector, significantly.

Measurement Model
The measurement model in the current form was first presented in [1] and is summarised here, since it is essential for the identification model. Solving the equations of motion ∂p/∂t = ∑ f at the platform frame {P}, Figure 2 right, gives external spatial forces and moments f Ext that act on the end-effector on the basis of the measured and processed sensor forces f * q : The sensor forces are transformed to Cartesian space using the geometric Jacobian P J P of the platform structure, Figure 2 right, and gravity forces P f G are included by with the gravitational acceleration B g 0 as well as mass m P , orientation B R P , and centre of gravity P x C of the platform. Underlines indicate 6 DoF vectors; upright letters constants. The spatial dynamic force vector P f Dyn divides in forces P f Dyn and moments P m Dyn : with inertia P I P = P I P,C + m P S( P x C ) S( P x C ), angular velocity B ω P , and acceleration Bω P of the platform in relation to the base frame {B}. Here, S(r) is the cross product matrix: With the platform position B o P , the centre of gravity acceleration Bẍ C is The platform pose x calculates from measured drive positions q in forward kinematics (FK) of the hexapod and the platform velocities v from measured drive velocitiesq in differential kinematics with the help of the geometric Jacobian of the hexapod J Hex (x) Numeric differentiation ofq and symbolic differentiation of J Hex give the acceleratioṅ v of the platformv = BöT P BωT P Finally, this leads to the measurement model in its canonical form, with P ω = B R T P B ω P , Here, f * q stands for compensated sensor values that already contain corrections and calculates to For the present setup, these are sensor zero offsets f q0 and pose-dependent influences f K (x), which address lateral sensor forces and moments from platform deformation, as described next.
According to the rigid body assumption of the end-effector, the force measurement is expected to be independent of the platform position, Equation (10). However, experiments show a not negligible position-dependent force deviations in a range of ±20 N (within a measuring range of ≈9 kN in XY) that are strictly systematic and mainly dependent on Xand Y-coordinates [1]. Ref. [2] evaluates the platform deformation that results from the 36 varying platform joint forces and moments, which are applied by the hexapod kinematic struts, as reasonable explanation for this behaviour through an analysis based on a finite element model (FEM). As a full elastic model has disadvantages in parametrisation and online calculation in the control system, the effect is compensated by a prediction model with the choice of where the hexapod pose is used to describe the current static load situation of platform forces. This approach produces two linear factors for each sensor, that can be identified by quasi-static measurements. Evaluations show that these 12 parameters do not represent a minimal model. Transforming K 1 to Cartesian, reveals a systematic force deviation that can be described in a reduced form as a second modelling approach Now, the effect is compensated for by only two parameters {k 0 , k 1 }. In fact, platform deformations resulting from platform joint forces not only occur at position changes but also at orientation changes. Neglecting these influences leads to the incorrect identification of mass and centre of gravity by quasi-static measurements, as shown below. As the orientation is part of the rigid body gravity model, Equation (3), a compensation approach will correlate with the parameters mass m P and centre of gravity P x C . Therefore, identification is not possible when the force offsets are modelled in any formulation using the orientation matrix, such as f K3 To separate the parameters, a different angle format can be used for the orientation error model, Equation (14), as an extension to Equation (12) Again, a transformation to Cartesian shows that the parameters can be reduced to the two most significant influences k 2 and k 3 For this notation, an axis based format-here the Gibbs vector-needs to be used. A change of the Jacobian due to platform deformation can be neglected.
In conclusion, by those four parameters, {k 0 , k 1 , k 2 , k 3 }, the pose-dependent force deviations are compensated for by more than factor 10 to a remaining total error of approximately ±2 N (3 s, 99.73%) in XY for the measurement at quasi-static poses in the workspace [1].
Alternatively, by inserting joints [12,15] or flexure hinges [13,25] into the structure, forces lateral to the sensors can be reduced mechanically. Nevertheless, they cannot remove longitudinal forces that result from a deformation of the fundamental frame, which here is the platform part that deforms due to forces applied by the hexapod struts. Further, joints and flexure hinges introduce backlash or reduce the platform stiffness respectively, which is not desired in machine tools.

Identification Model
Besides process force measurement and control, structure-integrated force sensors can be used for parameter identification: As mass, centre of gravity and inertia of the end-effector can change during machine usage, for example through a workpiece or tool change, a fast and easy procedure for identifying these parameters is required, which should last less than 5 min and can run in a production break. Other parameters, such as sensor offsets and error compensation factors, are only affected marginally by load changes but need to be determined at least once after system assembly. Generally, it is advantageous to identify all parameters at the same time to achieve the least global variance. If the system has been identified, a pure sensor drift can be tared without further measurements by adjusting only f q0 , Equation (11).
Since friction and stick-slip effects can be neglected for the sensor placement in the rigid end-effector without joints (setup 1b), an identification from quasi-static measurements is possible: Doing this, the 10 parameters, platform mass m P , centre of gravity P x C and force sensor offsets f q0 , as well as the pose-dependent compensation offsets K can be identified. This procedure is appropriate for quasi-static or continuous applications that do not require inertia, which cannot be identified statically. However, limited dynamics apply for most applications in force measurement and control to reduce measuring noise. On the other hand, a measurement during accelerated machine movement allows a parameter identification from dynamic excitations and, additionally, the identification of the platform inertia P I P,C . With this, both standard approaches to identification are discussed: the quasistatic identification from independent poses and the identification from dynamic excitation trajectories with dependent data points. The Jacobian, Equation (2), is not identified since it is, as a pure geometrical transformation, already known from mechanical drawings. In the rare case when it is unknown, it can be determined directly by geometric measurements far more accurately than indirectly from dynamic identification.
Below, a linear identification model of the form is introduced for the complete system. A partial identification model, for example, without inertia, can be achieved by elimination of the corresponding columns in the regressor matrix A and rows in u. First, the vector u of u unknowns is set to where matrices have been transformed to vectors, which is indicated by an overline: P I P = {I 11 , I 12 , I 13 , I 22 , I 23 , I 33 }. Rearranging the model from Equation (1) to express n virtual sensor forcesf q gives the non-linear identification model The identification procedure can be performed with the non-linear model, Equation (18), by calculating its local Jacobian and iterating over it or by transforming Equation (18) into a linear form analytically, which is advantageous, as no iterations are needed to solve the problem. The linear form of the rigid body part can directly be obtained from Equation (10) with the help of the operator L(ω) [22,24] Together with the measurement parameters, this gives the linear form for one dataset where 1 is the identity matrix andK represents compensation factors of f K . In the case of Equation (12)K isK 1 =L(x j ), whereL provides a transformation K 1 x =L(x) K 1 similar to the inertia linearisation. In the case of Equation (13)K is with the unknowns K 2 = {k 0 , k 1 }, and for K 3 , K 4 , and K 5 in a similar way. Stacking m measurementsf qj and the corresponding models A j to the overdetermined global system leads tof which can be solved by the use of a Gauss-Markov approacĥ with the weighting matrix W = cov( f q ) −1 [26]. The weighting matrix can be used to take varying characteristics of the measurements into account, regarding different units or the individual measuring accuracy of different sensors. As only forces are measured and identical sensors are used in the present setup (ALF256 5 kN, Althen GmbH Mess-und Sensortechnik, Kelkheim, Germany), task variable scaling is not necessary and W is set to identity, which reduces the problem to an ordinary least-squares regression.
In contrast, parameter scaling is still required to obtain meaningful and comparable singular values for the next steps model analysis and optimisation as well as improved convergence. With the scaling matrix H = diag(h 1 . . . h u ), the model changes to [24] f q = (AH)(H −1 u) =Ãũ + ε. (25) For the choice of H, column scaling is used, as it is simple and does not require a priori statistical data, with a l as lth column of A [27] An estimate of the global standard deviation of the regression problem can be obtained from the χ 2 statistic [24,26] the number of poses m, the number of forces measured per pose n, and the number of unknowns u. The covariance matrix as well as local standard deviations of the unknowns givê

Model and Identifiability Analysis
Generally, only those parameters that are modelled are identifiable, stimulated by the exciting trajectory (sensitivity), and that can be distinguished from each other (linear independence). Therefore, the optimisation of the parameter set and the optimisation of the exciting data points or trajectory parameters, respectively, are the two major steps in identification preparation. The parameters can be grouped in identifiable and unidentifiable parameters as well as parameters that are only identifiable in linear combination, expressed by a new parameter. QR or singular value decomposition ofÃ helps in analysing the identifiability of the parameter set [24].
For the end-effector, different models are proposed for identification with different parameter sets and properties (Table 2). Identification results for these models are shown in Table 3. Generally, the inertia P I P is not identifiable by quasi-static measurements. Model 0 contains the basic static rigid body parameters and sensor offsets but does not include error correction factors. Model 1, Equation (12), and model 2, Equation (13), include positiondependent error correction, which increases the measurement accuracy significantly [1], where model 2 describes this effect with fewer parameters, as described above.  Table 3; In squares: selected models for excitation optimisation.

Model & Parameters No. Comment
Quasi-static Identification Full error correction: correct m P and P x C . 5 m P , P x C , f q0 , K 2 , K 5 14 Full error correction: correct m P and P x C .
Dynamic Identification 6 m P , P x C , f q0 , P I P 16 No error correction. 7 m P , P x C , f q0 , P I P , K 2 18 XY error correction. 8 m P , P x C , f q0 , P I P , K 2 , K 4 36 Full error correction, parameters not minimal 9 m P , P x C , f q0 , P I P , K 2 , K 5 20 Full error correction. For all static models that neglect a modelling of orientation-dependent force errors (models 0-2), the parameters m * P and P x * C do not represent the physical platform mass and centre of gravity, which results in an empirical model rather than in a physical model. Still, the models are valid because these parameters include the mentioned effects. Changes of mass and centre of gravity, on the other hand, can be identified with high precision and the correct physical values [2]. This is important, as tool or workpiece parameters are represented with physical meaning.
Including orientation-dependent force errors in the models always leads to a correlation of the error model parameters to the rigid body parameters, mass and centre of gravity, for quasi-static measurements since both parameter sets can only be excited by inclinations. When modelled by rotation matrix, Equation (15), the parameters are not identifiable as the mathematical representation is too close to the rigid body gravity model in Equation (3).  (not identifable in any description based on R T g), centre: preferred quasi-static model 5 (includes models 0, 2), right: preferred dynamic model 9 (includes models 6, 7). The models 1, 4, 8 are not displayed due to too many parameters. Corr shows the pairwise linear correlation between 2 parameters, that points out linear dependencies. In R high values in the row behind the diagonal element indicate dependent parameters. The diagonal elements R ii of R that are less than the numerical zero (m · n · · max |R ii |) indicate unidentifiable parameters u i [24].
Including orientation-dependent force errors, based on a different orientation representation, leads to the model 4, Equation (14), and model 5, Equation (15), which allows for the identification of mass and centre of gravity with correct physical meaning, where model 5 describes the effect with less parameters.Still, the correlation between these parameters is a physical fact that cannot be completely overcome by a different modelling for quasi-static measurements so that the uncertainties stay high, as can be seen in Figure 3 centre.
Even though models 0-2 are valid and contain fewer parameters, a structural rigid body approach using parameters with physical meaning is advantageous for a second reason: as Table 2 shows, dynamic identification results in correct rigid body parameters for models 6 and 7 without orientation-dependent force error models. This is because a dynamic trajectory excites the parameters in a different way, while the platform joint forces that cause the force errors through platform deformation are foremost inducted by static hexapod strut forces-dynamic joint forces are much smaller and can be neglected, as separate simulations show. For the same reason, model 3 becomes identifiable for dynamic excitations when expanded by an inertia model (not displayed).
Finally, models 8 and 9 include all error compensation parameters for dynamic identification, where model 9 uses fewer parameters (Figure 3 right).
For a symmetric platform without accessories, some centre of gravity parameters { P x Cx , P x Cy } and moments of deviation { P I xy , P I xz , P I yz } may also be eliminated, which reduces the model by 5 more parameters. On the other hand, in practice, a workpiece can be applied at any position at the end-effector making these parameters necessary, which is why this reduction is not performed in the present case.
In conclusion, models 5 and 9 provide a full error correction with minimal parameters and are chosen as the basis for further steps in the evaluation of quasi-static and dynamic identification procedures, respectively.
In order to optimise poses and exciting trajectories, criteria are needed that can be defined on the basis of the regression matrix A and on the inverse Fisher information matrix as presented in Table 4. Table 4. Observability indicies and optimization criteria [24,28].

Quasi-Static Identification Measurements
For quasi-static identification measurements, several poses are approached in the workspace and, at each pose, a set of force samples is taken in a stationary condition. The ability to wait until dynamic effects extinguish, and to average over many samples, leads to a higher data point quality for quasi-static measurements when compared to dynamic measurements. At the same time, the data points are independent so that every data point is valid; the data do not need to be cut and, generally, fewer data points are acquired, which is advantageous concerning memory and processing time of the regression. Still, the measuring time is very high and the inertia cannot be identified. Depending on the model, a minimum of m = u/n poses is required, where an overdetermined measurement including approximately 10 times more poses than necessary aims to obtain regression results with acceptable uncertainty. The choice of the correct poses is crucial for good identification results, where different parameters are excited by different poses: The identification of {m P , P x C } and {k 2 , k 3 } requires large inclinations to excite varying gravity forces in Equation (3) and platform joint forces compensated by Equation (15), respectively. The identification of {k 0 , k 1 }, Equation (13), requires large position changes in XY directions.
Optimising the pose set by selecting the most sensitive poses avoids redundancies and, in total, reduces both poses and measuring time, significantly. Basically, two approaches are available to optimise the dataset: As the data points are independent, poses can be added randomly until the chosen criterion saturates. Alternatively, starting from a generated good dataset, poses can be removed as long as the chosen criterion does not change significantly. Here, the latter approach is performed using the criteria shown in Table 4. For every pose set, one pose is deleted temporarily and the criterion is calculated. After examining all poses, the pose with the smallest change of the selected criterion is found as pose with least influence and deleted permanently. Finally, this procedure repeats until the algorithm reaches the desired number of poses (Algorithm 1). Table 5 shows results for a reduction to 30 poses starting from a working dataset of 480 generated poses that include position changes from −300 mm to 300 mm in steps of 100 mm in XY and, at each position, orientation changes from 10°to 20°, 0°to 270°and −15°to 15°in steps of 10°, 90°and 15°, respectively, for modified Euler angles [29], Figure 4 left. Table 5. Optimisation of quasi-static identification procedure with model 5 through a reduction of the number of poses to 30 starting from algorithmically generated 480 poses; last row: normal probability density plots for the parameters.  Algorithm 1: Selection of optimal poses using criteria of Table 5.
// Create a set of valid poses for x ← x min to x max by x step do for y ← y min to y max by y step do for z ← z min to z max by z step do for a ← a min to a max by a step do for b ← b min to b max by b step do for c ← c min to c max by c step do if is_valid(pose(x,y,z,a,b,c)) then pose_list.add(pose(x, y, z, a, b, c)); // Reduce from n to r optimal poses red_pose_list ← pose_list; for i ← 1 to n − r do crit0 ← calc_criteria(red_pose_list); for j ← 1 to length(red_pose_list) do test_list ← red_pose_list; test_list.remove_pose(j); crit(j) ← calc_criteria(test_list); red_pose_list.remove(index_o f (min(crit0 − crit))); This optimisation reduces the measuring time from 21 min for 480 poses to 81 s for 30 poses (with approximately 2.6 and 2.7 s per pose, respectively), which meets the requirement of a fast measurement in production breaks and, at the same time, realises an over-determination factor of 12.8. Further, the algorithm can be run until A loses rank, which gives the minimum number of possible poses and, therefore, proves the selection of good poses through the algorithm, as the number is close to or reaches the theoretical minimum of 3 poses for all criteria (Table 5). A reduction to less than 30 poses is possible but not recommended since it increases the parameter uncertainties.
All criteria reduce the data set to 30 valid poses. As expected, poses with high absolute position values and large inclinations are preferred, whereas centred poses and poses with small inclinations are neglected ( Figure 4 centre and right). For all parameters and criteria, the standard deviations increase by approximately factor 5. The most sensitive parameters regarding the data reduction are {k 2 , k 3 } and {m P , P x C }, as the standard deviations are already higher without reduction. This matches the expectations from the correlation and QR plots in Figure 3, which indicate this correlation. Even when the condition criteria K shows the smallest total standard deviationŝ 0 of the optimised solutions, the plots reveal a displacement of the averages for several parameters when compared to the notoptimised solution. However, all optimised solutions are valid and do not differ very much. Besides validating the approach, this indicates that in the presented case, the previous step-meaning the model optimisation-is far more important than the used criteria for pose optimisation.

Dynamic Identification Measurements
Dynamic parameter identification measurement stands for the online acquisition of data points during a highly accelerated movement of the end-effector. Using a dynamic excitation trajectory enables the collection of many data points in a short time and the identification of the platform inertia that cannot be identified by quasi-static measurements. On the other hand, if an individual data point of the trajectory has a lower quality, a pre-processing of the data may be necessary and the regression needs more calculation time owing to more data points when compared to quasi-static measurements. Filtering options are also limited, the reconstruction of the necessary acceleration can be challenging and a synchronous measurement of force and acceleration needs to be ensured beforehand. With the use of 16 Bit analogue inputs EL3742 and EL3702 (Beckhoff, Verl, Germany), which are connected to a local EtherCAT-Client EK1101 (Beckhoff, Verl, Germany) at the end-effector platform, a fast and synchronous data acquisition of the forces is realised. The oversampling functionality allows a data acquisition rate up to 100 kSamples/s (10 µs) that is used for 100 times averaging during a task cycle of 1 ms. Further, strut positions and velocities are acquired synchronously from a Sercos II drive bus, where velocities are differentiated to accelerations numerically. Finally, filter times of force and acceleration are adjusted to obtain a synchronous data acquisition, where the differentiator delay needs to be taken into account. Details of data acquisition and measurement results, for example, during a milling process, are presented in [1].
For the exciting trajectory, two approaches are presented: a manually created trajectory and an optimised trajectory based on periodic functions. The model is so simple that the individual parameter sensitivity with regard to specific movements can be estimated by examining the equations, compare Equation (10), Equation (13) and Equation (15). A trajectory that is sufficient for parameter identification can be found considering the following aspects: the parameters mass m P and centre of gravity P x C are excited best by accelerations, which can be created by gravitational acceleration using inclinations ( B R T P B g 0 ), as in the quasi-static approach, or by highly accelerated dynamic movements in translational directions Pv , where the latter is advantageous as it reaches higher magnitudes and avoids correlation to the error correction factors {k 2 , k 3 }. The inertia requires high angular accelerations Pω and the error correction factors {k 0 , k 1 } require large position changes with low acceleration.
This leads to a trajectory, as shown in Figure 5 left, that contains successive excitations in all Cartesian directions {x, y, z, α, β, γ} followed by slow inclinations, a fast 6 DoF motion, and 3 circles placed in the 3 Cartesian planes, where the latter has a large diameter to excite {k 0 , k 1 }. The measured accelerations show clear excitations in the respective directions, except for slower or more complex movements, where they become noisy.
Optimising a trajectory, on the other hand, requires criteria that are defined in Table 4 and, as the data points are not independent, a proper description of the trajectory that is smooth, closed and limited with respect to the workspace borders. Periodic Fourier series, as presented by [30] for serial kinematics and applied for parallel kinematics by [20], meet these requirements and are used here as well.
For every DoF i = 1 . . . 6 in the workspace, a movement is defined as follows: [20,30] The analytical availability of the derivatives is advantageous, as they can be used to filter the usually noisy velocity or acceleration measurements. Depending on the angle format, v andv can be obtained fromẋ andẍ by the well known transformation v = Tẋ (37) v = Tẍ +Ṫẋ, with (38) on the example of Roll-Pitch-Yaw angles {Φ, Θ, Ψ}. Further, the trajectory must be kept within valid bounds during optimisation. Basically, the hexapod movement is limited by strut lengths, velocities and accelerations, passive joint angles and singularities. Owing to the specific configuration of the hexapod (e.g., eccentric joints), limited workspace coordinates are sufficient to avoid critical joint angles and singularities with certainty for the present kinematic, as long as they are not chosen very close to the workspace border. Velocity and acceleration limits can also be expressed in workspace coordinates x min ≤ẋ(t) ≤ẋ max (41) x min ≤ẍ(t) ≤ẍ max .
Altogether, this approach significantly improves optimisation speed as it avoids the calculation of the inverse kinematic, differential kinematic and validation algorithms inside the optimisation problem. At the same time, the remaining workspace stays large enough to house an optimised trajectory.
Looking at the parameters, the trajectory consists of 6(2N + 1) + 1 variables i x 0 , i a l , i b l and ω that can be optimised, where N is the order of the Fourier series and needs to be chosen beforehand. As the trajectory can be centred in the workspace due to machine symmetry, the parameters i x 0 are set to zero, which reduces the problem by 6 parameters. Generally, fewer parameters is advantageous for the convergence of the optimisation algorithm, where low order Fourier series produce very simple trajectories, which make an optimisation pointless, for example, circles for order 1. On the other hand, high order trajectories may easily converge to a local minimum and, further, may contain many direction changes that can break the acceleration constraints or require a reduction of the total velocity, which also reduces the excitation effect of the trajectory. For these reasons, the order should be at least 3 (37 parameters) and not higher than 7 (85 parameters). The presence of very small parameter values can indicate that the order can be reduced. Optimisation calculations, performed using model 9 for the criteria (A, D, K, N) of Table 4, different orders (1-7) and different solvers ('interior-point', 'active-set' and 'sqp' of Matlab R2019a (The Mathworks Inc., Natick, MA, USA) fmincon optimiser for constrained nonlinear minimisation), show the best results using the 'active-set' solver, D-or N-optimal criterion and third to fifth order trajectories, such as those presented in Figure 5 right. Even with these optimisations, which are performed with equal start values, can result in similar trajectories for different solvers and criteria, there are differences in performance and quality. In detail, the 'active-set' converges in most cases and needs fewer iterations than 'sqp', which sometimes does not converge-'interior-point' does not converge in general. Regarding the criteria, the D-optimal criterion, which maximises the volume of the parameter uncertainty hyper-ellipsoid, converges best due to least iterations and uncertainties, as also predicted by other authors for this identification task [17,20]. Further, the N-optimal criterion, which can be seen as a combination of minimising the condition number and maximising the minimum singular value [24,28], performs quite well, as it produces valid excitation trajectories within few iterations and with good parameter identifiability. Finally, A-and K-criterion can also lead to valid trajectories but usually need more iterations and contain weak parameters indicated by QR-analysis and, therefore, are not recommended. Furthermore, the algorithm requires good initial parameters, as it may converge to local minima, which is indicated by different trajectory results on randomly generated start values. On the other hand, even local minima results lead to trajectories that are sufficient for parameter identification. For this reason, random initial values are acceptable when they significantly differ from zero and are within the trajectory bounds. Consequently, the optimisation algorithm using the mentioned settings leads to a fast convergence (mainly depending on solver and criteria), a limited trajectory with significant acceleration (mainly depending on start values and criteria) and significant trajectory parameters (mainly depending on order). In contrast to the quasi-static approach, both optimising the model and choosing the algorithm and criteria are essential for good results.
When comparing the measured accelerations of both trajectories ( Figure 5), the manually planned trajectory shows higher magnitudes and far less noise. An explanation for this behaviour can be assumed in two characteristics of the optimised path approach. First, the optimised trajectory is planned as a continuous spatial path in all DoF and is composed of many intermediate points, whereas the manually created trajectory behaves more like a point-to-point motion that realises the successive Cartesian excitations and contains fewer intermediate points. This is significant because the path planner and interpolator in the numeric control (CNC) pre-process the path under the condition of minimal contour deviations including a velocity adjustment during the look-ahead procedure and jerk limitation depending on contour artefacts, such as non tangential block transitions. By activating spline interpolation, allowing contour deviations and improving dynamic parametrisation, a significantly higher dynamic is achieved. Nevertheless, to find a good parametrisation of the CNC path planner for the optimised trajectories containing simultaneous translations and rotations in 6 DoF can be challenging for parallel kinematics. For point-to-point movements instead, the path can be planned with maximal acceleration in all DoF neglecting contour deviations. Second, by the use of 7 segment motion profiles for the point-to-point movements, higher accelerations can be reached in general, when compared to the sin/cos movements of the optimised Fourier series. Table 6 compares identification results for both dynamic trajectories. Both identify similar rigid body parameters and sensor offsets, where the manually chosen trajectory shows smaller uncertainties in both,ŝ 0 andŝ j , when compared to the periodic trajectory. Further data processing, such as the selection of highly accelerated motion sections or the filtering of the accelerations based on discrete Fourier transformation (DFT-filtering), does not improve the results for the presented application. As a consequence of less measuring noise in accelerations and smaller uncertainties during identification, the manually chosen trajectory is preferred.

Validation and Comparison
As a final step, the identification procedures and exciting poses/trajectories are validated and compared by the application and identification of benchmark payloads ( Figure 6). In detail, weights with a total mass of 100 kg are applied to the end-effector in steps of approximately 20 kg and parameters are identified using the presented methods of quasistatic identification (model 5, with all poses, Q5M, and with 30 N-optimal poses, Q5N) and dynamic identification (model 9, manually designed, D9M, and optimised periodic trajectory, D9N). Table 7 shows the identification results for selected parameters of all measurements as well as relative deviations to the theoretical physical values that result from the application of the known benchmark payloads. For all approaches, the identified masses, centre of gravities and inertias properly change with the applied load and show a relative error that is in large part significantly smaller than the maximum of 3.9%, measured for P z C at 100 kg with the Q5N method. Further, the constant parameters f q0 remain constant with a maximal deviation of 2.5% (f q02 at 100 kg Q5N). The relative errors increase with the applied load, where the average relative error per step is nearly constant, for example, ∆m P = −0.5%/20 kg for Q5M. The quasi-static methods show a slightly lower accuracy in the identification of the centre of gravity when compared to dynamic methods, as already expected from model analysis in Section 5. Least deviations in mass (<2%) and sensor offsets (<1.1%) shows D9N, where the other methods show slightly higher deviations that are similar between them. For the dynamic trajectories, the global standard deviationŝ 0 increases significantly with the load (factor 5.4 and 5.8 between 0 kg and 100 kg for D9M and D9N, resp.), whereas the quasi-static method is not affected that much (factor 1.8, Q5M) by the applied loads, even after reducing the number of poses (factor 1.7, Q5N). A reduction of eigenfrequencies with higher applied mass may increase the signal noise during movement and cause this effect. Again, the periodic trajectory D9N leads to higher uncertainties in identification, for example,ŝ 0 = 46.6 at 100 kg, caused by higher acceleration noise, when compared to the manually created trajectory D9M, for example,ŝ 0 = 29.8 at 100 kg.
More details are shown in Table 7. For the error correction parameters {k 0 , k 1 } theoretical values, and therefore deviations, are not obtainable. Differences between the methods partly result from different experimental setups, as static and dynamic methods were measured separately, where load positioning may vary slightly between the successive runs. N-optimal quasi-static results are extracted from the full quasi-static measurement by deleting 450 non-optimal poses.
It should be noted that a correct physical representation of the parameters is not the main focus of the identification procedure. In fact, all model parameters together must allow an accurate force measurement after calibration. Details on accuracy and precision in practise for quasi-static, dynamic and process force measurements with the end-effector setup are presented in [1]. All identification methods show acceptable accuracy and are suitable for the desired fast re-calibration of the measurement model. Hereby, dynamic measurement allows the identification of all parameters in a very fast measuring process (12 s), whereas the slightly slower quasi-static measurement (81 s) shows less noise at load conditions but does not include the inertia.

Summary
Structure-integrated force measurement is an accurate and cost-efficient alternative to classic force measurement in 6 DoF. To compensate for machine internal influences, such as dynamics, a measurement model is essential. Depending on the measuring task, classic solutions using FT sensors at the end-effector also require or benefit from a measurement model. However, every model needs a parametrisation that, in this case, both determines the measuring accuracy and, furthermore, can change during machine usage. Following modelling and model analysis, two methods are presented to identify these measurement model parameters, which run in a fast and easy procedure, in an operational machine state and do not require any additional knowledge or setup-efforts from the machine user. Where quasi-static identification measurement is based on data acquisition in stationary conditions at several poses, dynamic identification measurement is based on highly accelerated movements while taking samples continuously. Both methods are applied to the hexapod-based end-effector with integrated force sensors and are optimised to find the most significant measuring poses and the best most exciting trajectory, respectively. By the application and re-identification of reference payloads both methods are validated and compared. In conclusion, all of the presented procedures allow for an accurate determination of the model parameters with relative deviations to the theoretical values that are, in most cases, significantly smaller than the maximum error of −3.9% ( P z C Q5N, 100 kg), where the quasi-static method shows fewer uncertainties at load conditions (ŝ 0 = 2.3 Q5M vs.ŝ 0 = 46.6 D9N, 100 kg) and the dynamic method is faster (12 s D9N and 20 s D9M vs. 81 s Q5N and 21 min Q5M) and allows inertia identification.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Abbreviations
The following abbreviations and symbols are used in this manuscript: Coordinate vector/tensor: top left reference frame, bottom right body framê (·),( ·) Estimated values, and measured values (·), (·) Spatial vector (6 × 1), and vectorized matrix (·) −T Transposed inverse f q , f q0 , f K , f * q Measured forces, sensor bias, correction values, and processed sensor forces f Ext , f Dyn , f G External, Dynamic, and Gravitational spatial forces and moments J P , J Hex (x) Jacobian of the sensor framework and Jacobian of the hexapod kinematic m P , P x C , P I P Mass, centre of gravity, and inertia of the platform q,q,q Position, velocity, and acceleration of the drives x, B o P , B R P Pose, position, and orientation of the TCP v, Bȯ P , B ω P Velocity, translational velocity, and angular velocity of the TCṖ v, Bö P , Bω P Acceleration, translational acceleration, and angular acceleration of the TCP M P , C P , g P Mass matrix, matrix of centrifugal and Coriolis terms, gravitational terms A, u