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

HAL
open science

# Impact-Aware Task-Space Quadratic-Programming Control

Yuquan Wang, Niels Dehio, Arnaud Tanguy, Abderrahmane Kheddar

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

## To cite this version:

Yuquan Wang, Niels Dehio, Arnaud Tanguy, Abderrahmane Kheddar. Impact-Aware Task-Space Quadratic-Programming Control. 2023. hal-02741682v5

HAL Id: hal-02741682

https://hal.science/hal-02741682v5

Preprint submitted on 6 Sep 2023

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.

Copyright
# Impact-Aware Task-Space Quadratic-Programming Control

Journal Title
XX(X):1–16
©The Author(s) 2023
Reprints and permission:
sagepub.co.uk/journalsPermissions.nav
DOI: 10.1177/ToBeAssigned
www.sagepub.com/

SAGE

Yuquan Wang, Niels Dehio, Arnaud Tanguy, and Abderrahmane Kheddar

## Abstract

Robots usually establish contacts at rigid surfaces with near-zero relative velocities. Otherwise, impact-induced energy propagates in the robot's linkage and may cause irreversible damage to the hardware. Moreover, abrupt changes in task-space contact velocity and peak impact forces also result in abrupt changes in robot joint velocities and torques; which can compromise controllers' stability, especially for those based on smooth models. In reality, several tasks would require establishing contact with moderately high velocity. We propose to enhance task-space multi-objective controllers formulated as a quadratic program to be resilient to frictional impacts in three dimensions. We devise new constraints and reformulate the usual ones to be robust to the abrupt joint state changes mentioned earlier. The impact event becomes a controlled process once the optimal control search space is aware of: (1) the hardware-affordable impact bounds and (2) analytically-computed feasible set (polyhedra) that constrain post-impact critical states. Prior to and nearby the targeted contact spot, we assume, at each control cycle, that the impact will occur at the next iteration. This somewhat one-step preview makes our controller robust to impact time and location. To assess our approach, we experimented its resilience to moderate impacts with the Panda manipulator and achieved swift grabbing tasks with the HRP-4 humanoid robot.

## Keywords

Intentional impact tasks, Impact-aware control, Optimization-based control.

## 1 Introduction

When (rigid) robots collide – intentionally or inadvertently – with a rigid surface with a fairly high relative speed, the induced forces are impulsive and the contact state is uncertain. The shock propagates through the robots’ linkages into the joints and can severely damage some parts of the hardware, e.g., the harmonic gears, weak linkages and/or torque sensors (if any).

A common remedy is to carefully plan contact transitions with near-to-zero relative speed. However, this strategy can not achieve specific tasks such as walking or jumping humanoids, hammering, and swift grabbing. For such tasks, the following robotic issues shall be improved concurrently: (i) design impact-resilient hardware (this is not tackled in this paper); and (ii) devise robust control strategies that switch the robot equations of motions and subsequent controllers following a transition policy called reset maps. The switching often requires a precise impact model and knowledge of additional parameters that depend on the environment and the robot, e.g., the impact localization on the robot (and on the environment surface), the contact normal, and the exact impact time. Acquiring these parameters *in-situ*, instantaneously, and reliably is not always possible in practice.

An impact event is instantaneous and too short for a robot to react effectively. The impact duration (i.e., time interval) depends on the particular contact properties and the robot controller. In our recent studies (Wang et al. 2022b), an impact lasted about 20 ms. This is the reason why, any impact-friendly control strategy shall act a *priori* and a *posteriori*.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** The impact-aware QP regulates the contact velocity in a modified search space to ensure that the post-impact state jumps are hardware-affordable.</figcaption></figure>

To tackle on-purpose impact tasks safely, we challenge the possibility of building on our existing task-space control

Y. Wang is with the Department of Advanced Computing Sciences, Maastricht University, Maastricht, The Netherlands.

A. Tanguy and A. Kheddar are with the CNRS-University of Montpellier, LIRMM, Montpellier, France. Y. Wang and N. Dehio were at this laboratory when this work was conducted.

A. Kheddar is also with the CNRS-AIST Joint Robotics Laboratory, IRL, Tsukuba, Japan.

N. Dehio is with KUKA Deutschland GmbH, Zugspitzstrasse 140, Augsburg, Germany.

**Corresponding author:**

Yuquan Wang Email: yuquan.wang@maastrichtuniversity.nl
framework* instead of devising an entirely new control scheme. There are many reasons for this choice: (i) the framework has been proven to be efficient in handling complex industrial scenarios (Kheddar et al. 2019) and multi-robot control (Bouyarmane et al. 2019); (ii) to our best knowledge, on-purpose impact objectives using the task-space quadratic programing (QP) formalism was not explicitly addressed; (iii) there is a relatively large community using it beyond the authors' circle. Our guiding quest is to try as much as possible to enhance task-space QP control formalisms to deal with impact without considerably changing its structure. In other words: can we envision handling impact tasks simply by adding or reformulating task objectives and constraints without introducing new decision variables? As we see later in this paper, the answer to this latter question is *yes, to some extent*.

Formulating hardware impact limits† as additional QP constraints is straightforward but insufficient. The main problem is dealing with unexpected state jumps that may damage the robot. Unexpected because, in practice, even if an impact is planned to occur, there will be uncertainties on both contact time and location. The feedback velocity and force jump compromise the constraints' feasibility. Consequently, we reformulate the usual constraints to account for such discontinuities. Then the QP is impact-aware and robust to state jumps by modifying the search space according to a one-step-ahead prediction of nearby intended (i.e., expected) impacts. As a result, the controller updates the optimized – and hence feasible – impact velocity reference in every control cycle. If the impact happened, the robot would fulfill both the hardware resilience and task-dictated constraints bounds.

We summarize our main contributions as follows:

C.1 bound the *post-impact states* (Sec. 4) with analytically-constructed convex sets (half-space represented polyhedra), assuming the impact is frictional in three dimensions;

C.2 bound a generic post-impact robot state with a closed-form *impact-aware template constraint*, (Sec. 5.1);

C.3 assess our *impact-aware control design* (Sec. 5.3) through experiments on the Panda manipulator and the humanoid robot HRP-4 (Sec. 6).

Our approach builds on top of our initial concept proposed in simulation for fixed-base robots in *Wang and Kheddar (2019)* and the preliminary extension to the floating-base robot in *Wang et al. (2019)*. This work focuses primarily on task-space control, specifically addressing the incorporation of various sets of post-impact states. The modeling part of our controller relies on the assumptions and impact models derived from the extensive experimental studies in (Wang et al. 2022a,b); where we reported hundreds of impact experiments to benchmark different initial poses, impact joint configurations, and impact velocities for (kinematic-controlled) robots that are typically controlled in joint velocity (or position).

# 2 Related Work

Second order equations of motion, which are commonly-used for robot control and planning, cannot effectively capture the impact dynamics. Thus, the robotics community has developed dedicated models to predict the impact-induced state jumps. These models are derived from non-smooth mechanics theory, see e.g., *Stronge (2000)*. Sec. 2.1 summarizes main models used in robotics starting from the seminal paper by *Zheng and Hemami (1985)* in the 1980s. Given the impact mechanics, Sec. 2.2 briefly reports on mitigating impacts by few ideas in absorbing shock through hardware design and a large part of control strategies in various robotic impact tasks.

## 2.1 Impact modeling

Because robots are commonly made from structures and links that are rigid; and because of the real-time constraints in robot control and planning, impulsive motions are modeled based on rigid-body dynamics in contrast to the computationally demanding approaches such as stress wave or finite element deformation models applied by *Pashah et al. (2008)*.

Recently, *Wang et al. (2022a)* found that for kinematic-controlled robots, it is more accurate to formulate the momentum-conservation equations in the task space and assume the robot to behave like a composite-rigid body during the impact event. Formulations relying on joint-space generalized momentum, e.g., *Zheng and Hemami (1985)*; *Lankarani (2000)*; *Khulief (2013)*, do not account for high-stiffness joints and apply for compliant robots, e.g., pure-torque controlled robots or partially actuated pendulums.

In the impulse space, the momentum conservation defines the *plane of compression*, and Coulomb's friction law defines the *friction cone*, see *Stronge (2000)*; their intersection points conservatively bound all the candidate impulses of an impact.

When the impulse is two-dimensional, post-impact contact modes are limited (and numerable), i.e., bounce, sliding, reverse sliding, and sticking. Checking the intersection points against the candidate contact modes leads to an analytical impact model. The state-of-the-art impact mechanics (Jia et al. 2019; Lankarani 2000; Khulief 2013; Wang and Mason 1992) refer to this strategy as Routh's graphical approach (Routh 1955). When legged robots undergo multiple impacts, post-impact velocity and sliding direction are ambiguous (Remy 2017). Halm and Posa (2019, 2021) recently extended Routh's impact model for multiple planar impacts by reformulating it as differential inclusions. The novel formulation guarantees the solution's existence. Casting the impact model as a linear complementarity problem, Halm and Posa (2021) developed a probabilistically-complete algorithm for approximating the post-impact velocity set.

However, many tasks can not escape from their full three-dimensional space dimension. Moving from 2D to 3D, makes it impossible to determine the post-impact tangential velocity and impulse without numerical integrations, as observed by *Stronge (2000)* and *Jia and Wang (2017)*.

<footnote>*mc_rtc https://github.com/jrl-umi3218/mc_rtc</footnote>

<footnote>†Very few robot providers disclose the hardware's impact-resilience bounds.</footnote>
We are interested in predicting the worst-case impulse to fulfill the hardware-affordable resilience bounds. In this paper, we analytically construct an impulse set in the three-dimensional space, whose interior covers all the candidate impulses that fulfill Coulumb's friction cone and task-space momentum conservation. Optimization-based controllers and in particular task-space formulated as a QP, can integrate the half-space-represented impulse set as part of the constraints set to steer the robot to intended impact objectives by iterative control.

## 2.2 Handling impacts in robotics

Shock-absorbing mechanisms can contribute to increasing the robot's resilience to impact and ease the control aspects, e.g., by providing more time for the controller to react to impacts. For example, Singh and Featherstone (2020) redesigned the quadruped leg's mechanical structure to cancel the shock propagating to the floating base. Battaglia et al. (2009) devised a soft cover with thickness and material that depends on (i) the required maximum pre-impact velocity, and (ii) on sensing plus field-bus transmission latencies. Dehio et al. (2022) attached soft pads on effectors, this allows the implementation of an impact-aware preview controller for dual-arm fast grabbing of objects. Recently, this has been extended by van Steen et al. (2023). In the bipedal walking domain, De Magistris et al. (2017); Pajon et al. (2017) mounted soft soles on humanoid feet for both leveraging impacts and casting terrains' irregularities, which allowed walking on gravels. Since hardware is not the focus of this study, note that there is an active research perspective that we do not cover on hybrid hardware/control codesign to increase the robots' resilience to impacts.

There are control strategies to avoid impacts. Pagilla and Yu (2001) modified the reference trajectories for zero reference velocity along the contact surface normal. The control design based on Zero-tilting Moment Point (ZMP) establishes contacts with near-zero velocity (e.g., Kajita et al. (2010)). Earlier, Grizzle et al. (2014) summarized that the impact-less reference trajectories are challenging to generate and inefficient to execute. Hypothetically, if we can leave aside the impact dynamics and assume that the associated state jumps remain hardware-affordable, Bombile and Billard (2022) proposed a sequence of time-invariant dynamical systems in a single control framework to continuously control the reach, grab, and toss motion; and even how to hit an object Khurana et al. (2021).

Impact plays a prominent role in robotic legged locomotion tasks, with the robot dynamics varying depending on the type and order of contact events. Complementarity problems uniformly represent a large range of (combinatorial) conditions that would require modeling using discontinuous functions, e.g., Stewart (2000). Combining the robot equations of motion and the complementarity conditions, the complementarity dynamical system (CDS) summarized in Hurmuzlu et al. (2004) offer a general conceptual tool to describe a walking robot's dynamics, see also Brogliato (2003). Posa et al. (2014) applied the complementarity Lagrangian models for trajectory optimization. However, as mentioned by Grizzle et al. (2014), effective control design for CDS was not established until at least 2014. Recently, novel control and

planning methods have emerged for CDS control issues, such as avoiding explicit contact sequencing by contact-implicit trajectory optimization (Manchester et al. 2019).

To avoid discontinuous task-space tracking errors, Yang and Posa (2021) projected the tracking objectives to a subspace invariant to the impact event, see more recent results in Yang and Posa (2023). The impact-invariant subspace is a generalization of a preliminary observation by Gong and Grizzle (2020), i.e., the angular momentum about the impact application point is invariant to the impacts.

Hu et al. (2007); Stanisic and Fernández (2012); Heck et al. (2016) model the dynamic contact transition with a mass-spring-damper approach. Yet, the continuity of the model cannot capture the non-smooth nature of impacts, thus failing to predict state jumps that could damage the hardware.

On the other hand, joint trajectory planners that are aware of the predicted state jumps can find the feasible impact motion. Konno et al. (2011) programmed the Fujitsu HOAP-2 humanoid robot to break wooden boards with a Karate chop. The impact motion consists of three-phase trajectories solved from nonlinear optimization. Konno et al. (2011) guaranteed standing stability by restricting the ZMP state jump. Rijnen et al. (2017) simulated an iCub robot impacting a wall employing two sets of pre-planned reference spreading trajectories (the first is pre-impact, the second is post impact) generated from an off-line task-space QP controller. Another controller switch robustly from the pre-impact trajectory to the post-impact one upon impact detection. Recently, van Steen et al. (2022) extended the reference spreading control to task-space trajectories. Konno et al. (2011); Rijnen et al. (2017); Shield et al. (2022); van Steen et al. (2022) share the same impulse prediction following Zheng and Hemami (1985).

Off-line trajectory planning is computationally expensive and non-reactive. Tracking the trajectories through contact-rich dynamics raises additional issues, e.g., impact timing, and invariant contact sequencing, see Johnson et al. (2016). Thus, the final state variation is discontinuous due to the perturbations. Pace and Burden (2017) summarized additional admissibility conditions that enforce continuous trajectory outcomes to be piecewise-differentiable with respect to the initial conditions. To prove the tracking error is infinitesimally contractive for mechanical systems, Burden et al. (2018) summarized conditions on the post-impact discretized system dynamics based on the stability index originally provided by Aizerman and Gantmacher (1958). Whereas, Burden et al. (2018) have not yet validated the contraction conditions in real robot experiments.

Within the scope of this paper, we focus on fulfilling the hardware-affordable resilience bounds before and right after a single impact event rather than generating periodic impacts. Hence, control design tools aiming for cyclic behaviors, e.g., Poincaré-map-dependent controllers (Grizzle et al. 2014) or reference spreading control (Rijnen et al. 2017), do not straightforwardly apply to our aim.

More than a decade ago, Zhang et al. (2004) and Abe et al. (2007) proposed task-space quadratic-programing (QP) based controllers for redundant manipulators and computer graphics. Task-space QP formulations are implemented for controlling redundant robots in recent years (Wang
et al. 2014; Righetti et al. 2013; Kuindersma et al. 2016; Liu et al. 2016; Bouyarmane and Kheddar 2018; Djeha et al. 2023). Like trajectory planners, the proposed impact-aware constraints embed the mapping between the impact-induced state jumps and the pre-impact contact velocities. The constraints modify the QP search space such that the optimized solution conforms to hardware limitations and post-impact feasibility, i.e., forward invariance.

# 3 Impact-Unaware QP-Control

In this section, we recall the task-space QP control formalism following Bouyarmane et al. (2019) and highlight impact-induced problems. Sec. 3.1 presents the equations of motion; Sec. 3.2 and 3.3 detail the impact-affected constraints for joint and centroidal spaces; Sec. 3.4 summarizes the QP controller without impacts.

## 3.1 Equations of motion

Given $n$ actuated joints, a robot has $(n+6)$ degrees of freedom (DOF); additional 6 DOF concern, if any, the floating-base configuration in $SE(3)$. Assuming $m_1$ sustained contacts, the equations of motion write:

$$M(q)\ddot{q} + N(q, \dot{q}) = B\tau + J^T W, \quad (1)$$

where $M(q) \in \mathbb{R}^{(n+6) \times (n+6)}$ is the generalized inertia matrix, $N(q, \dot{q}) \in \mathbb{R}^{(n+6)}$ gathers both Coriolis and gravitation vectors. We drop the dependency on $q$ and $\dot{q}$ in the rest of the paper for simplicity. $B \in \mathbb{R}^{(n+6) \times n}$ is the selection matrix for the actuated joints; $\tau \in \mathbb{R}^n$ is the joint torques vector. Furthermore, the matrix $J \in \mathbb{R}^{6m_1 \times (n+6)}$ and the vector $W \in \mathbb{R}^{6m_1 \times 1}$ vertically assemble $m_1$ Jacobians and wrenches of sustained contacts, respectively.

## 3.2 Joint space constraints

The equations of motion (1) fulfill the joint torque limits $\tau \le \tau_{t_k} \le \bar{\tau}$ by:

$$\begin{bmatrix} 1 \\ -1 \end{bmatrix} M \ddot{q}_{t_k} \le \begin{bmatrix} B\bar{\tau} \\ -B\bar{\tau} \end{bmatrix} + \begin{bmatrix} 1 \\ -1 \end{bmatrix} \left( J^T f_{t_{k-1}} - N_{t_{k-1}} \right). \quad (2)$$

Other constraints, such as collision avoidance, Cartesian space position constraints, actuated joint position or velocity limits... do not express directly in the robot state acceleration. We can rewrite them as a function of the QP decision variable $\ddot{q}_{t_k}$ by numerical derivation scheme, see examples in Djeha et al. (2020, 2023) for a closed-loop implementation.

## 3.3 Centroidal space constraints

In the humanoid case, Sugihara (2009) showed that the horizontal COM velocity $\dot{c}_{x,y} \in \mathbb{R}^2$ must remain within a convex 2D polygon $S_c$, i.e., the capture region, to keep the balance. $S_c$ depends on the sustained contacts and the COM position $c_{x,y} \in \mathbb{R}^2$. Employing proper matrix $G\dot{c}_{x,y}$ and vector $h\dot{c}_{x,y}$, the half-plane represented constraint $\dot{c}_{x,y} \in S_c$ writes:

$$G\dot{c}_{x,y} \dot{c}_{x,y} \le h\dot{c}_{x,y}. \quad (3)$$

Due to the kinematic and actuation limits, the robot controller shall minimize its angular momentum $L_c \in \mathbb{R}^3$

(Lee and Goswami 2012; Wiedebach et al. 2016). Suppose the angular momentum bounds $L_c \le \bar{L}_c$, we take the angular part $A_{\omega G}(q_{t_k}) \in \mathbb{R}^{3 \times (n+6)}$ from the centroidal momentum matrix $A_G(q)(q_{t_k})$ proposed by Orin et al. (2013) to formulate:

$$A_{\omega G}(q_{t_k}) \dot{q} \le \bar{L}_c. \quad (4)$$

These equations, being specific to humanoids, are the basis of a dedicated work (submitted elsewhere) on multi-legged balancing under impacts.

## 3.4 Impact-unaware whole-body QP controller

The optimization-based whole-body controller prioritizes multiple task objectives while imposing strict constraints:

$$\min_{\nu} \sum_{i \in \mathcal{X}_o} w_i \|e_i(\nu)\|^2$$

s.t. Joint space constraints:

Equations of motion and joint torque limits: (2),

Joint position, velocity, acceleration limits,

Centroidal space constraints:

COM velocity constraint: (3),

Angular momentum constraint: (4),

Other constraints, e.g., collision avoidance, vision... (5)

where $\sum_{i \in \mathcal{X}_o} w_i \|e_i(\nu)\|^2$ scalarizes multiple task objectives included in a set $\mathcal{X}_o$, $w_i$ weights the $i$-th task and $e_i(\nu)$ denotes the $i$-th task error. $e_i(\nu)$ is linear in terms of the QP decision variables $\nu$, see the details in Bouyarmane et al. (2019). The optimization variables are the generalized joint accelerations $\ddot{q}_{t_k}$ and the generating vectors of the contact wrench cone $f_{\lambda}(t_k)$, i.e., $\nu_{t_k} := \{\ddot{q}_{t_k}, f_{\lambda}(t_k)\}$. The QP updates the robot state and sensory feedback data in each control cycle.

When it occurs, impact instantaneously changes the contact velocity. As a consequence, all the constraints that are expressed in terms of velocity, e.g., joint velocity bounds, angular momentum, COM velocity... will have their value changed with a substantially high increment. As a result, the QP search space might shrink instantly to an empty set, rendering the QP infeasible for the subsequent control iteration.

Our main idea, is to prepare such constraints to handle and be robust to such an abrupt change of the velocity values. If possible, without changing much the QP control formalism and framework. The controller will then have the dual objective to meet at best desired possible pre-impact speeds while keeping the QP controller feasible for the next iteration (i.e., forward invariant). That is to say, make the QP aware of the possible change in the tasks and constraints due to post-impact speed. We therefore explicitly need models that constrain the post-impact state with an analytically-computed polyhedra in Sec. 4 and formulate impact-aware constraints in Sec. 5.

# 4 Modeling 3D Frictional Impact

Sec. 4.1 presents the impulse set that fulfills Coulomb's friction cone and the task-space momentum conservation.
The QP can integrate this set with an additional constraint. To numerically represent the impulse set, it's common to discretize Coulomb's friction cone as a polyhedron discretized-cone. Projecting the cone's vertices onto different spaces, we constrain other quantities: (1) contact velocities in Sec. 4.2, (2) joint velocities in Sec. 4.4.1, and (3) COM velocities and angular momentum in Sec. 4.5. Utilizing heuristics, we also constrain the peak impulsive force in Sec. 4.3 and the joint torque jump in 4.4.2.

We use the coordinate frames illustrated in Fig. 2 to define and exploit the impact model. Sec. 4.6 introduces the (1) impulse-to-velocity mapping in task space, and (2) the transforms that can uniformly represent the velocities and impulses in the same frame.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** The red, green and blue color indicates the x, y, and z axis respectively. The inertial frame $\mathcal{F}_O$ and the centroidal frame $\mathcal{F}_c$ share the same orientation. When the robot impacts, the force torque sensor, the end-effector frame $\mathcal{F}_e$, and the contact point frame $\mathcal{F}_p$ share the same translation $p_{Op} \in \mathbb{R}^3$ with different orientations.</figcaption></figure>

The whole-body QP can seamlessly integrate the polyhedra to enable impact awareness since all the polyhedra are linear with respect to the pre-impact contact velocity and, hence, to the QP optimization variables (joint accelerations).

## 4.1 Feasible impulse set

We assume the following statements are true during the impact event:

A.1 The impact force is dominant compared to other forces. Consequently, generalized forces such as the centrifugal forces or motor command torques are negligible (Stronge 2000, Chapter 8.1.1).

A.2 The impulse $\iota \in \mathbb{R}^3$ fulfills Coulomb's friction cone (Keller 1986), (Stronge 2000, Eq. 4.8), (Jia and Wang 2017):

$$ \sqrt{\iota_x^2 + \iota_y^2} \le \mu \iota_z. \qquad (6) $$

A.3 During the impact, the joint position changes are considered negligible, e.g., Zheng and Hemami (1985); Grizzle et al. (2014); Siciliano and Khatib (2016); Konno et al. (2011). Yet, it is well-known that substantial changes are expected in the joint velocities; see recent experiments for kinematic-controlled robots in Wang et al. (2022b).

A.4 The coefficient of restitution $c_r \in \mathbb{R}^+$ is configuration-dependent (Stoianovici and Hurmuzlu 1996; Wang et al.

2022b). We can not apply the material-dependent $c_r$ for a kinematic-controlled robot (Wang et al. 2022b, Fig. 5).

A.5 The contact area is tiny compared to the robot's dimension, hence, a 3D point contact model is appropriate (Chatterjee and Ruina 1998).

A.6 The impact induces significant impulsive forces and negligible moments (Stronge 2000; Chatterjee and Ruina 1998).

A.7 The impacting limb has a minimum of three degrees of freedom. The joint configuration is not singular.

The feasible impulse set is the intersection of:

1. the discretized Coulomb's friction cone, see Sec. 4.1.1;

2. the two planes of restitution given the coefficient of restitution $c_r \in [c_r, \bar{c}_r]$, see Sec. 4.1.2.

We analytically represent the impulse set as a polyhedron, whose interior points are candidate impulses fulfilling assumptions A.1–A.7.

### 4.1.1 Coulomb's friction cone
The impulse fulfills Coulomb's law of friction (Jia et al. 2019; Stronge 2000). We discretize the friction cone (in A.2) with $N_\mu$ vertices:

$$ \mathcal{X}_\mu := \{\iota \in \mathbb{R}^3, \lambda_\mu \ge 0 : \iota = K_\mu \lambda_\mu\}, \qquad (7) $$

where the matrix $K_\mu \in \mathbb{R}^{3 \times N_\mu}$ is constant and $\lambda_\mu \in \mathbb{R}^{N_\mu}$ are the impulse generators. Note that the same discretization applies to the static friction cone (Abe et al. 2007; Bouyarmane et al. 2019).

### 4.1.2 Planes of restitution
During the impact event, the contact velocity $v \in \mathbb{R}^3$ is determined by the sum of the pre-impact velocity $v^- \in \mathbb{R}^3$ and the net contact velocity jump $\Delta v \in \mathbb{R}^3$: $v = v^- + \Delta v$. The inverse inertia matrix (IIM) $W \in \mathbb{R}^{3 \times 3}$ defines the task-space impulse-to-velocity mapping $\Delta v = W\iota$ (Stronge 2000, Chapter 4.1) (Jia and Wang 2017, Sec. 2.1). Thus, the contact velocity writes:

$$ v = v^{-} + \Delta v = v^{-} + W\iota. \qquad (8) $$

While Stronge (2000) and Jia and Wang (2017) discuss impacts between pairs of rigid bodies, they do not provide a specific definition for the inverse inertia matrix ($W$) in the context of articulated robots. To bridge this gap, we introduce the $W$ in Sec. 4.6.2, based on assumption A.3.

The impact event consists of the compression and restitution phases. The compression ends when the normal contact velocity is zero, that is $v_n = \hat{n}^\top v = 0$, where $\hat{n} \in \mathbb{R}^3$ denotes the impact normal (often taken as the impacted surface normal). Equating the inner product of $\hat{n}$ and (8) to zero:

$$ v_n = v_n^{-} + \hat{n}^{\top} W \iota = 0, $$

we obtain the plane of compression in the 3D impulse space:

$$ \mathcal{X}_{\text{comp}} := \{\iota \in \mathbb{R}^3 : v_n^{-} + \hat{n}^{\top} W \iota = 0\}. \qquad (9) $$

The coefficient of restitution, $c_r$, is known to vary with the robot configuration A.4. This variability exists not only for a single rigid body with constant inertia, where $c_r$ differs
across different contact points (Stoianovici and Hurmuzlu 1996), but also for articulated robots (Wang et al. 2022a). Consequently, researchers have adopted different approaches to address this issue. Some have chosen to consider inelastic impacts (Rijnen et al. 2017). In our case, we empirically bound $c_r$ within the range [$c_r, \bar{c}_r$] based on our previous experiments (Wang et al. 2022a). This allows us to define two constraining planes similar to (9):

$$-c_r v_n \le v_n = v_n + \hat{n}^T W \mathbf{I} \le -\bar{c}_r v_n. \quad (10)$$

The negative sign arises because $v_n$ in the contact point frame $\mathcal{F}_p$ undergoes a transition from negative to positive at the end of restitution, see the details in Sec. 4.6.1. This convention follows the coordinate frame definitions established in impact mechanics, as outlined by Stronge (2000). Thus, (10) defines the set of impulses as:

$$\mathcal{X}_{\text{res}} := \{ \mathbf{I} \in \mathbb{R}^3 : -c_r v_n \le v_n + \hat{n}^T W \mathbf{I} \le -\bar{c}_r v_n \}. \quad (11)$$

### 4.1.3 Impulse set
The intersection of the friction cone (7) and the two planes (10-11) defines the impulse polyhedron:

$$\mathcal{X}_{\mathbf{I}} := \{ \mathbf{I} \in \mathbb{R}^3 : \mathbf{I} \in \mathcal{X}_{\text{res}} \cap \mathcal{X}_{\mathbf{I}} \}. \quad (12)$$

<figure><img src="image_5.png" /><figcaption>**Figure 3.** We visualize the impulse set $\mathcal{X}_\mathbf{I}$ (12) when the Panda robot in Fig. 2 impacts a rigid surface with contact velocity $v = [0.004, 0.054, -0.175]$ m/s (in the contact point frame). $\mathcal{X}_\mathbf{I}$ contains the interior of the intersection by the plane of restitutions (10-11) and Coulomb's friction cone (7).</figcaption></figure>

## 4.2 Contact velocity

Given the impulse set and the impulse-to-velocity mapping $W$, the contact velocity set is:

$$\mathcal{X}_v := \{ \Delta v \in \mathbb{R}^3 : \Delta v = W \mathbf{I}, \mathbf{I} \in \mathcal{X}_\mathbf{I} \}. \quad (13)$$

## 4.3 Contact force

We suppose that the peak contact force is bounded by:

$$f \le a \frac{\mathbf{I}}{\delta t},$$

where we estimate the impact duration $\delta t$ and the positive scalar $a > 0$ from experiment data. Given the impulse set (12), the peak contact force set is:

$$\mathcal{X}_f := \{ f \in \mathbb{R}^3 : f \le a \frac{\mathbf{I}}{\delta t}, \mathbf{I} \in \mathcal{X}_\mathbf{I} \}. \quad (14)$$

## 4.4 Joint-space projections

### 4.4.1 Joint velocities
According to the recently-proposed impulse-to-joint-velocity mapping by Wang et al. (2022b):

$$\Delta \dot{\mathbf{q}} = J^\top (JJ^\top)^{-1} W \mathbf{I} = \arg \min_{\Delta \dot{\mathbf{q}}} \|J \Delta \dot{\mathbf{q}} - W \mathbf{I}\|^2,$$

the joint velocity jumps span the following set:

$$\mathcal{X}_{\Delta \dot{\mathbf{q}}} := \{ \Delta \dot{\mathbf{q}} \in \mathbb{R}^n : \Delta \dot{\mathbf{q}} = J^\top (JJ^\top)^{-1} W \mathbf{I}, \mathbf{I} \in \mathcal{X}_\mathbf{I} \}. \quad (15)$$

### 4.4.2 Impulsive joint torques
Given the peak contact force set (14) and the mapping $J^\top \Delta \mathbf{f}$, the impulsive joint torque jumps are:

$$\mathcal{X}_{\Delta \gamma} := \{ \Delta \gamma \in \mathbb{R}^n : \Delta \gamma = J^\top \mathbf{f}, \mathbf{f} \in \mathcal{X}_f \}. \quad (16)$$

**Remark 4.1.** Impulses propagate through the mechanical linkages and have the potential to affect various components of the robot, including joint mechatronic implements such as gearboxes, torque sensors (if any), mechanical linkages, and covers. It is crucial to mitigate impacts to prevent damage to these components.

The exact origin of the limitation is not of utmost importance, because, ultimately, a conservative approach is taken by considering the minimum or maximum allowable shock. The proposed QP controller embeds such constraints concurrently to address intentional impacts. Stall torque limits are established based on the robot's equations of motion. These limits help ensure that a robot operates within safe operating conditions. The impact-aware constraint (30) imposes a distinct tolerance on the impulsive torque limits, allowing for control over the magnitude of the impact forces exerted on the robot. ■

## 4.5 Centroidal space projections

High-stiffness kinematics controlled robots behave as a composite rigid body during the impact event (Wang et al. 2022b). Thus, we have the COM velocity jumps $\mathcal{X}_{\Delta \dot{\mathbf{e}}}$ and angular momentum jumps $\Delta \mathcal{L}$. We ease the reading by leaving the impulse to COM velocity mapping (21) and the impulse to angular momentum mapping (22) in Sec. 4.6.

### 4.5.1 COM velocity
Given the impulse polyhedron (12) and the mapping (21), the COM velocity jump span the following set:

$$\mathcal{X}_{\Delta \dot{\mathbf{e}}} := \{ \Delta \dot{\mathbf{e}} \in \mathbb{R}^3 : \begin{bmatrix} \Delta \dot{\mathbf{e}} = \frac{1}{m} R_{cp} \mathbf{I}, \\ \mathbf{I} \in \mathcal{X}_\mathbf{I} \end{bmatrix} \}, \quad (17)$$

where $m \in \mathbb{R}^+$ denotes the robot mass and $R_{cp}$ denotes the rotation from COM to the contact point.
**4.5.2 Angular momentum** Similar to the derivation of (17), we have the set of angular momentum jump given the mapping (22):

$$ \chi_{\Delta L} := \{\Delta L \in \mathbb{R}^3 : \Delta L = \hat{p}_{cp} R_{cp} \iota, \text{ and } \iota \in \chi_{\iota}, \} \quad (18) $$

where $\hat{p}_{cp} \in \mathbb{R}^{3 \times 3}$ denotes the skew-symmetric matrix constructed with the translation: $p_{cp} = p - c$.

## 4.6 Implementation details

Sec. 4.6.1 summarizes the coordinate frame definitions, aligning with the notation used in impact mechanics (Stronge 2000; Jia and Wang 2017). Sec. 4.6.2 presents the equivalent inertial properties at the end-effector.

**4.6.1 The coordinate frame specifications** We express the contact velocity and impulse in the contact point frame $\mathcal{F}_p$, as depicted in Figure 2. As viewed in frame $\mathcal{F}_p$, the pre-impact normal contact velocity is negative $v_n^- < 0$ and the impulse is positive $\iota_n > 0$. Following the geometric notations established by Murray et al. (1994) and our earlier papers on contact modeling (Wang et al. 2022b,a), we specify the contact velocity as:

$$ v_{Oe}^p \in \mathbb{R}^3 $$

which corresponds to the end-effector velocity with respect to the inertial frame $\mathcal{F}_O$, but represented (viewed) in the contact point frame $\mathcal{F}_p$. We compute the associated Jacobian according to (1) the rotation from the end-effector frame $\mathcal{F}_e$ to the contact point frame $\mathcal{F}_p$: $R_{ep}^\top$, and (2) the translation rows of the body velocity Jacobian $^tJ_{Oe}^b \in \mathbb{R}^{3 \times (n+6)}$ as:

$$ ^tJ_{Oe}^p = R_{ep}^\top {}^tJ_{Oe}^b. \quad (19) $$

The superscript $^b$ indicates the body velocity (Murray et al. 1994, Page 55), i.e., velocity of the origin of the body coordinate frame relative to the inertial frame, as observed in the current body frame.

In line with Murray et al. (1994), we represent a six-dimensional velocity as $V_{Oe}^b = [v_{Oe}^b \ w_{Oe}^b]$, which stacks the linear velocity $v_{Oe}^b \in \mathbb{R}^3$ on top of the angular velocity $w_{Oe}^b \in \mathbb{R}^3$.

The adjoint transformation, as described in Murray et al. (1994, Eq. 2.57), allows the mapping of a six-dimensional velocity or a wrench (or momentum, impulse) from one coordinate frame ($\mathcal{F}_e$) to another ($\mathcal{F}_p$) as follows:

* velocity

$$ Ad_{g_{cp}}^{-1} = \begin{bmatrix} R_{cp}^{\top} & -R_{cp}^{\top} \hat{p}_{cp} \\ 0 & R_{cp}^{\top} \end{bmatrix} = \begin{bmatrix} R_{cp}^{\top} & 0 \\ 0 & R_{cp}^{\top} \end{bmatrix}, $$

* wrench (or momentum, impulse):

$$ Ad_{g_{pe}^{-1}}^{\top} = \begin{bmatrix} R_{pe} & 0 \\ \hat{p}_{pe} R_{pe} & R_{pe} \end{bmatrix} = \begin{bmatrix} R_{pe} & 0 \\ 0 & R_{pe} \end{bmatrix}. $$

It is important to note that $p_{ep} = 0$ because we assume $\mathcal{F}_p$ shares the same origin with $\mathcal{F}_e$ when the impact occurs. Therefore, the translation-dependent term $-R_{cp}^\top \hat{p}_{cp}$ disappears in (19). In the rest of Sec. 4.6, we use unique colors to distinguish body velocity $V_{Oe}^b$, contact velocity $v_{Op}^p$, and adjoint maps $Ad_{g_{cp}}^{-1}$.

**4.6.2 The inverse inertia matrix** The inverse inertia matrix $W \in \mathbb{R}^{3 \times 3}$ denotes the impulse-to-velocity mapping at the contact point. We derive $W$ with the composite-rigid-body inertia $crbI \in \mathbb{R}^{6 \times 6}$ under assumption A.3: The momentum conservation in the centroidal frame $\mathcal{F}_c$ involves negligible impulsive moment and other forces according to assumptions A.6 and A.1:

$$ crbI \Delta V_{Oe}^b = Ad_{g_{cp}^{-1}}^{\top} \begin{bmatrix} \iota \\ 0 \end{bmatrix} = \begin{bmatrix} R_{cp} & 0 \\ \hat{p}_{cp} R_{cp} & R_{cp} \end{bmatrix} \begin{bmatrix} \iota \\ 0 \end{bmatrix}. \quad (20) $$

The adjoint map $Ad_{g_{cp}^{-1}}^{\top} \in \mathbb{R}^{6 \times 6}$ transforms the momentum (or wrench) from frame $\mathcal{F}_p$ to $\mathcal{F}_e$. The body velocity $V_{Oe}^b \in \mathbb{R}^6$ is the same as the average velocity (Orin and Goswami 2008, Eq. 24). The composite-rigid-body inertia $crbI$ is block-diagonal:

$$ crbI = \begin{bmatrix} m \mathbb{1}_{3 \times 3} & 0 \\ 0 & \mathcal{I} \end{bmatrix}, $$

with $m$ represents the total mass, $\mathbb{1}_{3 \times 3} \in \mathbb{R}^{3 \times 3}$ is the identity matrix, and $\mathcal{I} \in \mathbb{R}^{3 \times 3}$ denotes the rotational inertia. Consequently, the first three rows of (20) define the COM velocity jump:

$$ \Delta v_{Oe}^b = \Delta \dot{c} = \frac{1}{m} R_{cp} \iota, \quad (21) $$

while the remaining rows of (20) define the jump in angular momentum:

$$ \Delta \mathcal{L} = \hat{p}_{cp} R_{cp} \iota. \quad (22) $$

We compute the velocity jump $\Delta V_{Oe}^p \in \mathbb{R}^6$ by left multiplying $crbI^{-1}$ and the velocity transform $Ad_{g_{cp}^{-1}}^{-1}$: to (20):

$$ \Delta V_{Oe}^p = Ad_{g_{cp}}^{-1} \Delta V_{Oe}^b = Ad_{g_{cp}}^{-1} crbI^{-1} Ad_{g_{cp}^{-1}}^{\top} \begin{bmatrix} \iota \\ 0 \end{bmatrix}. \quad (23) $$

Thus, the inverse inertia matrix $W$ is the $3 \times 3$ upper-left corner of $Ad_{g_{cp}^{-1}}^{-1} crbI^{-1} Ad_{g_{cp}^{-1}}^{\top}$:

$$ \Delta v_{Oe}^p = \underbrace{(\frac{\mathbb{1}_{3 \times 3}}{m} - R_{cp}^\top \hat{p}_{cp} \mathcal{I}^{-1} \hat{p}_{cp} R_{cp})}_{W} \iota. \quad (24) $$

Alternatively, we can find a more detailed derivation of $W$ in Wang et al. (2022a, Sec. IV.B).

# 5 Impact-Aware Control Design

In this section, Sec. 5.1 presents the impact-aware template for constraints; Sec. 5.2 substitutes the impulse polyhedra and robot states of interest into the template constraint.

Due to the modified search space, the impact-aware QP summarized in Sec. 5.3 is feasible (under our hypotheses and correct values of the bounds) if the impact occurs in the next control cycle, i.e., all the post-impact states fulfill hardware resilience bounds and the QP is robust to abrupt changes of the velocity.
## 5.1 Impact-aware constraints

We account for the post-impact state of constraints that are written in a generic quantity $d$ as follows

$$D d_{k+1}^{+} \leq \bar{d} \qquad (25)$$

with the matrix $D$ representing half-spaces and the vector $\bar{d}$ representing the upper and lower bounds. Substituting $\bar{d}_{k+1} \approx d_{t_k}$ and the prediction $d_{k+1}^{+} = d_{k+1}^{-} + \Delta d_{t_{k+1}}$, we obtain

$$D(d_{t_k} + \Delta d_{t_{k+1}}) \leq \bar{d}. \qquad (26)$$

We define a special Jacobian $\mathcal{J}_{\Delta d}$ as a mapping from the impulse to the state jump of different quantities:

$$\Delta d_{t_{k+1}} = \mathcal{J}_{\Delta d} \iota. \qquad (27)$$

The various special Jacobians are:

* joint velocities $\mathcal{X}_{\Delta \dot{q}}$ given (15):

$$\mathcal{J}_{\Delta \dot{q}} = J^{\top} (J J^{\top})^{-1} W;$$

* impulsive joint torques $\mathcal{X}_{\Delta \gamma}$ given (14) and (16):

$$\mathcal{J}_{\Delta \gamma} = J^{\top} a \frac{1}{\delta t};$$

* COM velocity $\mathcal{X}_{\Delta \dot{c}}$ given (17):

$$\mathcal{J}_{\Delta \dot{c}} = \frac{1}{m} R_{cp};$$

* centroidal angular momentum $\mathcal{X}_{\Delta \mathcal{L}}$ given (18):

$$\mathcal{J}_{\Delta \mathcal{L}} = \hat{p}_{cp} R_{cp}.$$

Thus, we reformulate (26) by substituting the impact-induced jump $\Delta d_{t_{k+1}}$ (27):

$$D \mathcal{J}_{\Delta d} \iota \leq \bar{d} - D d_{t_k}. \qquad (28)$$

Note that it is easy to adapt (28) to constrain the actuated joints $q$ only.

## 5.2 Constraining the post-impact states

The template constraint (28) adapts to different quantities by choosing the corresponding half-planes $D$ and the bounds $\bar{d}$. We enumerate the following:

Constraint 1: Joint velocity $d := \dot{q}$ with $D := [1, -1]^T$ and $\bar{d} := \dot{\bar{q}}$ yields:

$$\begin{bmatrix} 1 \\ -1 \end{bmatrix} \mathcal{J}_{\Delta \dot{q}} \iota \leq \begin{bmatrix} \bar{q} \\ -\dot{q} \end{bmatrix} - \begin{bmatrix} 1 \\ -1 \end{bmatrix} \dot{q}_{t_k}, \qquad (29)$$

where $\mathcal{J}_{\Delta \dot{q}}$ only includes the actuated joint rows.

Constraint 2: Impulsive joint torque $d := \Delta \gamma$ with $D := [1, -1]^T$ and $\bar{d} := \underline{\Delta \gamma}$ leads to:

$$\begin{bmatrix} 1 \\ -1 \end{bmatrix} \mathcal{J}_{\Delta \gamma} \iota \leq \begin{bmatrix} \underline{\Delta \gamma} \\ -\underline{\Delta \gamma} \end{bmatrix}, \qquad (30)$$

where $\Delta \gamma_{t_k} = 0$ since we assume the impact occurs at the time step $k+1$.

Constraint 3: Angular momentum $d := \mathcal{L}_c$ with $D := 1$ and $\bar{d} := \mathcal{L}_c$ yields:

$$\mathcal{J}_{\Delta \mathcal{L}_c} \iota \leq \mathcal{L}_c - \mathcal{L}_{c_{t_k}}. \qquad (31)$$

Constraint 4: Horizontal COM velocity $d := \dot{c}_{x,y}$ with $D := \mathcal{G}_{\dot{c}_{x,y}}$ and $\bar{d} := \mathbf{h}_{\dot{c}_{x,y}}$ generates:

$$\mathcal{G}_{\dot{c}_{x,y}} \mathcal{J}_{\Delta \dot{c}_{x,y}} \iota \leq \mathbf{h}_{\dot{c}_{x,y}} - \mathcal{G}_{\dot{c}_{x,y}} \dot{c}_{x,y}(t_k). \qquad (32)$$

## 5.3 Impact-aware whole-body QP controller

To seamlessly integrate the impact-aware constraints (29-32) with respect to the impulse, we need to:

1. add the impulse generators $\lambda_{\mu} \in \mathbb{R}^{N_{\mu}}$ to the QP decision variables $\nu$;

2. rewrite the pre-impact normal contact velocity $v_n^{-}$ by the current-step joint velocities $\dot{q}_{t_k}$ and accelerations $\ddot{q}_{t_{k+1}}$ (QP decision variable).

Assuming the impact occurs at the next control iteration, we approximate the pre-impact normal velocity as:

$$\begin{aligned} v_n^{-}(t_{k+1}) &= v_n(t_{k+1}) = J_{t_{k+1}} \dot{q}_{t_{k+1}} \\ &= (J_{t_k} + \dot{J}_{t_k} \Delta t + \ddot{J}_{t_k} \Delta t^2 + \dots)(\dot{q}_{t_k} + \ddot{q}_{t_k} \Delta t) \\ &= J_{t_k} \dot{q}_{t_k} + J_{t_k} \ddot{q}_{t_k} \Delta t + \underbrace{\dot{J}_{t_k} \dot{q}_{t_k} \Delta t + \dot{J}_{t_k} \ddot{q}_{t_k} \Delta t^2 + \dots}_{\approx 0}, \end{aligned}$$

where we read $\dot{q}_{t_k}$ from the current robot state, $\Delta t$ denotes the sampling period, and we define the Jacobian as: $J_{t_k} = \hat{n}^{\top} t J_{Oe}^p$ given (19).

Substituting $v_n^{-}(t_{k+1}) = J_{t_k} \dot{q}_{t_k} + J_{t_k} \ddot{q}_{t_k} \Delta t$ into (11), the following constraints specify the impulse polyhedron:

Friction cone: (7),

$$\lambda_{\mu} \geq 0;$$

Planes of restitution: (11),

$$\begin{aligned} &-(\bar{c}_r + 1)J_{t_k} \ddot{q}_{t_k} \Delta t - \hat{z}^{\top} W K_{\mu} \lambda_{\mu} \geq (\bar{c}_r + 1)J_{t_k} \dot{q}_{t_k}, \\ &-(c_r + 1)J_{t_k} \ddot{q}_{t_k} \Delta t - \hat{z}^{\top} W K_{\mu} \lambda_{\mu} \leq (c_r + 1)J_{t_k} \dot{q}_{t_k}. \end{aligned} \qquad (33)$$

Thus, we enhance the impact-unaware optimization-based controller (5) with the impulse polyhedron $\mathcal{X}_{\iota}$ (33) and the impact-aware constraints (29-32) (substituting $\iota = K_{\mu} \lambda_{\mu}$ from (7)); hence, enforcing consistency between the pre-impact state and the impulse, that is:

$$\min_{\nu, \lambda_{\mu}} \sum_{i \in \mathcal{X}_o} w_i \|e_i(\nu)\|^2$$

s.t. Impulse Polyhedron $\mathcal{X}_{\iota}$ given by (33),

Joint space constraints:

Post-impact joint velocity: (29),

Post-impact impulsive joint torque: (30),

Joint position, velocity, and torque (2),

Centroidal space constraints:

Post-impact angular momentum: (31),

Post-impact COM velocity: (32),

Other constraints, e.g., collision avoidance: (26). (34)
The improved formulation (34) shares the same objectives as the formulation (5) and modifies its search space such that all the post-impact states fulfill the corresponding bounds. Hence, (34) can decide on the highest feasible contact velocities if required.

**Remark 5.1.** We deactivate the impact-awareness upon impact detection by switching the QP controller from (34) to (5). Since the constraints associated with (34) are more conservative than (5), switching from (34) to (5) by full initialization does not lead to infeasible solutions.

Integrating (i.e., activating) impact-aware constraints by switching from (5) to (34) may lead to infeasible QP due to the actuation limits. Yet, this problem is common to any other usual inequality constraint that is triggered or inserted into the QP on the fly. Solutions exist for the case of one constraint, see e.g., Del Prete (2018); Djeha et al. (2020). Yet interaction among several constraints is still an open problem that we are currently investigating. ■

# 6 Experiments

We first assess our impact-aware control design with the Panda manipulator. All the manipulator's post-impact states fulfilled the proposed impact model. The impact-aware QP (34) exploited the manipulator's resilience bounds (the parameters are provided by our partner Franka Emika) to achieve the highest possible impact velocity; see the video‡.

The impact model is suitable for both the Panda and the humanoid robot HRP-4 since they are both stiff kinematic-controlled robots. Enforcing impact-aware constraints facilitated robot motions to become impact-friendly, thereby enabling the HRP-4 to grasp boxes with human-levele swift motions. The detailed experiment description and the highlights are:

**Experiment 1** The Panda manipulator impacted an ATI sensor with two distinct configurations. In configuration one (refer to Fig. 5(a)), the friction coefficient is 0.279, while in configuration two (refer to Fig. 6(a)), the friction coefficient is lower with a value of 0.114 due to a different contact surface material. Comparing the impulse set of configuration one in Fig. 3, the reduction in friction coefficient resulted in a narrower impulse set, as illustrated in Fig. 4. In both cases, the robot promptly pulled back the end-effector upon impact detection. To mitigate random effects and ensure consistent observations, we conducted the experiment 10 times for each configuration.

(H-1.1) All the impulses (measured by integrating the ATI force readings) fulfilled the impulse polyhedron proposed in Sec. 4.

(H-1.2) The impact-aware joint constraints are fulfilled.

(H-1.3) We empirically concluded that when the impulsive joint torque constraint (30) is active, the post-impact joint velocities are still away from the bounds.

(H-1.4) Given a high reference normal contact velocity, i.e., 0.50 m/s, the QP autonomously steered the robot to the highest impact-aware feasible contact velocity in real-time.

**Experiment 2** The HRP-4 robot swiftly grabbed two different boxes, see the snapshot in Fig. 10. The two arm's impacts are synchronized and symmetric.

(H-2.1) Grabbing the box swiftly without marking any stop or slow at contacts.

(H-2.2) The high contact velocities, i.e., 0.25 m/s, were well-tracked.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** We visualize the impulse set $\mathcal{X}_t$ when the Panda manipulator in Fig. 6(a) impacted at contact velocity $v = [-0.0001, -0.0117, 0.174]$ m/s (in the contact point frame). In comparison to Fig. 3, the reduction in friction coefficient from 0.279 to 0.114 leads to a narrower impulse set $\mathcal{X}_t$.</figcaption></figure>

## 6.1 Panda Manipulator Experiments

Sec. 6.1.1 details the steps to measure the post-impact data; Sec. 6.1.2 customizes the impact-aware QP for the Panda manipulator; and Sec. 6.1.3 analyzes the results.

### 6.1.1 Data acquisition
We applied the open-source dynamics model§ for the 7 DOF Panda manipulator from Franka Emika. It is worth mentioning that alternative models are also available¶.

We mounted a 3D printed semi-spherical end-effector on the Panda robot to meet the point-contact assumption A.5, see Fig. 2. The friction coefficient is $\mu = 0.279$ for configuration one (Fig. 5(a)) and $\mu = 0.114$ for the second configuration (Fig. 6(a)).

The impact event lasted for about 20 ms according to the recent experimental study we performed in Wang et al. (2022a). Subsequently, we measured the joint velocity

<footnote>‡ https://youtu.be/78xPQ_7qM4I</footnote>

<footnote>§ github.com/jrl-umi3218/mc_panda/blob/master/data/urdf/panda_default.urdf</footnote>

<footnote>¶ github.com/marcocognetti/FrankaEmikaPandaDynModel</footnote>

<footnote>|| github.com/StanfordASL/PandaRobot.jl</footnote>
changes in an interval (5 ms before and 15 ms after the impact) to be the joint velocity jumps $\Delta \dot{q}^o$.

We sampled the ATI-mini45 sensor readings at 25000 Hz without low-pass filtering and integrated the force measurements to compute the impulse in three dimensions.

**Coefficient of restitution** We construct the two planes of restitution (11) by choosing the upper and lower bound of the coefficient of restitution $c_r \in [c_r, \bar{c}_r]$.

According to the observations from hundreds of impact experiments (Wang et al. 2022a, Sec. VII.C), the kinematic-controlled robot impact is almost inelastic (the coefficient of restitution is close to zero $c_r \le 0.15$) if the contact velocity is greater than 0.1 m/s. We choose the upper bound $\bar{c}_r = 0.3$, conservatively higher than 0.15, and the lower bound $\bar{c}_r = 0.0$, which corresponds to inelastic impact.

**Impact detection** The virtual force sensor of the Panda robot cannot capture the impact dynamics. Timely impact detection is essential for accurately observing the post-impact states. We achieved 3 – 6 ms detection time by thresholding the joint torque error:

$$e_{\tau} = \sum_{i=5}^{6} |\tau_{i} - \tau_{i}^{\text{ref}}| \le \tau_{0},$$

where the motor torque $\tau_i$ and its reference $\tau_i^{\text{ref}}$ update at 1 ms. The threshold $\tau_0$ is 2.5 N·m.

**6.1.2 Controller formulation** A manipulator requires the following set of impact-aware constraints:

$$ \begin{array}{l} \text{Impulse polyhedron } \mathcal{X}_t \text{ (33),} \\ \text{Joint space constraints:} \\ \text{Post-impact joint velocity: (29),} \\ \text{Post-impact impulsive joint torque: (30).} \end{array} \qquad (35) $$

Thus, we customized the impact-aware controller (34) for the Panda manipulator as:

$$ \begin{array}{ll} \min_{\boldsymbol{v}, \boldsymbol{\lambda}_{\mu}} & \sum_{i \in \mathcal{X}_o} w_i \|e_i(\boldsymbol{v})\|^2 \\ \text{s.t.} & \text{Manipulator Impact-awareness (35),} \end{array} \qquad (36) $$

Joint position, velocity, and torque (2).

The impulsive joint torque constraint (30) depends on the peak contact force. According to the contact force profiles (Wang et al. 2022a, Fig. 6), we conservatively choose the impact duration $\delta t = 18$ ms and the positive scalar $a = 3$ to construct the peak contact force set $\mathcal{X}_f$ (14).

We assigned an unrealistically-high reference contact velocity, i.e., $v^{\text{ref}} = [0.5, 0, 0]$ m/s (represented in the inertial frame $\mathcal{F}_O$ in Fig. 2). The other tasks include the end-effector orientation and posture tasks; see the detailed QP formulation by Bouyarmane et al. (2019).

### 6.1.3 Result analysis

**The contact velocities** The impact-aware QP (36) solves the desired joint commands at each control iteration. Since the reference is too high to be precisely tracked, (36) modified the maximum feasible contact velocity as shown in Fig. 5(b) and Fig. 6(b) for the two configurations (H-1.4).

The light area in Fig. 5(b) and Fig. 6(b) (H-1.4) highlight the period during which the impact-aware constraints (35) modify the optimizer's search space. The active constraint is the impulsive joint torque constraint (30), i.e., the QP solution in Fig. 8 is close to the upper bound of 12 N·m. Due to the conservative choice of parameters, the measured impulsive joint torque $\Delta \gamma = J^\top a \frac{\partial \dot{q}}{\partial t}$ (16) of the sixth joint was lower than the corresponding QP solution.

**Fulfilling the Impulse Polyhedra** Fig. 5(c) and Fig. 6(c) illustrate that all measured impulses (from 2 × 10 experiments marked by red crosses) $\iota^o$ lie within the (transparent-green) polyhedron-represented impulse set $\mathcal{X}_t$ (12) (H-1.1).

Substituting $\iota^o$ into (14), (17), and (18), Fig. 5(d)-5(f) and Fig. 6(d)-6(f) illustrate that all the quantities are constrained by the corresponding polyhedron-represented sets $\mathcal{X}_f$ (14), $\mathcal{X}_{\Delta \dot{q}}$ (17), $\mathcal{X}_{\Delta t}$ (18).

**Fulfilling the impact aware constraints** We focus on the first configuration, see Fig. 5(a), as similar results are obtained for the second. According to the Panda robot descriptions, the joint velocity bounds are:

$$ \bar{\dot{q}} = \pm [2.175, 2.175, 2.175, 2.175, 2.61, 2.61, 2.61] \text{ rad/s} \qquad (37) $$

Fig. 7 plots the measured joint velocity jumps $\Delta \dot{q}^o$ from ten experiments. All the measurements are within the set $\mathcal{X}_{\Delta \dot{q}}$ (15) (H-1.2). The predicted joint velocity jumps are away from the bounds (37), i.e., the constraints (29) are not active (H-1.3).

We define the impulsive joint torque bounds:

$$ \overline{\Delta \gamma} = \pm [87, 87, 87, 87, 12, 12, 12] \text{ N·m.} $$

Substituting the measured impulse $\iota^o$, the predicted impulsive joint torque $\Delta \gamma = J^\top a \frac{\partial \dot{q}}{\partial t}$ (16) fulfilled the set $\mathcal{X}_{\Delta \tau}$ as can be seen in Fig. 8.

## 6.2 Swift box-grabbing with the HRP-4 Robot

The HRP-4 robot swiftly grabbed a 1.08 kg box and another 0.38 kg box. The target contact velocity is 0.25 m/s. The robot did not reduce speeds for establishing contacts in both cases. The impacts are synchronized for the two hands, and the contact locations are symmetric.

**6.2.1 Data acquisition** The HRP-4 is also high-stiffness kinematic-controlled (assumption A.3). The HRP-4's geometric size is significantly higher than the contact area (assumption A.5). Thus, we adopt the same impact model used in Experiment 1. We bounded the coefficient of restitution by $\bar{c}_r = 0.3$ and $\bar{c}_r = 0.0$. We mounted coarse material with friction coefficient $\mu = 0.79$ on the two palms.

The HRP-4 has two ATI sensors attached to the wrists. We can timely and reliably detect the collisions by thresholding the force measurement at 20 N.

We measured the impulses in Fig. 5(c) and Fig. 6(c) by sampling the force torque sensors at 25000 Hz. On the other hand, the HRP-4's control system runs at 200 Hz. To
<figure><img src="image_7.png" /><figcaption>(a) Panda impact configuration one</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Normal contact velocity</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) Impulse</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) Contact force</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(e) COM velocity</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(f) Angular momentum</figcaption></figure>

**Figure 5.** For the Panda impact configuration in Fig. 5(a), the impact-aware QP solver adjusted the normal contact velocity, as depicted in Fig. 5(b). We visualize the polyhedra that constrain the impulse (12), the peak contact force (14), the COM velocity (17), and the angular momentum (18). All the measurements (red crosses) satisfy the corresponding polyhedron-represented sets.

<figure><img src="image_13.png" /><figcaption>(a) Panda impact configuration two</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) Normal contact velocity</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(c) Impulse</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(d) Contact force</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(e) COM velocity</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(f) Angular momentum</figcaption></figure>

**Figure 6.** For the Panda impact configuration in Fig. 6(a), the impact-aware QP solver adjusted the normal contact velocity, as depicted in Fig. 6(b). We visualize the polyhedra that constrain the impulse (12), the peak contact force (14), the COM velocity (17), and the angular momentum (18). All the measurements (red crosses) satisfy the corresponding polyhedron-represented sets.
<figure><img src="image_19.png" /><figcaption>**Figure 7.** The measured and predicted joint velocity jumps of the Panda manipulator during different trials of *Experiment 1*. The inequality constraint (29) was inactive as the post-impact joint velocities remained well within the bounds defined by (37).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 8.** The measured and predicted impulsive joint torque jumps Δγ during various trials of *Experiment 1*. The inequality constraint (30) was active since the QP solution of joint 6 reached the bound. Due to the conservative tuning of parameters, the measured impulsive joint torque jump Δγ° was observed to be lower than the corresponding QP solution.</figcaption></figure>

study the effects of the low sampling rate on the impulse calculation, we compared the same impact experiments against different sampling rates; Fig. 11 collects four sets of contact force profiles corresponding to 25000 Hz, 1000 Hz, 500 Hz, and 200 Hz. As the mean impulses are similar, we concluded that it is acceptable to integrate force measurement sampled at 200 Hz to compute the impulses.

We placed two contact point frames at the two sides of the boxes. The frame axes are defined following a similar configuration as shown in Fig. 2, in alignment with the impact mechanics definitions (Stronge 2000). We visualize their orientations in Fig. 9. At the moment of impact, the origin of the two contact point frames shared the same translation with the corresponding end-effector.

### 6.2.2 Controller formulation
We formulated the impact-aware constraints (36) for the two arms independently. The impact prediction employed the centroidal momentum matrix considering the kinematic chain from the floating-base link to the end-effector. The chain** weighs 11.63 kg.

Given the inertial frame $F_O$, the floating-base frame $F_B$, and the end-effector frame $F_e$, the end-effector's body

<figure><img src="image_21.png" /><figcaption>**Figure 9.** We color x, y and z axes with red, green and blue, respectively. The two green polyhedra align with the two contact frames, representing the impulse sets $\mathcal{X}_i$ for the two arms during the second grabbing. It should be noted that the scale of the polyhedra does not correspond to the robot's geometry.</figcaption></figure>

velocity writes††:

$$V_{Oe}^{b} = Ad_{gBe}^{-1}V_{OB}^{b} + V_{Be}^{b},$$

where $V_{OB}^b \in \mathbb{R}^6$ parameterizes the floating-base velocity. Thus, the body velocity Jacobian $J_{Oe}^b$ for the kinematic chain with a floating-base joint is:

$$J_{Oe}^{b} = [Ad_{gBe}^{-1}, J_{Be}^{b}]. \quad (38)$$

We align $J_{Oe}^b$ to the contact point frame with (19) as before.

Reactively optimizing the following QP generates the impact-aware whole-body motion (prior to impact):

$$ \begin{array}{ll} \min_{\nu, \lambda_{\mu}} & \sum_{i \in \mathcal{X}_o} w_i \|e_i(\nu)\|^2 \\ \text{s.t.} & \text{Impact-awareness for the left arm (35),} \\ & \text{Impact-awareness for the right arm (35),} \\ & \text{Joint position, velocity, and torque (2),} \end{array} \quad (39) $$

where the impact-aware constraints (35) do not restrict the floating-base joint velocity or torque.

The QP (39) synchronized the two impacts to grab the box without exerting unnecessary rotating moments. The two arms followed pre-defined (off-line planned) trajectories considering the approximate location of the box. Upon impact detection, the QP activates the admittance task for the robot to firmly grab the box and then toss it to the ramp located on its left side. Following another pre-defined trajectory, the HRP-4 robot repeated the grabbing and tossing motion for the second box before resuming the initial configuration. The impact-aware constraints were fulfilled for the impacts.

<footnote>**The total weight of the HRP-4 robot is 34.62 kg without the batteries.</footnote>

<footnote>††Body velocity transform (Murray et al. 1994, Proposition 2.15).</footnote>
<figure><img src="image_22.png" /><figcaption>**Figure 10.** Snapshot of *Experiment 2*: the HRP-4 robot grabbed two boxes in a row with contact velocities of 0.25 m/s. The green-highlighted *Grab box* states applied the impact-aware QP (34). During the blue-highlighted *Toss box* states, the robot grabbed the box with admittance control.</figcaption></figure>

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

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

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

<figure><img src="image_26.png" /><figcaption>**Figure 11.** The measured force curves in Fig. 11(a)-11(d) were collected by repeating the same experiment: impacting the ATI force-torque sensor at 0.15 m/s by the Panda robot with the joint configuration shown in Fig. 2. Despite the variations in sampling frequencies, the mean impulses $\bar{\iota}$ (obtained by integrating the force over time) are similar. Consequently, we can accurately measure the impulses of the HRP-4 robot during experiments by utilizing force measurements sampled at 200 Hz.</figcaption></figure>

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

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

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

<figure><img src="image_30.png" /><figcaption>**Figure 12.** The HRP-4's *Experiment 2* includes four impacts (two arms with two grabbings each). The measured impulses in all situations are constrained within the impulse set $\mathcal{I}$. Figure 12(a) and 12(b) depict the impulse sets of the left and right arms during the first grabbing, while Figure 12(c) and 12(d) correspond to the second grabbing.</figcaption></figure>

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

<figure><img src="image_32.png" /><figcaption>**Figure 13.** The robot accurately tracked the reference contact velocities in *Experiment 2*. The light-blue areas in Fig. 13(a), 13(b) indicate the period during which the controller activated the impact-aware constraints (35).</figcaption></figure>
6.2.3 **Result analysis** The HRP-4 robot established the two contacts at 0.25 m/s, without slowing down or following pre-defined deceleration trajectories (H-2.1). We synchronized the two impacts (H-2.2) to grab the boxes.

We independently model the two synchronized impacts on the two palms, see Fig. 9. Given the ATI force sensor data, Fig. 12 shows that the measured impulses are within the predicted impulse sets regardless of different situations.

Fig. 13(a) and Fig. 13(b) plot the inertial frame end-effector velocities for the first and second grabbing. Contrary to *Experiment 1*, in both grabbing cases, the robot precisely tracked the reference contact velocities (as they are feasible w.r.t the embedded other constraints).

## 7 Conclusion and Future Work

We aim to enhance the task-space QP controller to deal with intended impact tasks. To the authors best knowledge, our paper is the first to integrate frictional impacts in three dimensions into such optimization-based controllers.

We construct the impulse polyhedron to cover all the candidate solutions that fulfill Coulomb's friction law and task-space momentum conservation. By projecting the polyhedra (half-space representation) into joint or task space, the controller gains awareness of post-impact states. As we are interested in a conservative solution, the polyhedra can accommodate model uncertainties by decreasing the angle of the friction cone or the restitution coefficients' bounds.

The impact-aware constraints, represented as convex polyhedra, modify the whole-body QP's search space according to the impact model and hardware-affordable resilience bounds. We assessed our approach with two robots: the HRP-4 humanoid and the Panda manipulator, both of which achieved high contact velocities without exceeding the hardware's resilience limits.

In future work, we plan to address certain aspects from a broader task-space QP control perspective, including robustness to handling interacting, potentially conflicting tasks, and the robust activation of constraints on the fly, regardless of their nature. We are currently obtaining promising results by formulating the QP as an MPC governor.

For aspects dealing more with impacts. We still have no solution if the critical bounds are not the correct ones (except being restrictive or conservative). Additionally, there is a need to devise a good shock propagation model for both fixed and floating based robots. Our observations have also highlighted the necessity to reevaluate the standard dynamic balance criteria for humanoid robots under external impacts. We are currently preparing a separate paper dedicated to this topic. The humanoid experiments also revealed that handling multiple non-synchronous impacts on a moving object (floating box) is another shortcoming. This problem can be easily understood from a simple 2D toy example: rapidly grasping a 2D floating box (i.e., a box that can move under external forces) with two points robots controlled under task-space QP. In such cases, where the robots do not impact the object simultaneously, the object's mobility on the side opposite to the first impact can generate a higher relative velocity than expected. Although, by the time we wrote

the paper, an interesting result was disclosed by Aydinoglu et al. (2022), it is important to note that their work did not integrate resilience constraints. We will thoroughly examine the scenario of multiple impacts on moving objects.

## Acknowledgment

This work is supported by the EU H2020 research grant GA 871899, I.AM. project. We thank Pierre Gergondet for his continuous support in setting up the mc_rtc controller.

## References

Abe Y, Da Silva M and Popović J (2007) Multi-objective control with frictional contacts. In: *ACM SIGGRAPH/Eurographics symposium on Computer animation*. pp. 249–258.

Aizerman MA and Gantmacher FR (1958) Determination of stability by linear approximation of a periodic solution of a system of differential equations with discontinuous right-hand sides. *The Quarterly Journal of Mechanics and Applied Mathematics* 11(4): 385–398.

Aydinoglu A, Sieg P, Preciado VM and Posa M (2022) Stabilization of complementarity systems via contact-aware controllers. *IEEE Transactions on Robotics* 38(3): 1735–1754.

Battaglia M, Blanchet L, Kheddar A, Kajita S and Yokoi K (2009) Combining haptic sensing with safe interaction. In: *IEEE/RSJ Int. Conf. on Intelligent Robots and Systems*. pp. 231–236.

Bombile M and Billard A (2022) Dual-arm control for coordinated fast grabbing and tossing of an object: Proposing a new approach. *IEEE Robotics and Automation Magazine* 29(3): 127–138.

Bouyarmane K, Chappellet K, Vaillant J and Kheddar A (2019) Quadratic programming for multirobot and task-space force control. *IEEE Transactions on Robotics* 35(1): 64–77.

Bouyarmane K and Kheddar A (2018) On weight-prioritized multitask control of humanoid robots. *IEEE Transactions on Automatic Control* 63(6): 1632–1647.

Brogliato B (2003) Some perspectives on the analysis and control of complementarity systems. *IEEE Transactions on Automatic Control* 48(6): 918–935.

Burden SA, Libby T and Coogan SD (2018) On contraction analysis for hybrid systems. *arXiv preprint arXiv:1811.03956*.

Chatterjee A and Ruina A (1998) A new algebraic rigid-body collision law based on impulse space considerations. *Journal of Applied Mechanics* 65(4): 939–951.

Dehio N, Wang Y and Kheddar A (2022) Dual-arm box grabbing with impact-aware model predictive control utilizing soft deformable end-effector pads. *IEEE Robotics and Automation Letters* 7(2): 5647–5654.

Del Prete A (2018) Joint position and velocity bounds in discrete-time acceleration/torque control of robot manipulators. *IEEE Robotics and Automation Letters* 3(1): 281–288.

De Magistris G, Pajon A, Miossec S and Kheddar A (2017) Optimized humanoid walking with soft soles. *Robotics and Autonomous Systems* 95: 52–63.

Djeha M, Gergondet P and Kheddar A (2023) Robust task-space quadratic programming for kinematic-controlled robots. *IEEE Transactions on Robotics*.
Djeha M, Tanguy A and Kheddar A (2020) Adaptive-gains enforcing constraints in closed-loop QP control. *IEEE Robotics and Automation Letters* 5(4): 6504–6511.

Gong Y and Grizzle J (2020) Angular momentum about the contact point for control of bipedal locomotion: Validation in a lip-based controller. *arXiv preprint arXiv:2008.10763*.

Grizzle JW, Chevallereau C, Sinnet RW and Ames AD (2014) Models, feedback control, and open problems of 3d bipedal robotic walking. *Automatica* 50(8): 1955–1988.

Halm M and Posa M (2019) Modeling and analysis of non-unique behaviors in multiple frictional impacts. In: *Proceedings of Robotics: Science and Systems*, volume 15. FreiburgimBreisgau, Germany.

Halm M and Posa M (2021) Set-valued rigid body dynamics for simultaneous frictional impact. *Preprint arXiv:2103.15714*.

Heck D, Saccon A, Van de Wouw N and Nijmeijer H (2016) Guaranteeing stable tracking of hybrid position-force trajectories for a robot manipulator interacting with a stiff environment. *Automatica* 63: 235–247.

Hu G, Makkar C and Dixon WE (2007) Energy-based nonlinear control of underactuated euler-lagrange systems subject to impacts. *IEEE Transactions on Automatic Control* 52(9): 1742–1748.

Hurmuzlu Y, Génot F and Brogliato B (2004) Modeling, stability and control of biped robots-a general framework. *Automatica* 40(10): 1647 – 1664.

Jia YB, Gardner M and Mu X (2019) Batting an in-flight object to the target. *The Int. Journal of Robotics Research* 38(4): 451–485.

Jia YB and Wang F (2017) Analysis and computation of two body impact in three dimensions. *Journal of Computational and Nonlinear Dynamics* 12(4): 041012 (16 pages).

Johnson AM, Burden SA and Koditschek DE (2016) A hybrid systems model for simple manipulation and self-manipulation systems. *The Int. Journal of Robotics Research* 35(11): 1354–1392.

Kajita S, Morisawa M, Miura K, Nakaoka S, Harada K, Kaneko K, Kanehiro F and Yokoi K (2010) Biped walking stabilization based on linear inverted pendulum tracking. In: *IEEE/RSJ Int. Conf. on Intelligent Robots and Systems*. pp. 4489–4496.

Keller JB (1986) Impact with friction. *Journal of Applied Mechanics* 53(1): 1–4.

Kheddar A, Caron S, Gergondet P, Comport A, Tanguy A, Ott C, Henze B, Mesesan G, Englberger J, Roa MA et al. (2019) Humanoid robots in aircraft manufacturing: The airbus use cases. *IEEE Robotics & Automation Magazine* 26(4): 30–45.

Khulief Y (2013) Modeling of impact in multibody systems: an overview. *Journal of Computational and Nonlinear Dynamics* 8(2): 021012(15 pages).

Khurana H, Bombile M and Billard A (2021) Learning to hit: A statistical dynamical system based approach. In: *IEEE/RSJ International Conference on Intelligent Robots and Systems*. pp. 9415–9421.

Konno A, Myojin T, Matsumoto T, Tsujita T and Uchiyama M (2011) An impact dynamics model and sequential optimization to generate impact motions for a humanoid robot. *The Int. Journal of Robotics Research* 30(13): 1596–1608.

Kuindersma S, Deits R, Fallon M, Valenzuela A, Dai H, Permenter F, Koolen T, Marion P and Tedrake R (2016) Optimization-based locomotion planning, estimation, and control design for the atlas humanoid robot. *Autonomous Robots* 40(3): 429–455.

Lankarani HM (2000) A poisson-based formulation for frictional impact analysis of multibody mechanical systems with open or closed kinematic chains. *Journal of Mechanical Design* 122(4): 489–497.

Lee SH and Goswami A (2012) A momentum-based balance controller for humanoid robots on non-level and non-stationary ground. *Autonomous Robots* 33(4): 399–414.

Liu M, Lober R and Padois V (2016) Whole-body hierarchical motion and force control for humanoid robots. *Autonomous Robots* 40(3): 493–504.

Manchester Z, Doshi N, Wood RJ and Kuindersma S (2019) Contact-implicit trajectory optimization using variational integrators. *The Int. Journal of Robotics Research* 38(12-13): 1463–1476.

Murray RM, Li Z, Sastry SS and Sastry SS (1994) *A mathematical introduction to robotic manipulation*. CRC press.

Orin DE and Goswami A (2008) Centroidal momentum matrix of a humanoid robot: Structure and properties. In: *IEEE/RSJ Int. Conf. on Intelligent Robots and Systems*. pp. 653–659.

Orin DE, Goswami A and Lee SH (2013) Centroidal dynamics of a humanoid robot. *Autonomous robots* 35(2-3): 161–176.

Pace AM and Burden SA (2017) Piecewise-differentiable trajectory outcomes in mechanical systems subject to unilateral constraints. In: *20th International Conference on Hybrid Systems: Computation and Control*. ASME, pp. 243–252.

Pagilla PR and Yu B (2001) A stable transition controller for constrained robots. *IEEE/ASME transactions on mechatronics* 6(1): 65–74.

Pajon A, Caron S, De Magistris G, Miossec S and Kheddar A (2017) Walking on gravel with soft soles using linear inverted pendulum tracking and reaction force distribution. In: *IEEE-RAS Int. Conf. on Humanoid Robots*. p. 432–437.

Pashah S, Massenzio M and Jacquelin E (2008) Prediction of structural response for low velocity impact. *Int. Journal of Impact Engineering* 35(2): 119–132.

Posa M, Cantu C and Tedrake R (2014) A direct method for trajectory optimization of rigid bodies through contact. *The Int. Journal of Robotics Research* 33(1): 69–81.

Remy CD (2017) Ambiguous collision outcomes and sliding with infinite friction in models of legged systems. *The Int. Journal of Robotics Research* 36(12): 1252–1267.

Righetti L, Buchli J, Mistry M, Kalakrishnan M and Schaal S (2013) Optimal distribution of contact forces with inverse-dynamics control. *Int. Journal of Robotics Research* 32(3): 280–298.

Rijnen M, de Mooij E, Traversaro S, Nori F, van de Wouw N, Saccon A and Nijmeijer H (2017) Control of humanoid robot motions with impacts: Numerical experiments with reference spreading control. In: *IEEE Int. Conf. on Robotics and Automation*. pp. 4102–4107.

Routh EJ (1955) *Dynamics of a system of rigid bodies*. Dover Publications, New York.

Shield S, Johnson AM and Patel A (2022) Contact-implicit direct collocation with a discontinuous velocity state. *IEEE Robotics and Automation Letters* 7(2): 5779–5786.
Siciliano B and Khatib O (2016) Springer handbook of robotics. Springer.

Singh BRP and Featherstone R (2020) Mechanical shock propagation reduction in robot legs. *IEEE Robotics and Automation Letters* 5(2): 1183–1190.

Stanisic RZ and Fernández ÁV (2012) Adjusting the parameters of the mechanical impedance for velocity, impact and force control. *Robotica* 30(4): 583–597.

Stewart DE (2000) Rigid-body dynamics with friction and impact. *SIAM review* 42(1): 3–39.

Stoianovici D and Hurmuzlu Y (1996) A critical study of the applicability of rigid-body collision theory. *Journal of Applied Mechanics* 63(2): 307.

Stronge WJ (2000) *Impact mechanics*. Cambridge university press.

Sugihara T (2009) Standing stabilizability and stepping maneuver in planar bipedalism based on the best com-zmp regulator. In: *Int. Conf. on Robotics and Automation*. pp. 1966–1971.

van Steen J, van den Brandt G, van de Wouw N, Kober J and Saccon A (2023) Quadratic programming-based reference spreading control for dual-arm robotic manipulation with planned simultaneous impacts. *arXiv. 2305.08643*.

van Steen JJ, van de Wouw N and Saccon A (2022) Robot control for simultaneous impact tasks via quadratic programming-based reference spreading. In: *American Control Conference*. IEEE, pp. 3865–3872.

Wang Y, Dehio N and Kheddar A (2022a) On inverse inertia matrix and contact-force model for robotic manipulators at normal impacts. *IEEE Robotics and Automation Letters* 7(2): 3648–3655.

Wang Y, Dehio N and Kheddar A (2022b) Predicting post-impact joint velocity jumps on kinematics controlled manipulators. *IEEE Robotics and Automation Letters* 7(3): 6226–6233.

Wang Y and Kheddar A (2019) Impact-friendly robust control design with task-space quadratic optimization. In: *Robotics: Science and Systems*, volume 15. Freiburg, Germany.

Wang Y and Mason MT (1992) Two-dimensional rigid-body collisions with friction. *Journal of Applied Mechanics* 59(3): 635.

Wang Y, Tanguy A, Gergondet P and Kheddar A (2019) Impact-aware multi-contact motion generation with a quadratic optimization controller. In: *IEEE-RAS Int. Conf. on Humanoid Robots*.

Wang Y, Vina F, Karayiannidis Y, Smith C and Ögren P (2014) Dual arm manipulation using constraint based programming. In: *19th IFAC World Congress*, volume 47. Cape Town, South Africa: Elsevier, pp. 311–319.

Wiedebach G, Bertrand S, Wu T, Fiorio L, McCrory S, Griffin R, Nori F and Pratt J (2016) Walking on partial footholds including line contacts with the humanoid robot atlas. In: *IEEE-RAS Int. Conf. on Humanoid Robots*. pp. 1312–1319.

Yang W and Posa M (2021) Impact invariant control with applications to bipedal locomotion. In: *IEEE/RSJ Int. Conf. on Intelligent Robots and Systems*. pp. 5151–5158.

Yang W and Posa M (2023) Impact-invariant control: Maximizing control authority during impacts. *arXiv. 2303.00817*.

Zhang Y, Ge SS and Lee TH (2004) A unified quadratic-programming-based dynamical system approach to joint torque optimization of physically constrained redundant manipulators. *IEEE Transactions on Systems, Man, and Cybernetics, Part B*

34(5): 2126–2132.

Zheng YF and Hemami H (1985) Mathematical modeling of a robot collision with its environment. *Journal of Field Robotics* 2(3): 289–307.