<figure><img src="image_1.png" /></figure>

HAL
open science

# Prediction of pose errors implied by external forces applied on robots: towards a metric for the control of collaborative robots

Vincent Fortineau, Vincent Padois, David Daney

► To cite this version:

Vincent Fortineau, Vincent Padois, David Daney. Prediction of pose errors implied by external forces applied on robots: towards a metric for the control of collaborative robots. IEEE International Conference on Robotics and Automation - ICRA 2024, May 2024, Yokohama, Japan. hal-04410256v2

HAL Id: hal-04410256

https://inria.hal.science/hal-04410256v2

Submitted on 9 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Prediction of pose errors implied by external forces applied on robots: towards a metric for the control of collaborative robots

Vincent Fortineau<fnref n="1" />, Vincent Padois<fnref n="1" />, David Daney<fnref n="1" />

***Abstract***— The presented work tackles the question of quantifying the pose deviations of robots subject to external disturbance forces. While this question may not be central for large robots perfectly rejecting disturbances through high controller gains, it is an important factor when considering collaborative settings where smaller robots may be deviated from their task because of unmodeled physical interactions. This is all the more true with human-robot collaboration where human capacities may fluctuate over time and have to be compensated by a proper adaptation of the robot control. To move forward in this direction, this work first derives a deviation prediction methodology and exemplifies it using three largely employed control approaches. The proposed prediction method is then validated using simulated and real robot experiments both in single and multiple robots cases. The obtained results constitute a stepping stone towards a quantitative metric for robots adapting their behaviour to human motor fluctuation.

## I. INTRODUCTION

Robotic manipulators are usually designed and optimised to operate in fully-predictable or static environments, typically standard manufacturing factories. When considering the potential use of robots in the more dynamic context of agile manufacturing, where humans play a central role, off-line computed and infinitely replayed trajectories are no longer a viable solution. One needs the ability to handle the possible physical interactions of the robot with its environment. This is especially true when considering control modes where such interactions are central, typically co-manipulation or more generally physical human-robot collaboration.

Multiple qualitative and quantitative metrics have been proposed to characterise robots' behaviour with physical interaction and collaboration in mind [1]. This characterisation and the understanding of the interaction behaviour of the couple wrench/trajectory in physical collaborations is important to help designing appropriate robots and their controllers.

Application-wise, one important characteristic is the deviation imposed on the end-effector of a robot by an externally applied wrench. For rigid robots, this deviation is mostly dominated by the controller architecture and gains. Relating these gains to the set of deviations potentially induced by bounded external forces can help to both characterise and modulate the dynamics of robots having to comply with some precision requirements. This raises several questions: given a bound on the external wrench, can we compute a bound on this deviation? Can we adapt the control gains to compensate for this deviation? In the case of the parallel manipulation of an object by two robots, or a robot and a human, can the drop in performance of one member of the dyad be compensated for by the other?

To tackle these questions, it is important to consider control paradigms where the question of the behaviour of the robot under external wrenches is central. This is the case of the impedance control strategy initially proposed in [2] as a way to unify both aspects through the control of the interaction behaviour. The impedance control method is still an active research topic for physical interaction between robots and their environments. There are several examples of impedance modulation in the literature to comply with task constraints [3], adapt to the requirements of a task with various contact phases [4], or in the case of bilateral teleoperation [5].

Yet, no metric actually considers explicitly the implied deviations. In terms of existing metrics, dexterity and transparency indices such as manipulability ellipsoids [6] provide, on the one hand, workspace capacities, but also performances indexes to optimise robot trajectories. On the other hand, stiffness, inertia and damping ellipsoids provides capability characteristics that can also help to modulate the behaviour of the robot in space [7]. Nevertheless, ellipsoid-based approaches remain qualitative approximations. Instead, the considered physical quantities can generally be bounded and the corresponding sets can be exactly represented by polytopes [8]. As an example, one advantage of polytopic representation lies in the ability to extend them to the collaborative context through well-defined intersection and summation operators, when it is a more complex exercise with ellipsoids [9].

In the presence of undesired external forces, [10] proposed the concept of residual force polytopes which qualifies the amount of forces a robot can sustain to avoid interference on its trajectory. [11] introduced feasible wrench polytopes, for legged robots, taking into account not only the capacities of the robots, but also features of the environments such as contact normals. Both these metrics allows to improve the robustness of a robot to external force disturbances. Yet, they neither explicitly define the formal relation between the applied force and the induced pose deviation nor do they experimentally verify the faithfulness of these proposed indicators on real-world robots.

In this work, this relation is explicited in the case of three different control strategies (Section II). It yields a polytopic representation which can be efficiently computed to characterise accurately the deviations, given bounds on the considered disturbance wrench (Section III). The obtained models of deviations are then validated both through simulations (Section IV) and experiments (Section V) on a 7 degrees of freedom (DoF) serial manipulator. The proposed models and experiments also consider the case of two robots in a parallel co-manipulation setting. Results are presented and analysed

<footnote n="1">¹Inria, Talence, France surname.lastname@inria.fr</footnote>
in Section VI. Section VII concludes the paper with some perspectives towards collaborative robotics.

## II. STATIC DEVIATION IN THE PRESENCE OF EXTERNAL WRENCHES

In this section, a model relating an external wrench, acting on a robot at the end-effector level, to the pose error it induces is derived. The transient behaviour induced by the apparition of this wrench is ignored and the final value theorem is used to compute the static pose error for the robot end-effector.

Let's first recall the inverse dynamic model of an n actuated joints serial fixed-base manipulator

$$ \ddot{q} = M^{-1}(q) (\tau_d - J(q)^t f - h(q, \dot{q}) - g(q)), \quad (1) $$

with $q \in \mathbb{R}^n$ the joint position, $\tau_d \in \mathbb{R}^n$ the controlled actuation torques, $f \in \mathbb{R}^6$ an external wrench applied, without loss of generality, on the end-effector of the robot, $M \in \mathbb{R}^{n \times n}$ the robot inertia matrix, $J \in \mathbb{R}^{6 \times n}$ the end-effector Jacobian matrix, $h \in \mathbb{R}^n$ and $g \in \mathbb{R}^n$ the joint torques respectively induced by velocity related inertial terms (namely Coriolis and centrifugal ones) and gravity. For the sake of brevity, the dependence on $q$ or $\dot{q}$ of all matrices is omitted in the rest of the paper.

Given this model and considering torque-based control approaches, three different types of control strategies are compared to explicit their influence on the pose deviations induced by $f$ which is considered as a disturbance. The first is a computed-torque control approach applied at the joint level [12, sec. 11.4.1.3]. The second one is its Cartesian space counterpart. The last one is a typical quasi-static Jacobian transpose controller [12, sec. 11.5].

### A. Computed torque control at the joint level

In the case of a desired joint behaviour, the commanded acceleration $\ddot{q}_d$ can be computed as the sum of a feedforward acceleration term $\ddot{q}^*$ with a feedback control one

$$ \ddot{q}_d = \ddot{q}^* + K_{D_q}(\dot{q}^* - \dot{q}) + K_{P_q}(q^* - q) \quad (2) $$

with, $K_{D_q} \in \mathbb{R}^{n \times n}$ and $K_{P_q} \in \mathbb{R}^{n \times n}$ diagonal gain matrices respectively the derivative and proportional control gain matrices. This control can be seen as a joint impedance control without inertia modification, with the gain matrices tuning the joint damping and stiffness behaviour. Feedback linearisation can be achieved by implementing this command using the estimated model of the robot:

$$ \tau_d = M\ddot{q}_d + h + g \quad (3) $$

Assuming minimal modelling error and injecting successively (2) in (3) and (1) yields

$$ \ddot{\epsilon}_q = M^{-1}J^t f - K_{D_q} \dot{\epsilon}_q - K_{P_q} \epsilon_q, \quad (4) $$

where $\epsilon_q = q^* - q$ is the joint space error.

Writing $\theta = (\epsilon_q^t, \dot{\epsilon}_q^t)^t$ and considering $f$ as a wrench input, (4) can be written in state-space representation

$$ \begin{aligned} \dot{\theta}(t) &= A\theta(t) + Bu(t) \\ y(t) &= C\theta(t) \end{aligned} \quad (5) $$

where the output $y(t)$ is chosen as $\epsilon_q$ and $u(t) = f$. For a step input $u(t)$ of amplitude $f_0$ and considering the static case ($\dot{\theta} = 0$), the application of the final value theorem to MIMO systems yields

$$ \lim_{t \to \infty} y(t) = -CA^{-1}Bf_0 \quad (6) $$

$K_{D_q}$ and $K_{P_q}$ being invertible by definition, the inverse of the state matrix $A$ always exists [13]. Therefore, in this case $\lim_{t \to \infty} \epsilon_q(t) = \hat{\epsilon}_q$ and

$$ \hat{\epsilon}_q = K_{P_q}^{-1} M^{-1} J^t f_0 \quad (7) $$

Note $x \in SE(3)$ the end-effector pose and $x^* = fk(q^*)$ the expected pose with $fk$ the forward kinematic model of the robot. Assuming a sufficiently small variation of the configuration $\epsilon_{qs}$, the static end-effector pose deviation $\hat{\epsilon} \in se(3)$ can be approximated by

$$ \hat{\epsilon} = J K_{P_q}^{-1} M^{-1} J^t f_0 \quad (8) $$

### B. Computed torque control at the Cartesian level

In the case of a desired end-effector behaviour, the desired acceleration $\ddot{x}_d$ can be computed from an end-effector feed-forward acceleration $\ddot{x}^*$, and a feedback control around the desired end-effector pose$^1 x^*$

$$ \ddot{x}_d = \ddot{x}^* + K_{D_x}(\dot{x}^* - \dot{x}) + K_{P_x}(x^* \ominus x) \quad (9) $$

This is analogous to a Cartesian impedance control (9) without inertia modification, with $K_{D_x} \in \mathbb{R}^{6 \times 6}$ and $K_{P_x} \in \mathbb{R}^{6 \times 6}$ diagonal gain matrices respectively tuning the end-effector damping and stiffness behaviour. The operator $\ominus$ represents the difference between two SE(3) elements as detailed in [12, sec. 3.3.3].

Regardless of the method used to compute $\tau_d$ (standard task space inverse dynamics [14] or its Quadratic Programming (QP) equivalent [15]) but assuming that it does not account for the unmodeled wrench disturbance $f$, the control torque that would account for this disturbance can be written $\tau_d + J^t f$. Thus, given the derived Jacobian mapping $\ddot{x} = J\ddot{q} + \dot{J}\dot{q}$, one can compute the task space acceleration as a function of the controlled joint torque using (1). In the case of $\tau_d$, this yields

$$ \ddot{x} = JM^{-1}(\tau_d - J^t f - h - g) + \dot{J}\dot{q}, \quad (10) $$

while when compensating for the external wrench, the acceleration would actually reach the desired one

$$ \ddot{x}_d = JM^{-1}(\tau_d - h - g) + \dot{J}\dot{q} \quad (11) $$

Subtracting (11) from (10), leads to $\ddot{x}_d = \ddot{x} + JM^{-1}J^t f$ which, injected in (9) yields

$$ \ddot{\epsilon} = JM^{-1}J^t f - K_{D_x} \dot{\epsilon} - K_{P_x} \epsilon, \quad (12) $$

with $\epsilon = x^* \ominus x$ and $\dot{\epsilon} = \dot{x}^* - \dot{x}$.

Then the static deviation $\hat{\epsilon}$ is obtained using the final value theorem (6), with a state space representation of (12), using the state $\theta = (\epsilon^t, \dot{\epsilon}^t)^t$ and a step input $f_0$

$$ \hat{\epsilon} = K_{P_x}^{-1} JM^{-1} J^t f_0 \quad (13) $$

<footnote>$^1 \ddot{x}$ is an abuse of notation for the end-effector twist as there is no representation of a pose $x \in SE(3)$ which times derivative is a twist. Yet, for the sake of conciseness, this notation is used throughout this paper.</footnote>
### C. Quasi-static control

Quasi-static control does not rely on the dynamic model of the robot, and is usually applied for near-zero acceleration and velocity contexts ($\ddot{x}^* = \dot{x}^* = 0$).

An equivalent actuation end-effector wrench $f_d$ is computed through a Cartesian proportional derivative feedback controller

$$f_d = K_{D_x}(-\dot{x}) + K_{P_x}(x^* \ominus x) \quad (14)$$

$f_d$ can then be mapped at the torque control level, including the compensation of gravity effects

$$\tau_d = \mathbf{J}^t \mathbf{f}_d + g \quad (15)$$

Injecting (15) in (1) leads to the induced joint space acceleration which in turn can be used to compute the task space acceleration through the derived Jacobian mapping $\ddot{x} = \mathbf{J}\ddot{q} + \dot{\mathbf{J}}\dot{q}$. Noting that $\ddot{x}^* = 0$ in the quasi-static case, $\ddot{\epsilon} = -\ddot{x}$ which leads to

$$\ddot{\epsilon} = \mathbf{J}\mathbf{M}^{-1}(\mathbf{J}^t(\mathbf{f} - \mathbf{K}_{D_x}\dot{\epsilon} - \mathbf{K}_{P_x}\epsilon) - \mathbf{h}) - \dot{\mathbf{J}}\dot{q} \quad (16)$$

Then the static deviation $\hat{\epsilon}$ can be obtained using the final value theorem (6) for a step input $f_0$, on a state space representation of (16) with the state $\theta = (\epsilon^t, \dot{\epsilon}^t)^t$. The quasi-static assumption leads to neglecting $\dot{\mathbf{J}}\dot{q}$ and $\mathbf{h}$ and $\hat{\epsilon}$ is given by

$$\hat{\epsilon} = (\mathbf{J}\mathbf{M}^{-1}\mathbf{J}^t\mathbf{K}_{P_x})^{-1}\mathbf{J}\mathbf{M}^{-1}\mathbf{J}^t\mathbf{f}_0 \quad (17)$$

Except in singular configurations, the apparent mass matrix $(\mathbf{J}\mathbf{M}^{-1}\mathbf{J}^t)^{-1}$ is invertible$^2$. $\mathbf{K}_{P_x}$ is also invertible by essence. Therefore, expression (17) simplifies to

$$\hat{\epsilon} = \mathbf{K}_{P_x}^{-1}\mathbf{f}_0 \quad (18)$$

### III. BOUNDING THE DEVIATIONS

In this section, bounds on the deviation are proposed given bounds on the disturbance wrench both in the case of a single arm and of multiple arms co-manipulating a load in a parallel way. In the latter case, a simple stiffness modulation approach is also proposed.

#### A. Polytopes of deviation

In a given configuration, the static deviations for all three control approaches (8), (13), (18) can be expressed linearly with respect to the disturbance wrench in the form $\hat{\epsilon} = \mathbf{H}\mathbf{f}_0$. $\mathbf{H} \in \mathbb{R}^{6 \times 6}$ is, apart from singularities, a positive definite matrix, homogeneous to a mechanical compliance. Given a bound on the disturbance wrench $\mathbf{f} \in [\mathbf{f}; \bar{\mathbf{f}}]$, one can compute the polytope (more specifically zonotope) bounding the deviation

$$\mathcal{P}_{\hat{\epsilon}} = \{\hat{\epsilon} = \mathbf{H}\mathbf{f} \mid \mathbf{f} \in [\mathbf{f}; \bar{\mathbf{f}}]\} \quad (19)$$

Conversely, given a maximum acceptable bound on deviations $\hat{\epsilon} \in [\hat{\epsilon}; \bar{\hat{\epsilon}}]$ one can compute the polytope bounding the maximum acceptable disturbance wrench

$$\mathcal{P}_f = \{\mathbf{f} = \mathbf{K}\hat{\epsilon} \mid \hat{\epsilon} \in [\hat{\epsilon}; \bar{\hat{\epsilon}}]\}, \quad (20)$$

with $\mathbf{K} = \mathbf{H}^{-1}$ the equivalent end-effector stiffness.

#### B. The case of multiple arms

One may want to connect multiple robotic arms by their end-effector, for example, to collaboratively carry a load. While this reduces the individual workspace of each arm, the end-point global stiffness is increased and corresponds to the sum of the stiffness of each individual arm. This yields the following deviation

$$\hat{\epsilon} = \left(\sum_{i}^{n_r} \mathbf{K}_i\right)^{-1} \mathbf{f} \quad (21)$$

with $\mathbf{K}_i$ is the apparent stiffness of each of the $n_r$ arms taken individually.

#### C. Stiffness modulation

In the case of two robots using the control law described in Section II-B, using (21), the link between the deviation and the disturbance wrench can be written

$$\mathbf{f} = (\Lambda_1 \mathbf{K}_{P_{x,1}} + \Lambda_2 \mathbf{K}_{P_{x,2}}) \hat{\epsilon} \quad (22)$$

with $\Lambda_i = (\mathbf{J}_i \mathbf{M}_i^{-1} \mathbf{J}_i^t)^{-1}$ the apparent mass matrix of the $i$th robot. As a consequence, to compensate for a variation $\Delta \mathbf{K}_{P_{x,1}}$ of the end-effector control gains of robot 1 and maintain a constant static deviation $\hat{\epsilon}$, robot 2 should modify its control gains considering the apparent mass $\Lambda_1$ of the first robot. Writing

$$\Lambda_1(\mathbf{K}_{P_{x,1}} - \Delta \mathbf{K}_{P_{x,1}}) + \Lambda_2(\mathbf{K}_{P_{x,2}} - \Delta \mathbf{K}_{P_{x,2}}) = \Lambda_1 \mathbf{K}_{P_{x,1}} + \Lambda_2 \mathbf{K}_{P_{x,2}} \quad (23)$$

yields

$$\Delta \mathbf{K}_{P_{x,2}} = -\Lambda_2^{-1} \Lambda_1 \Delta \mathbf{K}_{P_{x,1}} \quad (24)$$

While this result may not seem very useful for two robots, it becomes interesting when considering human-robot collaboration and the adaptation of the robot to human fatigue [4].

### IV. SIMULATED ANALYSIS

In this section, the static deviation $\hat{\epsilon}$ is numerically evaluated in simulation and compared to the theoretical results of Section II. First, the polytopic representation are illustrated. Then the results of a systematic comparison of the simulated and predicted error is presented for the case of the task space inverse dynamics control approach (Section II-B).

The prediction of the error based on the results of Section II is computed using the Pinocchio library [16]. The dynamic simulations of the Franka Emika Panda robot are performed using both the Pinocchio library and Gazebo.

In the remainder of the paper, the presented results are restricted, without loss of generality but for the sake of conciseness, to disturbance forces and deviations in terms of position only.

#### A. Spatial representation of static deviations

Fig. 1 shows, for each control strategy described in Section II, the displacement polytopes for 2 random configurations and arbitrarily chosen gains. The considered external force applied at the end-effector level is bounded in a unit cube.

As expected, the shapes, orientations and magnitudes obtained are configuration-dependant only for the first two control strategies.

<footnote>$^2$in the singular case, the apparent mass tends to infinity in which case no pose deviation is possible. This is well rendered by equation (17).</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1: Deviation polytopes in mm (blue), for a unit cubic input force, with each control strategy on two random configurations [1.91, 1.08, -0.13, -1.18, 0.03, 2.19, 1.87] (first row) and [-2.21, -1.64, -0.35, -1.81, -1.64, 1.67, 0.95] (second row). Joint control gain: $K_{P_q} = \text{diag}(1, 2, 3, 4, 5, 6, 7)$, end-effector control gain: $K_{P_x} = \text{diag}(1, 2, 3)$, quasi-static control gain $K_{P_x} = I_3$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2: Error profile (solid blue) of the deviation norm $||\hat{\epsilon}|| - ||\epsilon_s||$, for all tested configurations with both -20 N and 20 N simulated force inputs on the $\vec{x}$ axis. The interquartile range is coloured in blue, while the light blue dotted area represents 90% of the data between the 5% and 95% percentiles. Zoom ratios are indicated next to the zoomed windows. The force input starting front is signified with the vertical dashed purple line.</figcaption></figure>

Given the expressions of the static deviations for each control mode, one could compute configuration-dependant gains that yield similar deviations for all three modes. This is not illustrated here for the sake of brevity.

## B. Analysis of deviations across configurations and force amplitudes

In this section, a quantitative analysis of deviations across configurations and force amplitudes for the task-spaced inverse dynamics control mode is presented. A thousand random configurations are generated within the joint boundaries of the robot among which only 767 are kept, leaving out configurations with auto-collision, reaching joint limits or where the robot could not stably hold the pose. Two magnitudes of perturbation (20 N and -20 N) are tested in all directions ($\vec{x}, \vec{y}$ and $\vec{z}$), for a total of 6 sets of perturbations.

The median and interquartile range of the deviation error profile ($||\hat{\epsilon}(t)||_2 - ||\epsilon_s(t)||_2$), with $\epsilon_s$ the simulated deviation, is computed for all tested configurations and force disturbances. The median errors of the Euclidean displacement norms are lesser than 1 µm for the perturbations on all axes (0.16 µm, 0.20 µm and 0.10 µm, respectively for the force perturbations on $\vec{x}, \vec{y}, \vec{z}$), with interquartile range below 10 µm (4.61 µm, 5.67 µm and 3.59 µm). Fig. 2 illustrates these results for force disturbances applied along the x axis of the reference frame of

the robot.

This shows the adequacy of the predicted error with respect to the simulated behaviour for the static phase. The error profile is obviously greater during the transient phase (reaching more than 1 mm), since the estimation is not designed to predict this behaviour. Yet, this approximation can be considered valid for quasi-static movements.

## V. EXPERIMENTAL SETUP

To confront the predicted deviation to real-world experiments where the robot behaviour differs necessarily from the rigid body dynamic model of the system, two different experiments are conducted on a robot controlled using a computed torque control at the Cartesian level (Section II-B). First, force perturbations are applied to evaluate repeatedly several configurations assessing the reliability of the estimation. Then, two robots are connected rigidly and perturbed with external force to assess the estimation for multiple robots.

### A. Robot and associated controller

A 7 DoF Franka Emika Panda robot is used for these experiments. The external forces, when measured, are acquired using a 6-axis ATI Axia80-m20 sensor. The control strategy is implemented using the FrankaRos library.

In practice, computed torque control at the Cartesian level is implemented using a QP formulation [15] running at 1 kHz. At each control time step, the desired joint torque is computed by solving an optimisation problem in which the cost function aims at minimising the Cartesian acceleration error

$$ \tau_d = \arg \min_{\tau} \quad \begin{aligned}[t] & ||\ddot{\vec{x}}^* - J M^{-1} (\tau - h - g) - \dot{J} \dot{q}||_2^2 + \\ & \omega_{reg} ||\tau - g - k_d \dot{q}||_2^2 \end{aligned} \quad (25) \\ \text{s.t.} \quad B\tau \le b $$

where $\vec{x}^*$ is given by (9) and the linear inequality constraint is used to enforce joint-level position, velocity and torque limits. A regularisation cost ($\omega_{reg} \ll 1$) ensures the uniqueness of the solution for redundant robots without affecting the Cartesian task in practice.

### B. Single Robot experiments

Similarly to the validation of the predicted deviation using simulations, the deviation prediction error is compared to the deviation measured on the real robot for four different configurations related to the corner of the robot isocube in order to cover at best the robot's workspace. Five magnitudes of force disturbances are applied at the end-effector level: -25 N, -10 N, 10 N, 25 N and 40 N along the $\vec{x}, \vec{y}, \vec{z}$ axis of the robot reference frame.

To limit the influence on the results of non-linear effects notably coming from friction, a sequence of movement is defined between each measurement: (a) the robot moves to a neutral configuration $q_0$, which is defined as the mean between each joint limits, (b) the robot moves to the desired predefined configuration, (c) the data recording starts while the robot pauses for 1.5 s, (d) a force step of magnitude m and duration 2 s is generated on either one of the three axes: $\vec{x}, \vec{y}$,
<figure><img src="image_5.png" /><figcaption>Fig. 3: The three configurations of two panda arms connected through a mechanical interface. The first one (a) is a quasi-neutral configuration while the other two (b) and (c) are asymmetric configurations. A hook allows to add a weight jointly carried by both robots.</figcaption></figure>

$\tilde{z}$, (e) the robot pauses for 1.5 s and the recording stops. The sequence is repeated ten times for each unique combination of configuration, perturbation magnitude and axis.

To automate the modification of the external disturbance, an equivalent joint torque is added to the control torque $\tau_d$ rather than actually applying the force at the end-effector level. The equivalence of both methods has been validated experimentally (thus the use of a 6-axis FT sensor) but these results are not shown here for the sake of conciseness.

### C. Two robots experiments

To extend the previous experiments, two robots are connected by their end-effector using a fixed mechanical interface with a hook to add a load. The main idea of this experiment is, here again, to validate by experimentation the validity of the predicted deviation in the case of two robots carrying a shared load.

The robots are manually placed into three different configurations using gravity compensation control. The first one, shown on Fig. 3-a, is chosen to resemble the neutral configuration for both robots. The other two configurations (Fig. 3-b and -c) are defined so that the robots have a non-symmetric shape with an expected non-equal force distribution.

For each configuration, multiple control gains are tested. The gains tested are defined as a scalar multiplied by the identity matrix, so that only one value is attributed per gain matrix. The proportional gain $K_{P_x}$ is varied with the following paired values for each robot: (500, 500), (250, 250), (500, 250), (250, 500), (750, 250), (250, 750).

The measurement procedure is similar to the one described for the single robot case. However, in this case, the external disturbance is applied at the end-effector level using weights. Indeed, in the multiple robot case, the distribution of the added load between the two robots is harder to predict.

## VI. RESULTS

### A. Single robot experiments

The data collected from the experimental process described in Section V-B is compiled in the Table I. Despite the procedure, servoing errors are still observed. For a proper analysis, the servoing error appearing before the perturbation is excluded to disregard the effects of friction. Such effects can still be observed after the perturbation, since the new resting position is different from the initial one, as shown on Fig. 4 and 5.

<figure><img src="image_6.png" /><figcaption>Fig. 4: Average (solid lines) and standard deviations (colored area) of 10 measurements with estimated (dashed) deviation on each axis, for the 5 perturbations magnitudes on the $\tilde{z}$ axis (40 N, 25 N, 10 N, -10 N and -25 N, respectively in blue, orange, green, purple and grey).</figcaption></figure>

Over all data, the deviation prediction errors ($\hat{\epsilon} - \epsilon$) observed on each axis are about $0.14 \pm 2.12$ mm. Considering the whole displacement norm, these errors ($\|\hat{\epsilon}\|_2 - \|\epsilon\|_2$) are on average reaching $2.35 \pm 1.67$ mm, meanwhile the expected displacements $\|\hat{\epsilon}\|_2$ average $12.21 \pm 6.97$ mm.

It can be noted in Fig. 4 that the measured and predicted deviations are not perfectly symmetric, while the end-effector apparent mass $\Lambda(q) = (J M^{-1} J^t)^{-1}$, modulated by the inverse of the end-effector proportional gain $K_{P_x}$ is symmetric for a given configuration $q$. The non-symmetric behaviour is observed since the end-effector apparent mass is evaluated in the actual configuration $q = q_0 + \delta_q$ to improve the accuracy, rather than in the initial configuration $q_0$, with $\delta_q$ small joint variations induced by the interaction force.

Finally, while the obtained results validate the proposed deviation prediction model, the results tend to demonstrate an overestimation of the static deviation, compared to experiments with the robot (except in 5% of all cases where negative errors are all below 1 mm). This could be explained by a greater
<table><caption>TABLE I: Summary of the repeated experiments described in Section V-B. The average deviation norms $||\epsilon||$ are indicated for each configuration and each force perturbation magnitude, without the servoing error. Standard deviations are indicated as confidence intervals. The perturbation axis is specified in each cell. The expected results are indicated right after, in the sub-column $||\hat{\epsilon}||$, and the errors between the estimated and average measured deviations are indicated in the sub-column "err". The results indicated in bold refer to the experiment represented on Fig. 4.</caption><thead><tr><th rowspan="2">config.</th><th colspan="3">40 N</th><th colspan="3">25 N</th><th colspan="3">10 N</th><th colspan="3">-10 N</th><th colspan="3">-25 N</th></tr><tr><th>||ε|| ± std<br>(mm)</th><th>||ε̂||<br>(mm)</th><th>err<br>(mm)</th><th>||ε|| ± std<br>(mm)</th><th>||ε̂||<br>(mm)</th><th>err<br>(mm)</th><th>||ε|| ± std<br>(mm)</th><th>||ε̂||<br>(mm)</th><th>err<br>(mm)</th><th>||ε|| ± std<br>(mm)</th><th>||ε̂||<br>(mm)</th><th>err<br>(mm)</th><th>||ε|| ± std<br>(mm)</th><th>||ε̂||<br>(mm)</th><th>err<br>(mm)</th></tr></thead><tbody><tr><td>1</td><td>fx</td><td>24.94±0.11</td><td>27.6</td><td>2.7</td><td>15.18±0.13</td><td>17.4</td><td>2.2</td><td>6.06±0.10</td><td>7.1</td><td>1.0</td><td>5.25±0.09</td><td>7.1</td><td>1.8</td><td>15.40±0.16</td><td>17.8</td><td>2.4</td></tr><tr><td>1</td><td>fy</td><td>24.09±0.06</td><td>25.5</td><td>1.4</td><td>14.27±0.07</td><td>15.9</td><td>1.6</td><td>4.86±0.09</td><td>6.3</td><td>1.5</td><td>3.50±0.04</td><td>6.3</td><td>2.8</td><td>11.03±0.11</td><td>15.6</td><td>4.6</td></tr><tr><td>1</td><td><b>fx</b></td><td><b>7.22±0.08</b></td><td><b>8.3</b></td><td><b>1.0</b></td><td><b>4.60±0.07</b></td><td><b>5.3</b></td><td><b>0.7</b></td><td><b>1.56±0.03</b></td><td><b>2.1</b></td><td><b>0.6</b></td><td><b>0.65±0.05</b></td><td><b>2.2</b></td><td><b>1.6</b></td><td><b>3.25±0.06</b></td><td><b>5.6</b></td><td><b>2.4</b></td></tr><tr><td>2</td><td>fx</td><td>24.90±0.07</td><td>27.7</td><td>2.8</td><td>14.43±0.07</td><td>17.2</td><td>2.8</td><td>5.47±0.08</td><td>6.8</td><td>1.4</td><td>4.93±0.07</td><td>6.8</td><td>1.8</td><td>13.90±0.13</td><td>16.8</td><td>2.8</td></tr><tr><td>2</td><td>fy</td><td>22.44±0.10</td><td>21.7</td><td>-0.8</td><td>14.25±0.07</td><td>13.7</td><td>-0.6</td><td>5.78±0.08</td><td>5.5</td><td>-0.3</td><td>4.46±0.10</td><td>5.6</td><td>1.1</td><td>12.56±0.11</td><td>14.2</td><td>1.6</td></tr><tr><td>2</td><td>fx</td><td>18.35±0.13</td><td>22.6</td><td>4.3</td><td>10.34±0.11</td><td>14.4</td><td>4.0</td><td>3.72±0.04</td><td>5.8</td><td>2.1</td><td>5.32±0.14</td><td>6.0</td><td>0.7</td><td>14.00±0.15</td><td>15.1</td><td>1.1</td></tr><tr><td>2</td><td>fy</td><td>18.30±0.04</td><td>20.8</td><td>2.5</td><td>11.02±0.02</td><td>13.3</td><td>2.3</td><td>4.38±0.05</td><td>5.4</td><td>1.0</td><td>2.16±0.06</td><td>5.5</td><td>3.4</td><td>8.31±0.08</td><td>14.0</td><td>5.7</td></tr><tr><td>3</td><td>fy</td><td>15.80±0.08</td><td>19.6</td><td>3.8</td><td>9.09±0.05</td><td>12.1</td><td>3.0</td><td>2.87±0.05</td><td>4.8</td><td>1.9</td><td>2.71±0.04</td><td>4.7</td><td>2.0</td><td>8.08±0.05</td><td>11.5</td><td>3.5</td></tr><tr><td>3</td><td>fx</td><td>17.68±0.08</td><td>24.6</td><td>7.0</td><td>9.56±0.07</td><td>15.2</td><td>5.7</td><td>3.35±0.03</td><td>6.0</td><td>2.7</td><td>4.50±0.09</td><td>6.0</td><td>1.5</td><td>13.29±0.04</td><td>14.7</td><td>1.4</td></tr><tr><td>3</td><td>fx</td><td>20.12±0.07</td><td>21.0</td><td>0.9</td><td>12.72±0.05</td><td>13.3</td><td>0.6</td><td>5.10±0.13</td><td>5.4</td><td>0.3</td><td>1.16±0.04</td><td>5.5</td><td>4.3</td><td>7.36±0.13</td><td>13.7</td><td>6.4</td></tr><tr><td>4</td><td>fy</td><td>16.38±0.14</td><td>17.5</td><td>1.2</td><td>10.30±0.19</td><td>11.0</td><td>0.7</td><td>3.79±0.11</td><td>4.5</td><td>0.7</td><td>1.71±0.09</td><td>4.5</td><td>2.8</td><td>6.98±0.19</td><td>11.3</td><td>4.3</td></tr><tr><td>4</td><td>fx</td><td>22.51±0.12</td><td>27.9</td><td>5.4</td><td>12.68±0.03</td><td>17.5</td><td>4.8</td><td>4.51±0.04</td><td>7.0</td><td>2.5</td><td>4.25±0.08</td><td>7.0</td><td>2.7</td><td>14.20±0.06</td><td>17.4</td><td>3.2</td></tr></tbody></table>

dissipation of the energy injected into the robot, which is compatible with the fact that friction effects were neglected.

## B. Two robots experiments

This section describes the results of the experiment detailed in Section V-C, with a weight of 2 kg. To minimise the error induced by calibration issues between both robots, the deviation considered is computed from the average end-effector deviation of each robot $\epsilon_{1,2}$. The estimated deviation of the middle point between both robots $\hat{\epsilon}_{1,2}$ is computed using (22), with the sum of the end-effector forces estimated for each robot.

<figure><img src="image_7.png" /><figcaption>Fig. 5: Predicted (orange) and measured (blue) deviation norms of the middle point of both robots in a quasi-neutral configuration for three paired proportional end-effector stiffness gains (indicated in parenthesis on the graphics), when subject to an external force about 20 N on the $\bar{z}$ axis.</figcaption></figure>

The estimation error is then computed considering the difference between the estimated static deviation $\hat{\epsilon}_{1,2}$ and the actual deviation $\epsilon_{1,2}$. For each experience, the static estimation error is considered using the average value of the plateau phase of both estimated and measured deviations.

Over the three configurations, each repeated 6 times for all tested end-effector gains, the error between the measured and estimated static deviation reaches $0.12 \pm 0.22$ mm.

Fig. 5 shows the results of the deviation norm in the quasi-neutral configuration for two paired end-effector stiffness gains $K_{P_x}$, plus a compensation using the results of (24). As expected, when one of the robots reduces its gains, the deviation is increased. Moreover, the gain loss compensation

yields a deviation close to the original with a difference of $0.34 \pm 0.03$ mm over the measured plateau phase.

## VII. CONCLUSION

This paper shows the evaluation of the prediction of static deviations when a robot is subject to external forces. Despite the proposed linear approximation, the results show a good approximation with errors of the displacement norm below 2.5 mm on average. The deviation predicted is actually an overestimation in 95% of the tested configurations, which could be explained by static friction effects.

These results are extended to the multiple robots case and validated experimentally in the case of two robots. The application of these results to human-robot collaboration is not straightforward because of the rather complex research problem related to the online estimation of human arm stiffness.

Yet, this work provides a step towards a quantitative metric that could serve to analyse human-robot collaborative tasks and to modulate the robot behaviour as a function of the human state. In that respect, the proposed metric could be improved by considering explicitly the robot joint limits in the error model, as suggested in [17]. This is the main future direction of this work.

## REFERENCES

[1] J. Angeles and F. C. Park, “Performance evaluation and design criteria,” in *Springer Handbook of Robotics*, B. Siciliano and O. Khatib, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2008, pp. 229–244. DOI: 10.1007/978-3-540-30301-5_11.

[2] N. Hogan, “Impedance Control: An Approach to Manipulation: Part I—Theory,” *Journal of Dynamic Systems, Measurement, and Control*, vol. 107, no. 1, pp. 1–7, Mar. 1985. DOI: 10.1115/1.3140702.

[3] J. Lachner, F. Allmendinger, S. Stramigioli, and N. Hogan, “Shaping Impedances to Comply With Constrained Task Dynamics,” *IEEE Transactions on Robotics*, vol. 38, no. 5, pp. 2750–2767, Oct. 2022. DOI: 10.1109/TRO.2022.3153949.
[4] L. Peternel, N. Tsagarakis, D. Caldwell, and A. Ajoudani, "Robot adaptation to human physical fatigue in human-robot co-manipulation," en, *Autonomous Robots*, vol. 42, no. 5, pp. 1011-1021, Jun. 2018. DOI: 10.1007/s10514-017-9678-1.

[5] L. M. Doornebosch, D. A. Abbink, and L. Peternel, "Analysis of Coupling Effect in Human-Commanded Stiffness During Bilateral Tele-Impedance," *IEEE Transactions on Robotics*, vol. 37, no. 4, pp. 1282-1297, Aug. 2021. DOI: 10.1109/TRO.2020.3047064.

[6] T. Yoshikawa, "Dynamic manipulability of robot manipulators," in *1985 IEEE International Conference on Robotics and Automation Proceedings*, vol. 2, Mar. 1985, pp. 1033-1038. DOI: 10.1109/ROBOT.1985.1087277.

[7] A. Ajoudani, N. G. Tsagarakis, and A. Bicchi, "On the role of robot configuration in Cartesian stiffness control," in *2015 IEEE International Conference on Robotics and Automation (ICRA)*, May 2015, pp. 1010-1016. DOI: 10.1109/ICRA.2015.7139300.

[8] P. Chiacchio, Y. Bouffard-Vercelli, and F. Pierrot, "Evaluation of force capabilities for redundant manipulators," in *Proceedings of IEEE International Conference on Robotics and Automation*, vol. 4, Apr. 1996, 3520-3525 vol.4. DOI: 10.1109/ROBOT.1996.509249.

[9] P. Chiacchio, S. Chiaverini, L. Sciavicco, and B. Siciliano, "Task Space Dynamic Analysis of Multiarm System Configurations," en, *The International Journal of Robotics Research*, vol. 10, no. 6, pp. 708-715, Dec. 1991. DOI: 10.1177/027836499101000608.

[10] H. Ferrolho, W. Merkt, C. Tiseo, and S. Vijayakumar, "Residual Force Polytope: Admissible Task-Space Forces of Dynamic Trajectories," *Robotics and Autonomous Systems*, vol. 142, p. 103814, Aug. 2021. DOI: 10.1016/j.robot.2021.103814.

[11] R. Orsolino, M. Focchi, C. Mastalli, H. Dai, D. G. Caldwell, and C. Semini, "Application of Wrench-Based Feasibility Analysis to the Online Trajectory Optimization of Legged Robots," en, *IEEE Robotics and Automation Letters*, vol. 3, no. 4, pp. 3363-3370, Oct. 2018. DOI: 10.1109/LRA.2018.2836441.

[12] K. M. Lynch and F. C. Park, *Modern Robotics: Mechanics, Planning, and Control*. Cambridge University Press, 2017.

[13] T.-T. Lu and S.-H. Shiou, "Inverses of 2 × 2 block matrices," *Computers & Mathematics with Applications*, vol. 43, no. 1, pp. 119-129, Jan. 2002. DOI: 10.1016/S0898-1221(01)00278-4.

[14] O. Khatib, "A unified approach for motion and force control of robot manipulators: The operational space formulation," *IEEE Journal on Robotics and Automation*, vol. 3, no. 1, pp. 43-53, 1987.

[15] L. Joseph, V. Padois, and G. Morel, "Towards x-ray medical imaging with robots in the open: Safety without compromising performances," in *2018 IEEE International Conference on Robotics and Automation (ICRA)*, IEEE, 2018, pp. 6604-6610.

[16] J. Carpentier, G. Saurel, G. Buondonno, et al., "The Pinocchio C++ library – A fast and flexible implementation of rigid body dynamics algorithms and their analytical derivatives," en, Jan. 2019. DOI: 10.1109/SII.2019.8700380.

[17] A. Skuric, V. Padois, and D. Daney, "Approximating robot reachable space using convex polytopes," in *15th International Workshop on Human-Friendly Robotics*, Delft, Netherlands, Sep. 2022. DOI: 10.1007/978-3-031-22731-8_4.