A computational method for key-performance-indicator-based parameter identification of industrial manipulators

Abstract

We present a novel derivative-based parameter identification method to improve the precision at the tool center point of an industrial manipulator. The tool center point is directly considered in the optimization as part of the problem formulation as a key performance indicator. Additionally, our proposed method takes collision avoidance as special nonlinear constraints into account and is therefore suitable for industrial use. The performed numerical experiments show that the optimum experimental designs considering key performance indicators during optimization achieve a significant improvement in comparison to other methods. An improvement in terms of precision at the tool center point of 40% to 44% was achieved in experiments with three KUKA robots and 90 notional manipulator models compared to the heuristic experimental designs chosen by an experimenter as well as 10% to 19% compared to an existing state-of-the-art method.

1 Introduction

The total worldwide stock of operational industrial manipulators at the end of 2013 is estimated between 1,332,000 and 1,600,000 units [1]. Of all the different industries, the automotive industry requires the largest share of nearly 70,000 manipulators [1]. One application area of manipulators in automotive industry are flexible measurement systems (FMS) in assembly lines. During measurements with FMS, industrial manipulators collect measurement data for quality and process control, e.g. of automotive bodies. This working process requires a high tool center point (TCP) precision to detect production errors. The TCP defines the position and orientation of the working tool, which is attached to the last link of the manipulator.

1.1 Problem description

The TCP can be calculated using a geometric model of the manipulator by means of forward kinematics, e.g. the product of homogeneous transformations

$$ {}^{W}\mathrm {TCP}( \boldsymbol {p}, \boldsymbol {q}) = \prod _{i=1}^{k} {^{i-1}_{i}}\mathrm {T} ( \boldsymbol {p}_{i},q_{i}) = \prod_{i=1}^{k} \begin{bmatrix} {R}_{i}( \boldsymbol {p}_{i},q_{i}) & \boldsymbol {c}_{i}( \boldsymbol {p}_{i},q_{i}) \\ 0 & 1 \end{bmatrix} = \begin{bmatrix} {R}( \boldsymbol {p}, \boldsymbol {q}) & \boldsymbol {c}( \boldsymbol {p}, \boldsymbol {q}) \\ 0 & 1 \end{bmatrix} , $$

(1)

where ${^{i-1}_{i}}\mathrm {T} \in \mathbb {R}^{4 \times4}$ denotes the homogeneous transformation from the $(i-1)$th to the ith joint frame, where a frame is a Cartesian coordinate system attached to a joint. Each homogeneous transformation can be represented by means of an orthogonal matrix $R_{i} \in SO(3)$ describing the orientation and a vector $\boldsymbol {c}_{i}\in \mathbb {R}^{3}$ describing the position of the ith frame relative to the $(i-1)$th frame. The pose of the TCP is given by the matrix $R \in SO(3)$ describing the orientation and the vector $\boldsymbol {c}\in \mathbb {R}^{3}$ describing the position of the TCP in world coordinates, which is denoted by a 0 or W subscript if this is required for understanding.

The parameters $\boldsymbol {p}\in \mathbb {R}^{4\cdot k}$ describe the manipulator’s geometry such as the

• length of a link

• angle between two consecutive joint axes.

They are constant quantities, but due to manufacturing errors their exact values are unknown. In contrast, $q_{i} \in\mathbb{R}$ is a controllable value defining the pose of joint i. We denote by q either a single configuration of the manipulator, i.e. $\boldsymbol {q}\in \mathbb {R}^{k}$, or a set of n configurations, i.e. $\boldsymbol {q}\in \mathbb {R}^{n \cdot k}$, depending on the context. The homogeneous transformations ${^{i-1}_{i}}\mathrm {T}$, defining the position and orientation of link i with respect to link $(i-1)$, can be described in the Denavit-Hartenberg convention [2]

$$ {^{i-1}_{i}}\mathrm {T}_{\mathrm{DH}} ( \boldsymbol {p}_{i},q_{i}) = \mathrm {T}_{Rz}( \boldsymbol {p}_{i,0}+q_{i}) \mathrm {T}_{Tz}( \boldsymbol {p}_{i,1}) \mathrm {T}_{Tx}( \boldsymbol {p}_{i,2}) \mathrm {T}_{Rx}( \boldsymbol {p}_{i,3}) , $$

where

$$\begin{aligned}& \mathrm {T}_{Rz}(\theta_{i}) = \begin{bmatrix} \cos(\theta_{i}) & -\sin(\theta_{i}) & 0 & 0 \\ \sin(\theta_{i}) & \cos(\theta_{i}) & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} ,\qquad \mathrm {T}_{Tz}(d_{i}) = \begin{bmatrix} & 0 \\ \mathbb{I}_{3,3} & 0 \\ & d_{i} \\ 0 & 1 \end{bmatrix} , \\& \mathrm {T}_{Rx}(\alpha_{i}) = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & \cos(\alpha_{i}) & -\sin(\alpha_{i}) & 0 \\ 0 & \sin(\alpha_{i}) & \cos(\alpha_{i}) & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} ,\qquad \mathrm {T}_{Tx}(a_{i}) = \begin{bmatrix} & a_{i} \\ \mathbb{I}_{3,3} & 0 \\ & 0 \\ 0 & 1 \end{bmatrix} . \end{aligned}$$

$\mathrm {T}\in \mathbb {R}^{4}$ represents again a homogeneous transformation with the subscripts R, T specifying the type, i.e. rotation or translation, and the additional subscripts x, y, z denote the respective coordinate axis along or around the transformation takes place. Alternatively, the Hayati convention [3] can be used in the case two consecutive joints are parallel or almost parallel. Please note that the methodology derived in this article is not dependent on the chosen kinematic representation.

1.1.1 Parameter identification of industrial manipulators

The idea of parameter identification of manipulators is visualized in Figure 1. Parameter identification of manipulators and of more general industrial mechanisms, e.g. gantry manipulators or Gough-Stewart platforms, is required whenever differences between the model representing the ideal manipulator, e.g. (1), and the real manipulator exist. These differences are due to a mismatch between the parameters p of the model and the true parameters $\boldsymbol {p}^{\mathrm {true}}$ of the real manipulator causing a systematic error

$$ \mathrm {TCP}( \boldsymbol {p}, \boldsymbol {q}) - \mathrm {TCP} \bigl( \boldsymbol {p}^{\mathrm {true}},\boldsymbol {q}\bigr) \neq0 $$

(2)

between model-predicted and real TCP. Such deviations occur due to manufacturing tolerances, material failure, different environmental conditions like temperature and other effects. This mismatch can be reduced by exact identification of the mathematical model of the industrial manipulator, i.e. by adjusting the mathematical model to the real manipulator geometry. This is achieved by estimating the unknown parameters p of the model from measurement data, cf. [4]. Several measurement systems exist, which are used for the parameter identification procedure, e.g. laser trackers, laser modules, acoustic sensors, visual sensors, coordinate measuring machines and visual and automatic theodolites [5–7]. In this article, we only consider laser trackers as measurement devices in the numerical examples, however the actual methodology is applicable to any of the above mentioned measurement approaches.

Figure 1

Concept of parameter identification of industrial manipulators. Unknown geometric parameters p lead to deviations between model prediction and actual tool center point (TCP) pose. On the left, the actual TCP pose $\mathrm {TCP}( \boldsymbol {p}^{\mathrm {true}}, \boldsymbol {q})$ (in solid-black) and the predicted TCP pose $\mathrm {TCP}( \boldsymbol {p}, \boldsymbol {q})$ (in red-dashed) are shown before the parameter identification and they do not coincide. The identification tries to eliminate these deviations by estimating the true parameters of the real manipulator. On the right, the actual TCP pose (in solid-black) and the predicted TCP pose (in red-dashed) coincide after the parameter identification procedure (for $\boldsymbol {p}= \boldsymbol {p}^{\mathrm {true}}$).

The respective measurement system is mathematically described by the model response

$$ {h}: \mathbb {R}^{n_{p}} \times \mathbb {R}^{n_{q}} \to \mathbb {R}^{n_{h}}, $$

which describes the measurement device depending on the current geometric parameters and the measurement configuration determined by the configuration $\boldsymbol {q}^{\mathrm{meas}}$. Given a set of measurements $\boldsymbol {\eta}_{i} \in \mathbb {R}^{n_{h}}$, $i=1,\ldots,m$ from a measurement device, we further assume that parameters $\boldsymbol {p}^{\mathrm {true}}$ exist such that

$$ \boldsymbol {\eta}_{i} = {h}_{i} \bigl( \boldsymbol {p}^{\mathrm {true}}, \boldsymbol {q}_{i}^{\mathrm{meas}} \bigr) + \boldsymbol {\epsilon}_{i},\quad i=1, \ldots,m $$

with measurement errors

$$ \boldsymbol {\epsilon}_{i} \sim\mathcal{N}_{n_{h}} \bigl(0, { \Xi}_{i}^{2} \bigr),\quad i=1,\ldots ,m \text{ and } { \Xi}_{i} = \operatorname{diag}(\ldots, \sigma_{ij}, \ldots, j=1,\ldots,n_{h}). $$

Laser tracker measurement system

We briefly discuss a laser tracker measurement system as it is presented in Figure 2 to clarify the idea. In a laser tracker setup one or several reflectors are mounted on the end effector of the manipulator. A laser tracker then automatically tracks the reflector with its laser and measures the three dimensional position of the reflector relative to the coordinate system of the laser tracker. We assume that the origin 0 of the world frame is located in the base of the manipulator, that ${^{\mathrm{LT}}_{0}}\mathrm {T}$ is the homogeneous transformation from the laser tracker to the basis of the manipulator, that the manipulator has k joints and that the reflector position is ${^{k}}\boldsymbol {c}_{i}$ in the coordinate system of the TCP k. In several measurement configurations $\boldsymbol {q}_{i}^{\mathrm{meas}}$, $i=1,\ldots,m$ the laser tracker aims its laser at the reflector and measures the position of the reflector ${^{\mathrm{LT}}}\boldsymbol {c}_{i} \in R^{3}$ relative to the laser tracker coordinate system LT, e.g. by one of the following systems [5, 8, 9]. The model response relates the model position ${^{\mathrm{LT}}}\boldsymbol {c}_{i}$ of the reflector in the laser tracker coordinate system LT for the ith measurement configuration $\boldsymbol {q}_{i}^{\mathrm{meas}}$ and the TCP pose and is given by

$$ {h}\bigl( \boldsymbol {p}, \boldsymbol {q}_{i}^{\mathrm{meas}} \bigr) = {^{\mathrm{LT}}}\boldsymbol {c}_{i} = {^{\mathrm{LT}}_{0}} \mathrm {T} \cdot \mathrm {TCP} \bigl( \boldsymbol {p}, \boldsymbol {q}_{i}^{\mathrm{meas}} \bigr) \cdot{^{k}} \boldsymbol {c}_{i} \in\mathbb{R}^{3}. $$

(3)

Figure 2

A laser tracker measurement system for parameter identification of industrial manipulators. On the left, a laser tracker with a single reflector is presented. Here, the world frame is located in the foot of the manipulator. On the right, the figure shows the homogeneous transformations describing the manipulator and the laser tracker as well as the transition from the laser tracker coordinate system to the world coordinate system ${^{\mathrm{LT}}_{0}}\mathrm {T}$. Transition from joint $(i-1)$ to i is mathematically described by the homogeneous transformation ${^{i-1}_{i}}\mathrm {T}$. The reflector is mounted on the tool center point by the translation ${^{k}}\boldsymbol {c}_{j}$. Incident angle β between the reflector’s orientation vector ⁰ z in z-direction and the normalized laser beam ⁰ g from the reflector to the laser tracker is computed by $({^{0}}\boldsymbol {z}^{T} \cdot{^{0}}\boldsymbol {g}) \cdot\|{^{0}}\boldsymbol {g}\|$.

A measurement is considered feasible if the angle of incidence β between the reflector’s orientation vector ⁰ z in z-direction (also named central rotation axis) with length $\| {^{0}}\boldsymbol {z}\| =1$ and the laser beam ⁰ g from the reflector to the laser tracker with length $\|{^{0}}\boldsymbol {g}\| = 1$ is less than 30^∘, which can be described by the condition

$$ \cos \bigl(30^{\circ}\bigr) \le\cos(\beta) = \bigl({^{0}} \boldsymbol {z}^{T} \cdot{^{0}}\boldsymbol {g} \bigr) \cdot \bigl\Vert {^{0}} \boldsymbol {z} \bigr\Vert \cdot \bigl\Vert {^{0}}\boldsymbol {g} \bigr\Vert = \bigl({^{0}} \boldsymbol {z}^{T} \cdot{^{0}} \boldsymbol {g} \bigr) \cdot \bigl\Vert {^{0}} \boldsymbol {g} \bigr\Vert . $$

(4)

Parameter identification

The solution $\hat{ \boldsymbol {p}}\in \mathbb {R}^{n_{p}}$ of the nonlinear least-squares problem, given in the form of

$$ \min_{ \boldsymbol {p}} \sum_{i=1}^{m} \frac{1}{2} \bigl\Vert {\Xi_{i}}^{-1} \bigl( \boldsymbol {\eta }_{i} - {h}_{i} \bigl( \boldsymbol {p}, \boldsymbol {q}_{i}^{\mathrm{meas}} \bigr) \bigr) \bigr\Vert _{2}^{2} , $$

(5)

defines a local mapping $\hat{ \boldsymbol {p}} = g(\boldsymbol {\eta}, \boldsymbol {q}^{\mathrm{meas}})$ which yields an estimate for the true but unknown parameters $\boldsymbol {p}^{\mathrm {true}}$. Equation (5) can be written more compactly as

$$ \min_{ \boldsymbol {p}} \frac{1}{2} \bigl\Vert {f}_{1} ( \boldsymbol {p}, \boldsymbol {\eta}) \bigr\Vert _{2}^{2}, $$

(6)

where

$$\begin{aligned}& {f}_{1} := \Sigma^{-1}(\boldsymbol {\eta} - {h}) \in \mathbb {R}^{m\cdot n_{h}}, \\& \boldsymbol {\eta} :=\bigl[\boldsymbol {\eta}_{1}^{T}, \ldots, \boldsymbol { \eta}_{m}^{T} \bigr]^{T} \in \mathbb {R}^{m\cdot n_{h}}, \\& {h}:=\bigl[{h}_{1}^{T} \bigl( \boldsymbol {p}, \boldsymbol {q}_{1}^{\mathrm{meas}} \bigr), \ldots, {h}_{m}^{T} \bigl( \boldsymbol {p},\boldsymbol {q}_{m}^{\mathrm{meas}} \bigr) \bigr]b^{T} \in \mathbb {R}^{m\cdot n_{h}}, \\& {\Sigma} :=\operatorname{diag}({\Xi}_{i}, i =1,\ldots , m) \in \mathbb {R}^{m\cdot n_{h} \times m\cdot n_{h}} . \end{aligned}$$

Obviously, the systematic error (2) can be reduced by estimating $\boldsymbol {p}^{\mathrm {true}}$. Repeating the experiment leads to a different representation of measurements and thus to a different estimate $\hat{ \boldsymbol {p}}$. Since η is a random variable, $\hat{ \boldsymbol {p}}=g(\boldsymbol {\eta}, \boldsymbol {q})$ is also a random variable. One can define a confidence region

$$ CR(\hat{ \boldsymbol {p}},\alpha) = \bigl\{ \boldsymbol {p}: ( \boldsymbol {p}-\hat{ \boldsymbol {p}})^{T} {C}( \boldsymbol {p}-\hat{ \boldsymbol {p}}) \le\chi _{n_{p}}^{2} (1-\alpha) \bigr\} $$

containing the true parameter values to a certain probability $1-\alpha$. $\chi_{n_{p}}^{2}$ is the $\chi^{2}$ distribution with $n_{p}$ degrees of freedom. C is the variance-covariance matrix defined by

$$ {C}= \bigl({ {F}}_{1}^{T} { {F}}_{1} \bigr)^{-1} \in \mathbb {R}^{n_{p} \times n_{p}} \quad \text{with } { {F}}_{1}:=\frac{\mathrm{d} {f}_{1}}{\mathrm{d} \boldsymbol {p}} \bigl( \boldsymbol {p}, \boldsymbol {q}^{\mathrm{meas}} \bigr). $$

Please note the explicit dependence of the parameter covariance matrix on the measurement configurations $\boldsymbol {q}^{\mathrm{meas}}$, we will write ${C}(\boldsymbol {q}^{\mathrm{meas}})$ to emphasize this dependence when required.

1.2 Previous work on optimum experimental design

State-of-the-art methods for parameter identification of industrial manipulators estimate the geometric parameters $\boldsymbol {p}\in\mathbb {R}^{n_{p}}$ by solving the nonlinear least-squares problem (6). Fitting the model to the measurement data by means of parameter estimation does not imply the minimization of the respective uncertainties of the resulting parameters given by the variance-covariance matrix ${C}( \boldsymbol {q}^{\mathrm{meas}})$. The quality of the parameter estimates depends on the choice of measurement configurations $\boldsymbol {q}^{\mathrm{meas}}$. Therefore it is desired to use configurations, which provide maximal information gain and low statistical uncertainty. In that view, the measurement configurations $\boldsymbol {q}^{\mathrm{meas}}$ can be optimized by solving the optimum experimental design (OED) problem

$$\begin{aligned}& \min_{ \boldsymbol {q}^{\mathrm{meas}}} \phi\bigl( {C}\bigl( \boldsymbol {p}, \boldsymbol {q}^{\mathrm{meas}}\bigr)\bigr) \end{aligned}$$

(7a)

$$\begin{aligned}& \quad \text{s.t. } 0 \le\psi \bigl( \boldsymbol {p}, \boldsymbol {q}^{\mathrm{meas}} \bigr). \end{aligned}$$

(7b)

A cost function $\phi( {C}( \boldsymbol {p}, \boldsymbol {q}^{\mathrm{meas}}))$ of the covariance matrix is minimized subject to kinematic constraints $0 \le\psi( \boldsymbol {p}, \boldsymbol {q}^{\mathrm{meas}})$, which comprise

• manipulator specifications,

• measurement system specifications and

• collision avoidance.

The solution of the optimization problem provides measurement configurations with which a minimal uncertainty during parameter estimation is achieved. We summarize the state-of-the-art approach and give a comparative overview with respect to the methodology presented in this article in Table 1.

Table 1 An overview of relevant methods from the literature applied by other authors

In [10] and [11] the quality of a given set of measurement configurations is defined by five so-called observability indices (OIs) which were first proposed by Borm and Menq in [12] and [13]. They are based on the sensitivities w.r.t. geometric parameters and configurations ${F}_{1}:=\frac{\mathrm{d} {f}_{1}}{\mathrm{d} \boldsymbol {p}}( \boldsymbol {p},\boldsymbol {q}^{\mathrm{meas}})$ from the parameter estimation problem (6). In [10] the five quality criteria for parameter identification of industrial manipulators are related to those of OED problems [14]. Each of the quality criteria can be used as a cost function in the optimization problem (7a)-(7b).

The optimization of measurement configurations with respect to an OI achieving minimal parameter errors during parameter identification of industrial manipulators is an active field of research. Many solutions have been proposed which are based on derivative-free and gradient-based methods. Borm and Menq [12, 13] present the optimization of measurement configurations with respect to the objective

$$ \max_{ \boldsymbol {q}^{\mathrm{meas}}} \phi _{01} = \max_{ \boldsymbol {q}^{\mathrm{meas}}} \frac{\sqrt [r]{\mu_{1} \cdot\ldots \cdot\mu_{r}}}{\sqrt{r}}, $$

(8)

where $\mu_{1}, \ldots, \mu_{r}$ are the singular values of ${F}_{1}$ in descending order. They solve the optimization problem by a steepest descent and a successive quadratic line search method. In the optimization no mechanical joint constraints are considered. In [15] a conjugate gradient type method is used to maximize the product of singular values. In [16] and [17] the optimization problem

$$ \min_{ \boldsymbol {q}^{\mathrm{meas}}} \phi _{02} = \min_{ \boldsymbol {q}^{\mathrm{meas}}} \frac{\mu _{1}}{\mu_{r}} $$

(9)

is solved by a genetic algorithm and simulated annealing. During the optimization joint constraints are considered. The same objective is used in [18] to compute optimized TCP poses in a constraint search space with the function fmincon from MATLAB (The Mathworks Inc.).

Daney studied in [19] the objective

$$ \min_{ \boldsymbol {q}^{\mathrm{meas}}} \phi _{03} = \max_{ \boldsymbol {q}^{\mathrm{meas}}} \operatorname {det} \bigl( {F}_{1}^{T} {F}_{1} \bigr) $$

(10)

by using tabu search. Restrictions to the sliding joints of the manipulator are considered during computation. In [20] a further optimality criterion for OED is defined called $D^{*}$-criterion, which is adapted to the special structure of the information matrix. With regard to this criterion, optimal measurement configurations are selected from a set of configurations. In [21] the DETMAX algorithm is used to select optimal measurement configurations from a set of given configurations due to the objective $\max_{\boldsymbol {q}^{\mathrm{meas}}} \phi _{01}$. The same algorithm is used in [22] selecting optimized measurement configurations with respect to the five existing OI.

The approach of nonlinear experimental design for explicit key performance indicators (KPIs) applied to a chemical process has been discussed in [23].

1.3 Contribution of the article

All presented formulations for the objective function in (7a)-(7b) only consider the statistical uncertainties of the geometric parameters of the kinematic model. We present a new optimization problem formulation in which additionally the TCP position or other KPIs are considered in the objective. Thereby, the optimization problem (7a)-(7b) is extended with the statistical error of TCP positions given by a number of pre-defined working configurations. The OED problem is solved by a sequential quadratic programming (SQP) method. The avoidance of collisions has to be considered in the optimization method to be applicable in industrial practice. This is achieved by introducing additional nonlinear constraints to the optimization problem. The segments of the manipulator are approximated by capsules. The avoidance of collisions between two manipulator segments is described by a single nonlinear constraint in the optimization problem.

1.4 Organization of the article

In Section 2, we explain the methodology of OED for parameter identification of industrial manipulators in detail. First, we introduce the general formulation of KPIs. Second, we formulate the OED problem for KPIs. Afterwards, we introduce a collision avoidance strategy to generate collision-free measurement configurations. In Section 3, we present numerical results from two simulation studies. First, the KPI approach is presented for a simulation of real-world flexible measurement systems (FMS) example with three different KUKA robots (KR15, KR300 and KR500) considering collision avoidance. Second, for a set of 90 different manipulator geometries the proposed methodology is statistically compared against a heuristic and a state-of-the-art method in simulation. Afterwards, the results are discussed. We finish with a conclusion and give an outlook to future work.

2 Parameter identification of manipulators with key performance indicators

In the previous section, we introduced the principle method for parameter identification of manipulators, i.e. parameter estimation using nonlinear least-squares and a kinematic model of the manipulator geometry. Furthermore, we presented linear error propagation and the quantification of the uncertainties of the geometric parameters p by means of the variance-covariance matrix C, which takes the error propagation from measurements to parameters into account, i.e.

$$ \boldsymbol {\eta} \mapsto\hat{ \boldsymbol {p}} = g \bigl(\boldsymbol {\eta}, \boldsymbol {q}^{\mathrm{meas}} \bigr) \quad \text{and}\quad \bigl(\hat{ \boldsymbol {p}}, \boldsymbol {q}^{\mathrm{meas}} \bigr) \mapsto\phi \bigl(C(\hat{ \boldsymbol {p}}) \bigr). $$

(11)

In Section 1.2, we presented the state of the art of optimum experimental design (OED) for parameter identification of manipulators based on the minimization of so-called observability indices (OIs). This approach yields optimal measurement configurations $\boldsymbol {q}^{\mathrm{meas}}$ maximizing the OI and with this the statistical uncertainty of the parameter estimates $\hat{ \boldsymbol {p}}$ given by ${C}(\hat{ \boldsymbol {p}}, \boldsymbol {q}^{\mathrm{meas}})$. However, most often we are not interested in an estimate of the geometric parameters of the manipulator with a low statistical error but in having a high tool center point (TCP) precision.

In the following, we will present a mathematical concept in which the TCP precision as the main quantity of interest is explicitly considered. In the remainder of this article we will denote this main quantity of interest as key performance indicator (KPI). Primarily, we are interested in the TCP precision and, secondarily, in the precision of the geometric parameters.

In this context we use the example of flexible measurement systems (FMS) to derive the methodology of optimum experimental design for key performance indicators. In FMS an industrial manipulator measures automotive bodies, requiring a high TCP precision. An example setup is shown in Figure 3.

Figure 3

An exemplary setup of a $\pmb{360^{\circ}}$ SIMS flexible measurement system. Image courtesy of Hexagon Manufacturing Intelligence.

The manipulator takes measurements in the positions determined by the working configurations $\boldsymbol {q}_{i}^{\mathrm{work}} \in\mathbb{R}^{k}$ $\forall i=1,\ldots,t$ in which a high TCP precision is desired. These TCP positions are introduced as additional $\boldsymbol {s}\in \mathbb {R}^{n_{s}}$ variables

$$ \boldsymbol {s}= \mathrm {TCP} \bigl( \boldsymbol {p}, \boldsymbol {q}^{\mathrm{work}} \bigr) $$

(12)

depending implicitly on the parameters p and the working configurations $\boldsymbol {q}^{\mathrm{work}}$. Formulation (12) can be used to quantify the dependence of the KPI on parametric uncertainty by means of the explicit mapping $\boldsymbol {\eta} \mapsto \boldsymbol {p}\mapsto \boldsymbol {s}$ or the implicit mapping $\boldsymbol {\eta} \mapsto( \boldsymbol {p}, \boldsymbol {s}( \boldsymbol {p}))$.

2.1 Formulation of key performance indicators

In the following, we derive how the error propagation from the measurements to parameters as described by Equation (11) can be extended to specific quantities of interest, the so-called key performance indicators (KPIs). For this, the error propagation from the data η to the KPIs s is quantified by $\phi(C( \boldsymbol {s}))$.

For given configurations $\boldsymbol {q}=( \boldsymbol {q}^{\mathrm{meas}}, \boldsymbol {q}^{\mathrm{work}})^{T}$ and related data η, the least-squares problem (6) is augmented by nonlinear constraints ${f}_{2} \in \mathbb {R}^{n_{2}}$ to

$$\begin{aligned}& \min_{ \boldsymbol {v}} \frac{1}{2} \bigl\Vert {f}_{1} \bigl( \boldsymbol {v}, \boldsymbol {q}^{\mathrm{meas}}, \boldsymbol {\eta } \bigr) \bigr\Vert _{2}^{2} \end{aligned}$$

(13a)

$$\begin{aligned}& \quad \text{s.t. } 0 = {f}_{2} \bigl( \boldsymbol {v}, \boldsymbol {q}^{\mathrm{work}} \bigr) \end{aligned}$$

(13b)

with variables $\boldsymbol {v}=( \boldsymbol {p}, \boldsymbol {s})^{T}$, where ${f}_{1} \in \mathbb {R}^{m}$ is the least-squares objective of the parameter estimation problem (13a)-(13b) and ${f}_{2}$ represents the relation between parameters p, configurations $\boldsymbol {q}^{\mathrm{work}}$ and KPIs s. Specific formulations are for example the investigated precision of static TCP poses but also the precision of TCP path trajectories could be described. This general formulation additionally includes explicit KPI formulations like Equation (12), i.e. ${f}_{2} = \{ \mathrm {TCP}( \boldsymbol {p}, \boldsymbol {q}_{i}^{\mathrm{work}}) - \boldsymbol {s}_{i} \}_{i=1}^{t}$. In the Appendix, the transformation of the constrained least-squares problem (13a)-(13b) to an equivalent unconstrained problem under certain assumptions is derived.

Analogously to (6), the solution $\hat{ \boldsymbol {v}} = (\hat{ \boldsymbol {p}}, \hat{ \boldsymbol {s}})^{T}$ of the constrained nonlinear least-squares problem (13a)-(13b) is subject to statistical uncertainty due to the stochasticity of the measurement data η. The propagated uncertainty for both, the parameter and the KPI estimates $\hat{ \boldsymbol {p}}$, $\hat{ \boldsymbol {s}}$, can be considered for the minimization in OED.

In the following, we give the required definitions and assumptions to formulate the OED problem for KPIs. Assume in the following discussion that

$$ {F}_{1} := \frac{\mathrm{d}}{\mathrm{d} \boldsymbol {v}}{f}_{1}( \boldsymbol {v}, \boldsymbol {q}, \boldsymbol {\eta}) \in \mathbb {R}^{m \times n_{v}} \quad \text{and}\quad {F}_{2} := \frac{\mathrm{d}}{\mathrm{d} \boldsymbol {v}} {f}_{2}( \boldsymbol {v}, \boldsymbol {q}) \in \mathbb {R}^{n_{2} \times n_{v}}, $$

(14)

satisfy the regularity conditions

$$ (\mathrm{CQ}) \quad \operatorname{rank}( {F}_{2}) = n_{2} \quad \text{and}\quad (\mathrm{PD}) \quad \operatorname{rank}( {J}) = n_{v}, $$

(15)

${J}= ( {F}_{1}^{T}, {F}_{2}^{T})^{T} \in \mathbb {R}^{(m + n_{2}) \times n_{v}}$.

As for the unconstrained least-squares problem, one can describe a confidence region of the solution $\hat{ \boldsymbol {v}}$. The statistical uncertainty of the measurements $\boldsymbol {\eta} \in \mathbb {R}^{m}$ is propagated to the estimated parameters. We denote the mapping from measurement data η to estimated quantities $\hat{ \boldsymbol {v}}$, including the KPIs, by

$$\begin{aligned} \begin{aligned} g_{v}:{}& \mathbb {R}^{m} \to \mathbb {R}^{n_{v}} \\ &\boldsymbol {\eta} \mapsto\hat{ \boldsymbol {v}} = g_{v}(\boldsymbol {\eta}; \boldsymbol {q}) , \end{aligned} \end{aligned}$$

(16)

the solution of the constrained least-squares problem (13a)-(13b) with $\boldsymbol {q}= ( { \boldsymbol {q}^{\mathrm{meas}}}, { \boldsymbol {q}^{\mathrm{work}}} )$. It is well-known that for a normally distributed random variable $X \sim\mathcal {N}(\mu, \Sigma^{2})$ an affine mapping $Y=g(X) = AX + b$ yields another normally distributed random variable $Y \sim\mathcal {N}(A\mu+ b, A \Sigma^{2} A^{T})$. For nonlinear functions one can linearize the function and use the linear error propagation approximation. Thus, for $g_{v}$ the random variable $V \sim\mathcal {N} ( g_{v}(\boldsymbol {\eta}; \boldsymbol {q}), \frac{\mathrm{d} g_{v}}{\mathrm{d} \boldsymbol {\eta}}(\boldsymbol {\eta}; \boldsymbol {q}) \Sigma^{2} \frac{\mathrm{d} g_{v}}{\mathrm{d} \boldsymbol {\eta}}(\boldsymbol {\eta}; \boldsymbol {q})^{T} )$ is obtained.

Lemma

Under the assumptions (CQ), (PD), it holds that

$$\begin{aligned} \Delta \boldsymbol {v}&= \frac{\mathrm{d} g_{v}}{\mathrm{d} \boldsymbol {\eta}}(\boldsymbol {\eta}; \boldsymbol {q}) \Delta \boldsymbol {\eta} \end{aligned}$$

(17)

$$\begin{aligned} & = \bigl[ \textstyle\begin{array}{@{}c@{\quad}c@{}} \mathbb{I} & 0 \end{array}\displaystyle \bigr] \begin{bmatrix} {F}_{1}^{T} {F}_{1} & {F}_{2}^{T} \\ {F}_{2} & 0 \end{bmatrix} ^{-1} \begin{bmatrix} {F}_{1}^{T} \Sigma^{-1} \\ 0 \end{bmatrix} \Delta \boldsymbol {\eta}, \end{aligned}$$

(18)

where $\mathbb{I}\in \mathbb {R}^{n_{v} \times n_{v}}$, ${F}_{1} \equiv {F}_{1}(\boldsymbol {\eta}, \boldsymbol {q})$, ${F}_{2} \equiv {F}_{2}( \boldsymbol {q})$ and $\Delta \boldsymbol {v}:= \boldsymbol {v}- \hat{ \boldsymbol {v}}$, $\Delta \boldsymbol {\eta} :=\boldsymbol {\eta} - \hat{\boldsymbol {\eta}}$ denote the respective differences between current and expected values.

Proof

See [24]. □

In the following, $\mathbb {E}(\cdot)$ denotes the expected value of a random variable. From this Lemma it follows that

$$\begin{aligned} {C}_{v} & = \mathbb {E}\bigl( \Delta \boldsymbol {v}\Delta \boldsymbol {v}^{T} \bigr) \\ & = \bigl[ \textstyle\begin{array}{@{}c@{\quad}c@{}} \mathbb{I} & 0 \end{array}\displaystyle \bigr] \begin{bmatrix} {F}_{1}^{T} {F}_{1} & {F}_{2}^{T} \\ {F}_{2} & 0 \end{bmatrix} ^{-1} \begin{bmatrix} {F}_{1}^{T} {F}_{1} & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} {F}_{1}^{T} {F}_{1} & {F}_{2}^{T} \\ {F}_{2} & 0 \end{bmatrix} ^{-T} \begin{bmatrix} \mathbb{I} \\ 0 \end{bmatrix} . \end{aligned}$$

The formula can be shown to be equivalent to Equation (19):

Proposition

The variance-covariance matrix ${C}_{v} :=\mathbb {E}( ( \boldsymbol {v}- \mathbb {E} \boldsymbol {v}) ( \boldsymbol {v}- \mathbb {E} \boldsymbol {v})^{T} )$ satisfies

$$ {C}_{v} = \begin{bmatrix} {C}_{p} & {C}_{ps} \\ {C}_{sp} & {C}_{s} \end{bmatrix} = \big[ \textstyle\begin{array}{@{}c@{\quad}c@{}} \mathbb{I} & 0 \end{array}\displaystyle \bigr] \begin{bmatrix} {F}_{1}^{T} {F}_{1} & {F}_{2}^{T} \\ {F}_{2} & 0 \end{bmatrix} ^{-1} \begin{bmatrix} \mathbb{I} \\ 0 \end{bmatrix} . $$

(19)

Proof

See [24]. □

2.2 Collision avoidance

In the following, we present a collision avoidance strategy required for the practical realization of OED for KPIs in industrial applications. The optimized measurement configurations $\boldsymbol {q}^{\mathrm{meas}}$ are required to be free of collisions between the manipulator and its environment. This is achieved by introducing a collision avoidance strategy to the formulation of the OED problem.

We give an overview of existing collision avoidance approaches. Afterwards, our approach is presented in a general manner for the fact that the collision avoidance formulation is not restricted to flexible measurement systems. In Section 3.2 the collision avoidance technique is then applied on a specific experimental setting.

Collision avoidance is a major research field in which different approaches have already been successfully applied. In [25] the manipulator segments and obstacles are approximated by the union of convex polyhedra and the collision avoidance criterion between two polyhedra is based on Farkas’ Lemma. In [26] the body of a humanoid manipulator is described by ‘strictly convex hulls’. These hulls result from slightly blowing up the usual convex hulls through patches of spheres and tori. As a consequence, the gradient of the proximity distance function becomes continuous, which is useful in a derivative-based optimization framework.

The shapes of industrial manipulators are less complex than those of humanoid robots. Therefore, we follow the approach from [27], where the industrial manipulator as well as the obstacles in the working environment are approximated by spheres and capsules. The approximation process is exemplarily visualized in Figures 4(a) and 4(b). Hereby, each link c of the manipulator is described by appropriate capsules $b_{i} = \operatorname{Capsule}({^{W}_{c}}b_{i}^{A}, {^{W}_{c}}b_{i}^{B}, r_{i})$, $i \in\mathcal {R}$, where $\mathcal{R}$ is the index set of all manipulator specific capsules. Each capsule consists of a radius $r_{i}$ and its end points ${^{W}_{c}} b_{i}^{A}, {^{W}_{c}}b_{i}^{B} \in \mathbb {R}^{3}$ for the manipulator segment c in the world coordinate system W given by

$$ {^{W}_{c}} b_{k}^{A} = \prod _{i=1}^{c} {^{i-1}_{i}}\mathrm {T}( \boldsymbol {p}_{i}, \boldsymbol {q}_{i}) \cdot{^{c}}b_{k}^{A} \quad \text{and}\quad {^{W}_{c}} b_{k}^{B} = \prod_{i=1}^{c} {^{i-1}_{i}} \mathrm {T}( \boldsymbol {p}_{i}, \boldsymbol {q}_{i}) \cdot {^{c}}b_{k}^{B}. $$

Similarly, each object is approximated by a capsule $b_{j} = \operatorname{Capsule}({^{W}}b_{j}^{A},{^{W}} b_{j}^{B}, r_{j})$, $j \in\mathcal{O}$ with the index set of all environment specific capsules $\mathcal{O}$.

Figure 4

Visualization of the capsule approximation. Each joint $c_{i}$, $i \in\{0,\ldots,4\}$ in column (a) is described by an appropriate capsule $b_{i}$, $i \in\{0,\ldots,4\}$ illustrated in column (b). In column (c) the distance $d_{P}$ between the two line segments $b_{1}$ and $b_{3}$ is shown.

A configuration is collision-free when the distance between each pair of capsules is larger than zero. The distance between two capsules $b_{1}$ and $b_{3}$ of two different joints $c_{1}$ and $c_{3}$ can be computed by

$$ d( b_{1}, b_{3} ) = d_{P} ( b_{1}, b_{3} ) - ( r_{1} + r_{3} ) \ge0 $$

(20)

with

$$ d_{P} ( b_{1}, b_{3} ) = \sqrt {\sum_{i=1}^{3} \bigl[ \bigl({^{W}_{c1} } b_{1i}^{B} - {^{W}_{c1} } b_{1i}^{A} \bigr)t - \bigl({^{W}_{c3} } b_{3i}^{B} - {^{W}_{c3} } b_{3i}^{A} \bigr)u - {^{W}_{c3} } b_{3i}^{A} + {^{W}_{c1} } b_{1i}^{A} \bigr]^{2} } , \quad t,u \in\mathbb{R}, $$

(21)

which computes the distance of the two respective line segments, cf. Figure 4(c). Lumelsky introduced an algorithm to compute the distance $d_{P}$ between two line segments in [28]. The three dimensional coordinates of the endpoints of line segments serve as input values for the algorithm. First, three special cases are considered preventing computational difficulties due to small nominators and analyzing if one or both segments degenerate into points or if the line segments are parallel. If none of the cases holds, the minimal distance is computed using a parameterized formulation of the line segments. Following this, a nonlinear constraint like (20) for each pair of considered capsules is formulated in the optimization problem which guarantees optimized, collision-free configurations.

If the manipulator consists of n segments each approximated by a single capsule, the maximal number of constraints for collision avoidance can be computed by the formula

$$ \text{max \# of constraints} = \frac{n (n -1)}{2} - (n-1). $$

(22)

Two neighboring segments do not have to be compared as the limitation of joint movement automatically prevents their collision.

Not all segments can be described by just one suited capsule. If the first or last joint is described by k capsules, then additionally $(k-1)\cdot(n-2)$ constraints have to be added to (22). If one of the interior segments is described by k capsules this results in $(k-1)\cdot(n-3)$ additional constraints. The total number of constraints can be decreased by further investigation of the geometry of the manipulator and motion to find segments which will never collide and therefore do not have to be considered in the collision avoidance strategy. A study in [29] has shown, that the blow up from inlying to enclosing capsules only has a slightly negative influence on the optimization performance. As the description of an industrial manipulator by capsules is just an approximation of reality, the use of enclosing capsules guarantees by higher probability that the optimized measurement configurations are collision-free in a practical application without losing much performance.

Analyzing the derivative of (20) w.r.t. the control variables q yields that the first fraction has a singularity if the distance between two line segments becomes zero. But this will never be the case as in each constraint the distance of line segments must be equal or larger than the sum of the capsules’ radii.

$$\begin{aligned} \frac{\mathrm{d} d( b_{1}, b_{2} )}{\mathrm{d} \boldsymbol {q}} =& \frac{1}{2 \sqrt{\sum_{i=1}^{3} [ ({^{W}_{c}} b_{i1}^{B} - {^{W}_{c}} b_{i1}^{A} )t - ({^{W}_{c}} b_{i2}^{B} - {^{W}_{c}} b_{i2}^{A})u - {^{W}_{c}} b_{i2}^{A} + {^{W}_{c}} b_{i1}^{A} ]^{2} }} \\ &{} \cdot\frac{1}{\frac{\partial}{\partial \boldsymbol {q}}\sum_{i=1}^{3} [ ({^{W}_{c}} b_{i1}^{B} - {^{W}_{c}} b_{i1}^{A} )t - ({^{W}_{c}} b_{i2}^{B} - {^{W}_{c}} b_{i2}^{A})u - {^{W}_{c}} b_{i2}^{A} + {^{W}_{c}} b_{i1}^{A} ]^{2}} \\ &{}\cdot\frac{1}{\frac{\partial}{\partial \boldsymbol {q}} ({^{W}_{c}} b_{i1}^{B} - {^{W}_{c}} b_{i1}^{A} )t - ({^{W}_{c}} b_{i2}^{B} - {^{W}_{c}} b_{i2}^{A})u - {^{W}_{c}} b_{i2}^{A} + {^{W}_{c}} b_{i1}^{A}}. \end{aligned}$$

(23)

2.3 Optimum experimental design for key performance indicators

We now assemble the problem formulation of OED for KPIs considered in this article. A function of the variance-covariance matrix ${C}( \boldsymbol {q}^{\mathrm{meas}})$ is minimized to find optimal measurement configurations for the parameter identification procedure. Depending on the choice of the objective, either the parametric uncertainty of the geometric parameters and/or the KPIs is minimized. The explicit dependence of ${C}( \boldsymbol {q}^{\mathrm{meas}})$ can be exploited. Collision avoidance is considered via constraints. This results in the following OED problem

$$\begin{aligned}& \min_{ \boldsymbol {q}^{\mathrm{meas}}} \operatorname{tr} \bigl(P^{T} {C}_{v}( \boldsymbol {q}) P \bigr) \end{aligned}$$

(24a)

$$\begin{aligned}& \quad \text{s.t. } 0 \le\psi \bigl( \boldsymbol {p}, \boldsymbol {q}^{\mathrm{meas}} \bigr) \textstyle\begin{cases} q_{i}^{\mathrm{lo}} \le q_{i}^{\mathrm{meas}} \le q_{i}^{\mathrm{up}}, & i = 1,\ldots,k, \\ ({^{0}}\boldsymbol {z}^{T} \cdot{^{0}}\boldsymbol {g}) \cdot\|{^{0}}\boldsymbol {g}\| \le\cos (30^{\circ}), \\ d (b_{i},b_{j}) \ge0, & i,j \in\mathcal{R}, \{i,j\} \notin\mathcal{I}, \\ d (b_{i},b_{j}) \ge0, & i \in\mathcal{R}, j \in\mathcal{O} \end{cases}\displaystyle \end{aligned}$$

(24b)

$$\begin{aligned}& \quad \text{with } {C}_{v} = \bigl[ \textstyle\begin{array}{@{}c@{\quad}c@{}} \mathbb{I} & 0 \end{array}\displaystyle \bigr] \begin{bmatrix} {F}_{1}^{T} {F}_{1} & {F}_{2}^{T} \\ {F}_{2} & 0 \end{bmatrix} ^{-1} \begin{bmatrix} \mathbb{I} \\ 0 \end{bmatrix} , \end{aligned}$$

(24c)

$$\begin{aligned}& \hphantom{\quad \text{with }} {F}_{1} = \Sigma^{-1} \frac{\mathrm{d}}{\mathrm{d} \boldsymbol {v}} {h}\bigl( \boldsymbol {p}, \boldsymbol {q}^{\mathrm{meas}} \bigr)\quad \text{and} \end{aligned}$$

(24d)

$$\begin{aligned}& \hphantom{\quad \text{with }} {F}_{2} = \frac{\mathrm{d}}{\mathrm{d} \boldsymbol {v}} \mathrm {TCP} \bigl( \boldsymbol {p}, \boldsymbol {q}^{\mathrm{work}} \bigr). \end{aligned}$$

(24e)

P is a matrix projecting onto a subset of p, of s or of p and s with the possibility to scale each parameter independently. $\mathcal{I}$ is a set of manipulator segment indices, which do not have to be compared for collision avoidance. The variance-covariance matrix ${C}_{v}$ results from the constrained parameter estimation problem (13a)-(13b) with derivatives $F_{1}$ and $F_{2}$ as introduced in Section 2.1. As objective function (24a) we choose the arithmetic mean of the diagonal element of the projected variance-covariance matrix, the so-called A criterion, corresponding to the average of the half-axes of the confidence ellipsoid, see [30]. The constraints $\psi( \boldsymbol {p}, \boldsymbol {q}^{\mathrm{meas}})$ include restrictions for the movement of each joint i, the incident angle between laser beam and the orientation vector of the TCP in z direction and the collision avoidance between manipulator segments and between manipulator segments and obstacles. In problem (24a)-(24e) the control variables $\boldsymbol {q}^{\mathrm{work}}$ are fixed defining the TCP pose during the work process and only the measurement configurations $\boldsymbol {q}^{\mathrm{meas}}$ are optimization variables. As problem (24a)-(24e) does not depend on the data η, optimal experimental designs can be computed before the experiments are carried out.

3 Numerical results

The task of a manipulator in flexible measurement systems is to detect small production errors. Hence a high precision at the tool center point (TCP) of the manipulator is required. The precision is improved by sequentially identifying the manipulator’s parameters in a fixed time interval, e.g. after each work cycle of the manipulator. The work cycle consists of one working configuration $\boldsymbol {q}^{\mathrm{work}} \in \mathbb{R}^{6}$ after which the re-identification of the parameters is performed with 60 measurement configurations. The TCP of the working configuration is considered as key performance indicator (KPI) inside the optimum experimental design (OED) formulation as the manipulator requires a high TCP precision during its work cycle. A general setup for flexible measurement systems is defined in Section 3.2 and is used as an in-silico experiment to introduce and discuss different optimum experimental designs that provide the best improvement of the TCP precision in comparison to the configurations chosen by a heuristic or randomly. The designs differ in the consideration of collision avoidance and in different choices of variance-covariance matrix projections

$$ {C}_{p} = \bigl[ \textstyle\begin{array}{@{}c@{\quad}c@{}} \mathbb{I}_{1} & 0 \end{array}\displaystyle \bigr] {C}_{v} \begin{bmatrix} \mathbb{I}_{1} \\ 0 \end{bmatrix} \quad \text{and}\quad {C}_{s} = \bigl[ \textstyle\begin{array}{@{}c@{\quad}c@{}} 0 & \mathbb{I}_{2} \end{array}\displaystyle \bigr] {C}_{v} \begin{bmatrix} 0 \\ \mathbb{I}_{2} \end{bmatrix} , $$

(25)

with $\mathbb{I}_{1} \in\mathbb{R}^{n_{p} \times n_{p}}$ and $\mathbb{I}_{2} \in\mathbb{R}^{n_{s} \times n_{s}}$ considered in the cost function of the optimization problem. Covariance matrix ${C}_{p}$ only considers the statistical uncertainties of the geometric parameters, while ${C}_{s}$ considers the uncertainties of the KPI only. The OED problems are solved with the parameter values presented in Table 2 and Table 3. Please note that no parameter estimations are performed. In both tables, we highlight the parameters whose uncertainties are considered in the optimization. The different tasks in the in-silico experiment are depicted by Figure 5, which gives an overview of the work flow of the parameter identification process. In the following all tasks are performed to compute optimal sets of configurations for the current parameter identification.

Figure 5

Overview of the parameter identification procedure under consideration. The figure shows the general approach to parameter identification of industrial manipulators and mechanisms using KPIs. We assume the manipulator has to perform a measurement task on certain TCP positions s requiring a high precision. These positions are considered as KPIs for the parameter identification procedure. The required precision is achieved by combining a (re-)identification of the parameters with the actual work cycle of the manipulator, i.e. whenever the manipulator is idle it can sequentially identify itself. Both the measurement and parameter identification task take place in the same scene. From scene and manipulator model the necessary obstacle descriptions are required to setup the collision avoidance for the OED problem. The OED problem will find an optimal set of configurations q for the optimization problem considering the collision avoidance and the chosen KPI from a given set of initial configurations. Whenever the manipulator performs an identification task, the internal kinematic model is updated with the current optimal parameter values $\boldsymbol {p}^{\mathrm {true}}$.

Table 2 Denavit-Hartenberg and Hayati zero offset values with $\pmb{\lambda\in\{ 650,700,\ldots, 1\text{,}050 \} }$ and $\pmb{\mu\in\{ 600, 650,\ldots , 1\text{,}050 \} }$

Table 3 Position and orientation of laser tracker and reflector

3.1 Software

The optimization problem is solved with the software package VPLAN [30], developed in the research group of the authors. The evaluation of derivatives is performed by ADIFOR [31]. ADIFOR is able to evaluate the derivatives of the vector-valued objective function (6). The OED problem (24a)-(24e) is formulated as a nonlinear optimization problem and solved by the sparse SQP solver SNOPT7 [32].

3.2 Experimental setup

The experimental setup of the flexible measurement system example is visualized in Figure 6. The industrial manipulator is mounted on a long object next to a wall and the laser tracker for parameter identification is located in the left corner in the line of sight of the manipulator. Directly in front of the manipulator a car is positioned, whose body has to be measured and which acts as an obstacle for the choice of possible measurement configurations. The computation of optimum experimental designs is divided into two experiments. Firstly, we study the performance of our approach on three different types of KUKA robots (KR15, KR300, KR500). Afterwards we quantify our approach by using 90 notional industrial manipulators. The 90 manipulator models are generated by the combination of ten different lengths $\lambda\in\{650~\mbox{mm}, 700~\mbox{mm},\ldots, 1\text{,}100~\mbox{mm}\}$ of the manipulator segment two, and nine different lengths $\mu\in\{600~\mbox{mm}, 650~\mbox{mm}, \ldots, 1\text{,}000~\mbox{mm}\}$ of manipulator segment four. The second and fourth manipulator segments are responsible for the manipulators height such that we have a variety of small to large manipulators. In each measurement configuration the total number of collision avoidance constraints is 186. In the following the 90 simulated parameter identifications are denoted as 90 experiments. The control and geometry of the manipulators are modeled according to the product of homogeneous transformations (1) introduced in Section 1. The Denavit-Hartenberg and Hayati parameters of the notional manipulators, i.e., their zero offset values, are presented in Table 2. We underline the geometric parameters that are considered in the OED problems. In this example, a laser tracker measurement system is used with one reflector which is attached to the TCP of the manipulator. Its position ${^{k}} \boldsymbol {c}_{j}$ in the TCP coordinate system k, cf. Table 3, is not exactly known and will be estimated as well. The pose of the laser tracker is summarized in Table 3 and will also be identified during the parameter identification. The measurement error $\boldsymbol {\epsilon}_{i}$ of the laser tracker depends on the precision with which the laser tracker can determine the wavelength of light in the measurement environment (cf. [34]). In [34] an expanded uncertainty ($k=2$) of $1~\mu\mbox{m} + 0.24~\mu\mbox{m}$ on reference length measurements up to 60 m is achieved using a laser tracker in their tape tunnel facility. Please note that we are performing an in-silico experiment to draw the attention to the performance of our new method. Therefore, we assume

$$ \boldsymbol {\epsilon}_{i} \sim\mathcal{N}_{3} \left ( \left [ \begin{matrix}0 \\ 0 \\ 0 \end{matrix} \right ] , \left [ \begin{matrix}1&0&0\\ 0&1&0\\ 0&0&1 \end{matrix} \right ] \right ) ,\quad i = 1, \ldots, 60 . $$

Figure 6

Experimental setup of parameter identification of industrial manipulators using laser tracker systems and capsule-based collision avoidance. The figure is created with the free available visualization software MeshUp [33].

The ranges of manipulator motions are defined by the upper and lower bounds of the manipulator specific control variables $\theta_{1}, \ldots, \theta_{6}$ listed in Table 4. As we are using a laser tracker system the only system dependent constraint is the restricted angle of incidence between the sensor surface and the laser beam, cf. Table 4.

Table 4 Range of motion of joints and admitted incident angle between sensor surface and laser beam

The constraints which are needed during optimization to avoid collisions are discussed in the following. The work space of the manipulator is constrained by a wall behind the manipulator imposing constraints

$$ \bigl({^{W}_{c} b_{i}^{A}} \bigr)_{x} - r_{i} \ge-450.0, \qquad \bigl({^{W}_{c} b_{i}^{B}} \bigr)_{x} - r_{i} \ge-450.0 $$

$\forall i \in\mathcal{R}$ in x direction. The manipulator is not supposed to hit the ground. Therefore, we impose constraints

$$ \bigl({^{W}_{c} b_{i}^{A}} \bigr)_{z} - r_{i} \ge0.0, \qquad \bigl({^{W}_{c} b_{i}^{B}} \bigr)_{z} - r_{i} \ge0.0 $$

$\forall i \in\mathcal{R}$ that guarantee positive z positions for all manipulator bodies. Auto-collision is avoided by constraints between several manipulator bodies. The distance between the manipulator bodies and the object to which the manipulator is mounted should also be positive. The car, whose body should be measured, is modeled by two capsules, which are not allowed to collide with the manipulator parts.

The laser beam is also modeled by a thin capsule whose distance is compared to the manipulator and the car bodies to guarantee a feasible measurement configuration such that the laser can actually track the reflector. The values of the collision avoiding constraints are shown in Table 5.

Table 5 Values of the constraints to avoid collisions between manipulator bodies, laser beam and object bodies

The radii of the capsules are chosen in a way, that the capsules nearly surround their appropriate manipulator segments. The surrounding capsules are used as a conservative and robust approach to guarantee that the optimized measurement configurations are applicable in reality. Each simulated parameter identification consists of 60 measurement configurations from a heuristic approach and from the optimization method introduced in Section 2.3. The 60 heuristic measurement configurations are randomly and uniformly chosen from the motion range of the industrial manipulator. In the experiment two sets of measurement configurations are computed differing by the choice of variance-covariance matrix projections in the cost function of optimization.

3.3 Computation of optimum experimental designs

3.3.1 KUKA robots (KR15, KR300, KR500)

Firstly, three types of KUKA robots (KR15, KR300, KR500) [35] are considered in simulation with the experimental setting introduced in Section 3.2. Due to the different mechanical shapes also the number of capsules describing the manipulator segments differ and accordingly the number of collision avoidance constraints are 133 for the KR15, 284 for the KR300 and 116 for the KR500. Heuristic and optimized measurement configurations for the parameter identification of three KUKA robots are considered. The measurement configurations are optimized with the two different variance-covariance matrix projections, cf. (25) and with consideration of collision avoidance. The resulting average statistical uncertainties of the TCP precision by either using heuristic or optimized measurement configurations for parameter identification are shown in Table 6.

Table 6 Key performance indicator (KPI) average variances (mm ² ) of three KUKA robots with heuristic and optimized experimental designs for parameter identification

Figure 7 gives a visual representation of the three dimensional TCP’s uncertainty and their two dimensional projections of the linearized 95% confidence regions for every two Cartesian coordinate axes resulted from parameter identification of the KUKA KR15 by using heuristic or optimized measurement configurations with the cost function $\phi( {C}_{s})$.

Figure 7

Error ellipsoid defining the linearized 95% confidence region of the tool center point (TCP) of the working configuration. Error ellipsoid defining the linearized 95% confidence region (with $\chi_{3}^{2} (0.95) = 7.815$) of the TCP of the working configuration of a KUKA KR15 with its two dimensional projections resulting from heuristic $q_{h}$ and optimized $q_{s}$ measurement configurations used in parameter identification. The center of the ellipsoid is the TCP position with its parameter values presented in Table 3. The half-axes $v_{x}$, $v_{y}$, $v_{z}$ of the confidence ellipsoids are the eigenvalues of the covariance matrices ${C}_{s}( \boldsymbol {q}_{h})$ (blue) and ${C}_{s}( \boldsymbol {q}_{s})$ (red) with $\boldsymbol {q}_{s} :=\min_{ \boldsymbol {q}} \operatorname {tr}( {C}_{s}( \boldsymbol {q}))$.

3.3.2 Variety of 90 industrial manipulators

For a quantitative analysis, simulated parameter identifications are performed with 90 notional industrial manipulators. Figure 8 presents a comparison of the average variances of the TCP between the new approach using a KPI and (a) the heuristic approach and (b) the state-of-the-art optimization approach. During optimization collision avoidance is considered. Higher values on the x-axis mean a better result of our approach using a KPI. We also computed optimized experimental designs without consideration of collision avoidance during optimization. 50% to 70% of the measurement configurations lead to collisions so that the designs are not applicable in practice. The minimum and maximum value of the averaged three dimensional variances of the TCP from the 90 different experiments are shown in Table 7.

Figure 8

Kernel density estimations of the position precision of the tool center point (TCP) between the new approach using key performance indicator and (a) the heuristic approach and (b) the state-of-the-art optimization approach. The Histogram presents the difference of TCP precision of working configuration $\boldsymbol {q}^{\mathrm{work}}$ from the new approach using KPI and the heuristic/state-of-the-art optimized measurement configurations $\boldsymbol {q}_{x}^{\mathrm{meas}}$, where $x \in\{{h}, \boldsymbol {p}\}$ for 90 parameter identifications differing in the industrial manipulators’ models. Higher values on the x-axis mean a better result of the new approach. The optimized measurement configurations $\boldsymbol {q}_{p}$ and $\boldsymbol {q}_{s}$ result by solving the optimization problems: $\boldsymbol {q}_{p}^{\mathrm{meas}} :=\min_{ \boldsymbol {q}^{\mathrm{meas}}} \operatorname {tr}( {C}_{p}( \boldsymbol {q}^{\mathrm{total}}))$ and $\boldsymbol {q}_{s}^{\mathrm{meas}} :=\min_{ \boldsymbol {q}^{\mathrm{meas}}} \operatorname {tr}( {C}_{s}( \boldsymbol {q}^{\mathrm{total}}))$.

Table 7 Minimum and maximum value of the averaged three dimensional variances of the tool center point (mm ² ) from the 90 experiments with heuristic and optimized experimental design for parameter identification

In Table 8 the heuristic and optimized average variances of the TCP of the working configuration for the 90 experiments due to the four different experimental designs and the use of the collision avoidance method in the optimization are shown.

Table 8 Key performance indicator (KPI) average variances (mm ² ) of 90 experiments with heuristic and optimized experimental design for parameter identification

4 Discussion

4.1 KUKA robots (KR15, KR300, KR500)

In the simulated parameter identification of the three different KUKA robots, the optimized experimental designs $\boldsymbol {q}_{p}^{\mathrm{meas}}$ and $\boldsymbol {q}_{s}^{\mathrm{meas}}$ provide higher TCP accuracies in contrast to the heuristic experimental design (improvement of 32% to 44%). The largest improvement of the precision of the TCP is achieved by only using the covariance matrix ${C}_{s}$ of the KPI in the cost function for optimization (KR 15: 44%, KR 300: 40%, KR 500: 43%). A comparison between the state-of-the-art optimization approach and our new approach implies that the new problem formulation (24a)-(24e), which takes the TCP as a KPI into account, yields a higher TCP precision (KR 15: 17%, KR 300: 10%, KR 500: 12%). Moreover, our formulation of optimizing measurement configurations can be used for industrial manipulators with a low payload to heavy duty models.

4.2 Variety of 90 industrial manipulators

In order to show that our parameter identification method is not limited to some specific industrial manipulators, we performed a case study with 90 notional industrial manipulators. By analyzing Table 8 and Figure 8 one recognizes that the two optimum experimental designs reduce the statistical uncertainty of the 90 manipulators’ TCP accuracies by 36% and more on average in comparison to the heuristic set of configurations. The collision avoidance method, which is essential for the computation of practicable measurement configurations, has an increasing effect on the optimization performance by 6% to 9%. Table 8 indicates that the largest improvement (43%) of the TCP precision is achieved by only using the statistical uncertainty of the TCP position as cost function. Moreover, the case study stresses that the new approach using KPIs is superior to the state-of-the-art optimization approach providing a 19% higher TCP precision on average. The case study underlines the benefit of the TCP precision achieved through our problem formulation and shows that the improvement is not limited on specific KUKA models but on a range of small to large industrial manipulators.

5 Conclusion

This article introduces a new problem formulation for the computation of optimal and collision-free measurement configurations for parameter identification of industrial manipulators. The novelty lies in the fact, that the precision of the tool center point (TCP) is directly considered in the optimization problem as a key performance indicator (KPI). The approach is verified by the simulated parameter identification of three different KUKA robots and also by the quantitative results with 90 notional manipulator geometries. In the experiments an improvement of 40% to 44% of the precision of the TCP is achieved in contrast to the heuristic approach and 10% to 19% improvement compared to an existing state-of-the-art method. For the computation of collision-free configurations required in practice a collision avoidance method is introduced, which provides a minimal number of nonlinear constraints in the optimization problem. In the experiments a laser tracker system is used for parameter identification but the approach is also applicable to other measurement systems. Furthermore, the approach is not limited to a specific observability index (OI) as cost function. As our approach yields a higher TCP precision with the same number of measurement configurations it is also possible to reduce the needed number of configurations to achieve a certain TCP precision in a shorter re-identification time interval compared to the heuristic approach.

6 Outlook

Despite the fact that our problem formulation with key performance indicators (KPIs) improves the tool center point (TCP) precision by 40%, the parameter identification procedure is not adaptive at the moment and cannot incorporate an online estimation of the TCP precision. The next step to be investigated would be an online approach of parameter estimation and optimum experimental design (OED) for parameter identification of manipulators, which monitors the current precision and sequentially identifies itself using optimal measurement configurations calculated on-the-fly when necessary. We assume that with this approach a fast re-identification of parameters with a guaranteed TCP precision during the work cycles of the manipulator can be realized leading to more flexibility in re-identification and a higher throughput. Additionally, the online parameter identification has the advantage that the parameter identification can be stopped when a certain level of TCP precision is reached which results in shorter parameter identification intervals.

At the moment, the improvements are achieved by considering a kinematic model of the manipulator only. However, the improvement could be increased by incorporating a higher level of detail of the manipulator into the model. The higher level of detail is achieved by introducing further parameters, e.g. non-geometric errors like joint mobilities or elasticities. This would require the computation of forces acting on the links and joints, which can be achieved by using a dynamic model of the manipulator described by differential equations. As shown in [30], the same approach for differential equation models can be used when the required derivative information is available. However, this approach is yet to be investigated for the dynamics of multi-body systems acting as dynamic model for the parameter identification procedure.

Another possible extension is the definition and consideration of further KPIs in the optimization problem. The article demonstrates that the use of TCP precision as KPI achieves a significant precision improvement after parameter identification over the heuristic parameter identification approach as well as the existing state-of-the-art problem formulations not using KPIs. Future work will investigate the possibilities of not only taking one TCP pose into account but additional poses of important working configurations. Furthermore, the definition of KPIs is not limited to TCPs, such that working path trajectories of welding manipulators could be defined as KPIs as well. This way, the range of possible application of the proposed approach will be increased.

Appendix

1.1 A.1 A comparison between constrained and unconstrained least-squares

In the Appendix, we discuss an alternative approach to tackle key performance indicators (KPIs). Due to the (CQ) condition, the constrained least-squares problem (13a)-(13b) can be reformulated as an unconstrained least squares problem by parametrization of the s variables. By application of the implicit function theorem there exists an unique solution

$$ \boldsymbol {s}= \Phi( \boldsymbol {p}) := \bigl\{ \boldsymbol {s}\in \mathbb {R}^{n}_{s} | 0 = {f}_{2}( \boldsymbol {p}, \boldsymbol {s}) \bigr\} . $$

and hence (13a)-(13b) can be written as

$$ \min_{ \boldsymbol {p}\in \mathbb {R}^{n_{p}}} \frac{1}{2} \bigl\Vert {f}_{1} \bigl( \boldsymbol {p}, \Phi( \boldsymbol {p}); \boldsymbol {\eta}, \boldsymbol {q}\bigr) \bigr\Vert _{2}^{2} . $$

The Jacobian of $\tilde{{f}}_{1}( \boldsymbol {p}; \boldsymbol {\eta}, \boldsymbol {q}) := {f}_{1}( \boldsymbol {p}, \Phi( \boldsymbol {p}); \boldsymbol {\eta}, \boldsymbol {q})$ is

$$\begin{aligned} \tilde{ {F}}_{1} &= \frac{\mathrm{d} \tilde{{f}}_{1}}{\mathrm{d} \boldsymbol {p}} = \frac {\partial {f}_{1}}{\partial \boldsymbol {p}} + \frac{\partial {f}_{1}}{\mathrm{d} \boldsymbol {s}}( \boldsymbol {p}, \boldsymbol {s}; \boldsymbol {\eta}, \boldsymbol {q}) \bigg|_{s = \Phi( \boldsymbol {p})} \frac{\boldsymbol {d} \Phi}{ \boldsymbol {d} \boldsymbol {p}} \\ &= \frac{\partial {f}_{1}}{\partial \boldsymbol {p}} + \frac{\partial {f}_{1}}{\partial \boldsymbol {s}}( \boldsymbol {p}, \boldsymbol {s}; \boldsymbol {\eta}, \boldsymbol {q}) \bigg|_{ \boldsymbol {s}= \Phi( \boldsymbol {p})} \biggl( \frac{\partial {f}_{2}}{\partial \boldsymbol {s}} \biggr)^{-1} \frac{\partial {f}_{2}}{\partial \boldsymbol {p}} . \end{aligned}$$

For notational simplicity we define

$$\begin{aligned}& \bigl[\textstyle\begin{array}{@{}c@{\quad}c@{}} A & B \end{array}\displaystyle \bigr] := {F}_{1} = \begin{bmatrix} \frac{\partial {f}_{1}}{\partial \boldsymbol {p}} & \frac{\partial {f}_{1}}{\partial \boldsymbol {s}} \end{bmatrix}, \\& \bigl[\textstyle\begin{array}{@{}c@{\quad}c@{}} E & F \end{array}\displaystyle \bigr] := {F}_{2} = \begin{bmatrix} \frac{\partial {f}_{2}}{\partial \boldsymbol {p}} & \frac{d {f}_{2}}{\partial \boldsymbol {s}} \end{bmatrix} \end{aligned}$$

and thus

$$ \tilde{ {F}}_{1} = A + B F^{-1} E . $$

The covariance matrix of the unconstrained parameter estimation problem is therefore

$$\begin{aligned} \tilde{ {C}}_{p} &= \bigl(\tilde{ {F}}_{1}^{T} \tilde{ {F}}_{1} \bigr)^{-1} \\ &= \bigl( \bigl(A + B F^{-1} E \bigr)^{T} \bigl(A + B F^{-1} E \bigr) \bigr)^{-1} . \end{aligned}$$

To obtain the covariance matrix of the KPI s one can apply linear error propagation and finds

$$ \tilde{ {C}}_{s} = \bigl(F^{-1} E \bigr) \tilde{ {C}}_{p} \bigl(F^{-1} E \bigr)^{T} . $$

The following two propositions show that $\tilde{ {C}}_{p} = {C}_{p}$ and $\tilde{ {C}}_{s} = {C}_{s}$.

Proposition

Let (CQ) and (PD) be satisfied and let the columns of $Z \in\mathbb {R}^{(n_{v}-n_{2}) \times n_{v}}$ span the null-space of ${F}_{2}$, i.e.,

$$ {F}_{2} Z = 0. $$

Then it holds that

$$ \big[\textstyle\begin{array}{@{}c@{\quad}c@{}} \mathbb{I} & 0 \end{array}\displaystyle \bigr] \begin{bmatrix} {F}_{1}^{T} {F}_{1} & {F}_{2}^{T} \\ {F}_{2} & 0 \end{bmatrix} ^{-1} \begin{bmatrix} \mathbb{I} \\ 0 \end{bmatrix} = Z^{T} \bigl(Z {F}_{1}^{T} {F}_{1} Z^{T} \bigr)^{-1} Z . $$

Proof

$$\begin{aligned}& \begin{bmatrix} {F}_{1}^{T} {F}_{1} & {F}_{2}^{T} \\ {F}_{2} & 0 \end{bmatrix} ^{-1} \begin{bmatrix} \mathbb{I} \\ 0 \end{bmatrix} = \begin{bmatrix} C \\ D \end{bmatrix} \end{aligned}$$

(26a)

$$\begin{aligned}& \quad \Leftrightarrow\quad \mathbb{I} = {F}_{1}^{T} {F}_{1}C + {F}_{2}^{T} D, \end{aligned}$$

(26b)

$$\begin{aligned}& \hphantom{\quad \Leftrightarrow\quad\ }0 = {F}_{2} C. \end{aligned}$$

(26c)

From (26c) and the condition (CQ) one concludes that there exists a K of matching dimensions satisfying ${C}= ZK$. Inserting this identity into (26b) and multiplication by $Z^{T}$ from the left one obtains

$$ Z^{T} {F}_{1}^{T} {F}_{1} Z K = Z^{T} . $$

From condition (PD) it follows that $Z^{T} {F}_{1}^{T} {F}_{1} Z$ is non-singular and thus

$$\begin{aligned}& K = \bigl(Z^{T} {F}_{1}^{T} {F}_{1} Z \bigr)^{-1} Z^{T} \\& \quad \Leftrightarrow\quad {C}= Z \bigl(Z^{T} {F}_{1}^{T} {F}_{1} Z \bigr)^{-1} Z^{T} \end{aligned}$$

by multiplication of Z from the left. This is the desired result. □

Proposition

Let ${C}_{p}$, ${C}_{s}$, $\tilde{ {C}}_{p}$ and $\tilde{ {C}}_{s}$ be defined as above. It holds that

$$ \tilde{ {C}}_{p} = {C}_{p} \quad \textit{and} \quad \tilde{ {C}}_{s} = {C}_{s} . $$

Proof

We apply the result of the latter Proposition with the specific matrix

$$ Z^{T} = \begin{bmatrix} \mathbb{I} \\ - F^{-1} E \end{bmatrix} $$

yielding

$$\begin{aligned} {C}_{p} & = \bigl(Z {F}_{1}^{T} {F}_{1} Z^{T} \bigr)^{-1} \\ & = \left( \bigl[\textstyle\begin{array}{@{}c@{\quad}c@{}} \mathbb{I} & (-F^{-1} E)^{T} \end{array}\displaystyle \bigr] ^{T} \begin{bmatrix} A^{T} A & A^{T} B \\ B^{T} A & B^{T} B \end{bmatrix} \begin{bmatrix} \mathbb{I}\\ (-F^{-1} E)^{T} \end{bmatrix} \right)^{-1} \\ & = \left( \bigl[\textstyle\begin{array}{@{}c@{\quad}c@{}} \mathbb{I} & (-F^{-1} E)^{T} \end{array}\displaystyle \bigr] ^{T} \begin{bmatrix} A^{T} A - A^{T} B F^{-1} E \\ B^{T} A - B^{T} B F^{-1} E \end{bmatrix} \right)^{-1} \\ & = \bigl(A^{T} A - \bigl(F^{-1} E \bigr)^{T} B^{T} A + \bigl(F^{-1} E \bigr)^{T} B^{T} B F^{-1} E \bigr)^{-1} \\ & = \bigl( \bigl(A + B F^{-1} E \bigr)^{T} \bigl(A + B F^{-1} E \bigr) \bigr)^{-1} \\ & = \tilde{ {C}}_{p} \end{aligned}$$

and

$$\begin{aligned} {C}_{v} & = Z^{T} \bigl(Z {F}_{1}^{T} {F}_{1} Z^{T} \bigr)^{-1} Z \\ & = \begin{bmatrix} {C}_{p} & - {C}_{p} (F^{-1} E)^{T} \\ -F^{-1} E {C}_{p} & (F^{-1} E) {C}_{p}(F^{-1} E)^{T} \end{bmatrix} . \end{aligned}$$

Projecting onto the s variables results in

$$\begin{aligned} {C}_{s} & = \bigl[\textstyle\begin{array}{@{}c@{\quad}c@{}} 0 & \mathbb{I} \end{array}\displaystyle \bigr] {C}_{v} \begin{bmatrix} 0 \\ \mathbb{I} \end{bmatrix} \\ & = \bigl(F^{-1} E \bigr) {C}_{p} \bigl(F^{-1} E \bigr)^{T} \\ & = \tilde{ {C}}_{s} . \end{aligned}$$

 □