# Optimal 3D arm strategies for maximizing twist rotation during somersault of a rigid-body model

François Bailly, Eve Charbonneau, Loane Danès, Mickael Begon

► To cite this version:

François Bailly, Eve Charbonneau, Loane Danès, Mickael Begon. Optimal 3D arm strategies for maximizing twist rotation during somersault of a rigid-body model. Multibody System Dynamics, 2020, 52 (2), pp.193-209. 10.1007/s11044-020-09759-5. hal-05379211

HAL Id: hal-05379211

https://inria.hal.science/hal-05379211v1

Submitted on 24 Nov 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimal 3D Arm Strategies For Maximizing Twist Rotation During Somersault of a Rigid-Body Model

François Bailly<sup>a,*</sup>, Eve Charbonneau<sup>a</sup>, Loane Danès<sup>b</sup> and Mickael Begon<sup>a</sup>

**Abstract**— Looking for new arm strategies for better twisting performances during a backward somersault is of interest for the acrobatic sports community while being a complex mechanical problem, due to the non-linearity of the dynamics involved. As the pursued solutions are not intuitive, computer simulation is a relevant tool to explore a wider variety of techniques. Simulations of twisting somersaults have mainly been realized with planar arm motions. The aim of this study was to explore the outcomes of using 3D techniques, with the demonstration that increasing the fidelity of the model does not increase the level of control complexity on the real system. Optimal control was used to maximize twists in a backward straight somersault with both types of models. A multi-start approach was used to find large sets of near-optimal solutions. The robustness of these solutions was then assessed by modeling kinematic noise during motion execution. The possibility of using quaternions for representing orientations in this numerical optimization problem was discussed. Optimized solutions showed that 3D techniques generated about two additional twists compared to 2D techniques. The robustness analysis revealed clusters of highly-twisting and stable 3D solutions. This study demonstrates the superiority of 3D solutions for twisting in backward somersault, a result that could help acrobatic sports athletes to improve their twisting performance.

**Keywords** – Rigid Body Dynamics, Biomechanics; Optimal Control; Numerical Optimization, Sports Performance

## I. INTRODUCTION

Twist rotations during aerial somersaults in acrobatic sports have been frequently studied in biomechanics [1], [2], [3]. Apart from contact techniques which involve a twisting angular momentum at takeoff, there exist two main categories of techniques which have sparked researchers interest, namely aerial and cat-twist techniques [4]. Understanding the motor strategies underlying these approaches is an interesting challenge because they involve complex sequences of actions subject to highly non-linear dynamics. Pure aerial techniques–as opposed to cat twists–are of particular interest in acrobatic sports since they can be performed using only arm motions, leaving the rest of the body in a straight configuration, which is often preferred for aesthetic reasons, and simpler to model in simulation. In most of the acrobatic/artistic activities (trampoline, diving, gymnastics, freestyle skiing, etc.), the number of twists during somersaults is considered as a part of the sports performance. Looking for innovative arm strategies which enhance the twisting performance during a somersault

<figure><img src="image_2.png" /><figcaption>Fig. 1: Model definition with the 10 DOFs being: antero-posterior ($q_{tx}^{ra}$), transverse ($q_{ty}^{ra}$), and longitudinal ($q_{tz}^{ra}$) translations of the root, antero-posterior axis ($q_{rx}^{r}$, tilt), transverse axis ($q_{ry}^{r}$, somersault), and longitudinal axis ($q_{rz}^{r}$, twist) rotations of the root, right arm ($q_{el}^{ra}$) and left arm ($q_{el}^{la}$) abduction/adduction (elevation), right arm ($q_{cp}^{ra}$) and left arm ($q_{cp}^{la}$) abduction/adduction change of plane.</figcaption></figure>

is therefore both an enticing dynamics problem and a relevant matter for the acrobatic sports community.

Computer simulation is a relevant tool for innovating in sports and especially in disciplines in which new techniques can increase the risk of acute injuries, which is precisely the case in acrobatic sports. In this context, numerical trajectory optimization methods are frequently used as they are efficient means of finding innovative techniques, in compliance with a set of sport-related specifications. There are essentially four main parameters to be considered in order to discuss the results of scientific studies on sports performance that made use of trajectory optimization: the numerical model of the body, the parametrization of the problem, the cost function paired with a set of constraints and the nature of the solver.

Concerning the model first, the many studies in which computer simulation was used to investigate twist rotation in acrobatic sports using arms techniques [5], [3], [6] only modeled arm lowering/rising in the frontal plane (i.e., abduction/adduction), ignoring the change in plane of elevation. Because of several non-linear effects due to

<footnote><sup>a</sup> Laboratoire de Simulation et Modélisation du Mouvement, Faculté de Médecine, Université de Montréal, Laval, QC, Canada</footnote>

<footnote><sup>b</sup> AgroParisTech, Paris, France</footnote>

<footnote>* francois.bailly@umontreal.ca</footnote>
rotations, putting aside the possibility of changing the plane of abduction/adduction does not necessarily simplify the movement from a motor control point of view. Indeed, the fact that the body is moving freely in space implies that the gravity is not always pushing on the arm in the same orientation, and that Coriolis effect pushes the arm in a direction perpendicular to its trajectory. Therefore, abducting the arm during a twisting somersault requires additional joint torques compared to abducting the arm in a static position [7]. In order to numerically quantify this additional amount of torque and to demonstrate the improvements brought by 3D motions of the arms for generating twists during somersault (see Figs. 4, 5), two models were selected and compared in this study. In accordance with the literature, the first one included a free floating base and one degree of freedom (DoF) at each shoulder (6 + 2 = 8 DoFs total), whereas the second one includes two DoFs at each shoulder (6 + 4 = 10 DoFs total).

Once a model has been selected, comes the question of the optimization problem parametrization. In particular, what are the variables chosen for the optimization and how to discretize all the continuous functions that appear in the problem statement? The goal of this phase, called transcription, is to transform a specific trajectory optimization problem into a generic non-linear program (NLP) that will be solved using a dedicated algorithm. A first family of methods comes from numerical optimal control and are known as *direct* methods. They consist in straightforwardly choosing the state and/or the control as optimization variables at a given number of points along the trajectory and they rely on explicit or implicit integration for enforcing the dynamics between these points. For instance, *direct multiple shooting* was applied to the aerial phase of twisting somersaults in diving giving satisfying results [2], however the goal was not to maximize performance, but to minimize the joint torques for given sports performances in terms of somersaults and twists. *Direct collocation* methods have also shown their efficiency in several studies investigating human motion [8], [9], but their convergence might be challenging, as both the state and the control must be provided as initial guesses. Other choices can be made, as in [10], where the optimization variables are instants at which a switch in the motor strategy occurs, using polynomials function (4th, 5th order) in-between. Thus, this kind of parametrization requires an a priori knowledge about entire pieces of the state trajectory which is not suitable in the quest of innovative strategies. In [11], the optimization variables are the coefficients of fourth order polynomial approximations of the states, with linking conditions to enforce the continuity of the controls. Avoiding large variations of joints torques at discretization points is not a major concern of the present study, in which the search for innovative motor strategies is favored. Thus, a classical direct multiple shooting optimal control problem (OCP) was implemented, with optimization variables being the joints torques.

Then, concerning the cost functions in sports performance optimization studies, various objectives were already proposed: minimal time of execution [12], minimal joint torques [2], maximal execution score [13], etc. Since the

aim of the present work was to compare two families of techniques (2D versus 3D) for generating twist rotations, we were interested in the best performance possible regarding this particular task. Thus, the main objective function was to maximize the number of twists during a straight backward somersault. As detailed further (see section II-B), additional penalties were added to this main objective, in order to regularize the problem and to yield more realistic solutions. The proposed OCP was constrained by acrobatic-sports-related initial and terminal requirements such as the take-off and landing positions and velocities of the segments. Kinematic and kinetic path constraints were also added to reflect realistic positions, velocities and torques of the numerical avatar. When parameterized with Euler angles (termed as somersault, tilt and twist [14]), rotations of the root segment during somersault combined with twists are likely to yield gimbal lock resulting in numerical instabilities. Thus, in a first approach, the second Euler angle of the sequence (i.e. tilt) was bounded in order to avoid such a situation. In order to investigate whether this restriction might prevent the solver from finding more optimal solutions, the root segment was also parametrized using quaternions and both approaches are discussed in section II-C.

Regarding the nature of the solver, lots of options exist in the literature for solving NLPs and some of them are available as third party software. In several studies, the simulated annealing algorithm was employed to optimize the movement and find a global optimum [10], [15]. This method has shown its efficiency, but is restrictive as only a few parameters can be optimized ([16]). For instance, it is well suited for optimizing a reduced set of high-level parameters such as the key moments at which the motor strategy should be changed [10]. Sequential quadratic programming and quasi-Newton algorithms have been tested on biomechanical problems with convergence issues due to their high sensitivity to variable scaling [17]. Interior point methods led to satisfying results ([18], [19]), their main limiting factor being an important computational time. The aim of the present study being to explore new techniques, the initial conditions and scaling factors are unknown (especially for the control), thus the resolution method should not be too sensitive to them. To explore a wide variety of techniques, the problem's complexity has to be preserved, therefore making strong assumptions about the mathematical properties of the solutions or reducing the number of parameters to be optimized are not reasonable compromises to make. Because the goal of this study was not to optimize in real time, reasonably longer computational time could be afforded. For these reasons, the interior point method was more appropriate than simulated annealing, sequential quadratic programming or quasi-Newton algorithm. Having in mind that gradient based optimizers do not guarantee convergence towards a global extremum, when optimizing human movements, a series of local extrema might represent satisfying compromises leading to more robust techniques without compromising a lot the performance. This is particularly the case in this study, where the robustness analysis is not part of the
<figure><img src="image_3.png" /><figcaption>Fig. 2: Snapshots of the most-twisting 2D technique, without penalization on the hands trajectories. The final number of twists achieved by this technique is 2.94. For a visualization of the corresponding arms trajectories, see Fig. 8 in Appendix.</figcaption></figure>

optimized objective, but rather studied a posteriori, as detailed in Sec. II-E. Therefore, in this context, considering near-optimal solutions gives an opportunity to explore a wider range of techniques. To maximize the amount of solutions found, multiple optimizations can be run in parallel with slight variations in the initial guesses and in the settings of the problem. Here, a multi-start strategy was used to generate several techniques as in [20]. Aside from technical considerations, another aspect of our problematic comes from the fact that even elite athletes can experience little variations in the execution of their skills from one trial to another [21]. To make sure that the proposed optimal techniques can actually be executed by athletes, it is needed to evaluate their robustness, as techniques needing infinite precision are of no interest to the sports community. To simulate athletes technique variations, noise can be added to the positions of the model's segments [22]. Thus, after optimization, Gaussian noise was added to each joint in order to modify the optimal techniques and evaluate the impact of this disturbance on the performance criterion.

The purpose of this study was to optimize 3D arm techniques for generating twist rotations during a backward somersault and compare the results with state-of-the-art 2D techniques. For this purpose, different optimization problems were designed and solved using a multiple-shooting formulation and a NLP solver. The paper is organized as follows: first, the numerical models of the body are described (Sec. II-A), followed by a presentation of the OCPs settings, with an emphasis on the choice of the cost functions and regularization terms (Sec. II-B). After that, we provide an original discussion on the choice between Euler angles and quaternions for the parametrization of the root segment (Sec. II-C). Then, the multi-start approach and a method for computing the robustness of the solutions post-optimization are described (Sec. II-D, II-E). Afterwards, we present an analysis of the torques induced by non-linear effects during

somersault, which justifies the choice of our model (Sec. II-F). The contribution of this work is twofold: first it demonstrates that 3D arm techniques outperform 2D ones in terms of twist performance, secondly it shows that for a same level of performance, 3D arm techniques are simpler than 2D ones, which should be of interests to acrobatic sports athletes and coaches.

## II. METHODS

### A. Skeletal Model

The dynamic models of the body were developed in Biorbd [23], a RBDL-based [24] C++ library according to the spatial vector algebra described in [25] to compute the equations of motion. The first model was composed of three rigid bodies (8 DoFs total) connected with inelastic joints: the root segment that includes the head, trunk and legs, the right upper-limb and the left upper-limb, both with extended elbows (Fig. 1). The root segment was modeled with six DoFs to allow 3D translations and 3D rotations of the whole body in space, and each shoulder was modeled with one DoF for abducting and adducting the arms. It is further referred to as the *2D model*, because of the planar motion of the arms. The second model was based on the first one, with the addition of one extra DoF per shoulder (for changing the plane of elevation), bringing the total number of DoFs to 10. It is further referred to as the *3D model*. Masses and inertia matrices of the models were computed from 95 anthropometric measurements of one participant [14].

### B. Optimal control formulation

The 2D (resp. 3D) model was driven by 2 (resp. 4) torque actuators corresponding to the DOFs of the shoulders. In this study, the models being in free fall conditions, only gravity and non-linear forces coming from rotations are considered. The model dynamics is given by the following equation, including the motion of the free-floating base:
<figure><img src="image_4.png" /><figcaption>Fig. 3: Snapshots of the most-twisting 3D technique, without penalization on the hands trajectories. The final number of twists achieved by this technique is 4.90. For a visualization of the corresponding arms trajectories, see Fig. 9 in Appendix</figcaption></figure>

$$ \begin{bmatrix} 0_{6 \times 1} \\ \tau \end{bmatrix} = M(q)\ddot{q} + N(q, \dot{q}) + G(q), \quad (1) $$

where $\tau$ is the joint-torque vector, M is the mass matrix. $q, \dot{q}$ and $\ddot{q}$ are the base and joints positions, velocities and accelerations vectors, all in generalized coordinates. $N(q, \dot{q})$ represents the Coriolis and centrifugal effects and $G(q)$ represents the effects of gravity. Let $\dot{x} \triangleq [q(t), \dot{q}(t)]$ be the state of the system. Let $\mathbf{u} \triangleq [\tau]$ be the control applied to the system. Let us denote the time evolution of a variable, i.e. its trajectory, by underlining it ($\underline{x} = [\underline{x}(0), \dots, \underline{x}(T)]$). The forward dynamics which solves $\dot{x} = f(\underline{x}, \mathbf{u})$, can be derived from Eq. (1), as it amounts to compute $\ddot{q}$ as:

$$ \ddot{q} = M(q)^{-1}[\tau - N(q, \dot{q}) - G(q)]. \quad (2) $$

In order to find an optimal control law which maximizes the number of twists while satisfying the dynamics and the constraints of the problem, an optimal control problem (OCP) of the following form was set up:

$$ \min_{\underline{x}, \underline{u}, T} \quad C(\underline{x}, \underline{u}) = \int_{0}^{T} \ell(\underline{x}, \underline{u}) dt + m(\underline{x}(T)) \quad (3a) $$

$$ \text{s.t.} \quad \forall t, \quad \dot{x} = f(x, u) \quad (3b) $$

$$ \forall t, \quad \underline{u} \in \mathcal{U} \quad (3c) $$

$$ \forall t, \quad \underline{x} \in \mathcal{X} \quad (3d) $$

$$ x(0) = x_0 \quad (3e) $$

$$ x(T) = x_T \quad (3f) $$

where Eq. (3a) is the cost function to be minimized, with a Lagrange running term $\ell(\underline{x}, \underline{u})$ and a Mayer terminal term $m(\underline{x}(T))$ (with $T$, the duration of the motion, also an optimization variable). The first component of the Mayer term was the pelvis twist rotation at the end of the motion ($q_{rz}^r(T)$).

Hence, minimizing it amounts to maximize the number of negative (i.e. right) twists of the model. The second component of the Mayer term was the initial somersault velocity which was expected to be enough to yield an entire somersault with twist rotations but not too large, in order to obtain realistic solutions. The Lagrange term $\ell(\underline{x}, \underline{u})$ was either a regularization terms on the control variable $\mathbf{u}$ ($\ell_{UHP}$) or a combination of the latter and of a penalization term on the distance traveled by each hand ($\ell_{PHP}$). The solutions resulting from the first penalization are further called "Unpenalized Hand Path" (UHP) techniques, whereas the ones resulting from the second type of penalization are further called "Penalized Hand Path" (PHP) techniques. Weighting scalars were used to modulate the relative penalty of the Mayer and Lagrange terms (see Tab. (II)). Specifically, let the time $t$ range from 0 to $T$ (takeoff to landing), $\alpha, \beta_{qi}, \gamma$ and $\delta$ be the weighting coefficients. Let $N$ be the number of shooting nodes and $\tau_{ij}$ the control of the $i^{th}$ DoF at the $j^{th}$ node. Let $P_k$ be the $k^{th}$ hand path in space. The following equation presents the discretized cost function $C$ that is to be minimized:

$$ C_{UHP} = \alpha q_{rz}^r(T) + \gamma \dot{q}_{ry}^r(0) + \sum_{i=1}^{N_{\tau}} \beta_{qi} \sum_{j=0}^{N} \tau_{ij}^2, \quad (4) $$

$$ C_{PHP} = \alpha q_{rz}^r(T) + \gamma \dot{q}_{ry}^r(0) + \sum_{i=1}^{N_{\tau}} \beta_{qi} \sum_{j=0}^{N} \tau_{ij}^2 + \delta \sum_{k=1}^{2} \int_{P_k} ds, \quad (5) $$

where $N_{\tau} = 2$ (resp. 4) for the 2D (resp. 3D) model. $\beta_{qi}$ are coefficients for tuning the penalization of the control depending on the DoF considered. Finally, the $\int_{\gamma}$ symbol in the last term denotes a line integral along a path, which, in our case, represents the 3D distance traveled by each hand. In short, the Lagrange terms are regularizing the optimization at the kinematic (hand path) and kinetic (joint torques) levels for producing biomechanically meaningful solutions.
Eq. (3b) is the forward dynamics equation of the system. Eq. (3c) enforces that the control stays within a specified domain $\mathcal{U}$, typically $\mathbf{u}_{min} < \underline{\mathbf{u}} < \mathbf{u}_{max}$. Eq. (3d) enforces that the joints positions stay within their kinematic domain $\mathcal{X}$, typically for $i \in [0..N]$, $\dot{\mathbf{q}}_{min}^i < \underline{\mathbf{q}}^i < \dot{\mathbf{q}}_{max}^i$, with N the number of degrees of freedom of the model. Eqs. (3e) and (3f) constrain the initial and terminal states of the system, respectively. During the whole simulation the positions, velocities and torques of the arms were bounded to anatomically feasible movements (see Tab. I). The constraints applied to the model at the end point of the simulation were to finish the backward straight somersault at same level as the takeoff, with similar tilt ($q_5$) and arms above the head:

* Somersault rotation of the pelvis ($q_{ry}^r(T)$) in $360 \pm 10^\circ$

* Tilt of the pelvis ($q_{rx}^r(T)$) in $0 \pm 15^\circ$

* Arms abduction ($q_{el}^a(T)$ and $q_{el}^{la}(T)$) in $160 \pm 5^\circ$

Other constraints at the beginning and during the motion are listed in Tab. IV in appendix. All these constraints are derived from trampoline, enforcing the technique to be purely aerial ($\dot{q}_{rz}^r(0) = 0$) and ensuring successful landings consisting in upright position with arms above the head, for athletes to be prepared for the next movement of the routine. This optimal control problem was discretized using a direct multiple-shooting formulation, and then solved using CasADi [26] connected to the non-linear programming solver IPOPT [27].

<table><caption>TABLE I: Kinematic and torque constraints applied to the 3D and 2D models for solving the OCP. Arm position bounds reflect feasible range of motion [28], arm velocities and torques bounds are arbitrary high. However as torques are minimized, it is hypothesised that velocities and torques should not get to nonphysiological values.</caption><thead><tr><th rowspan="3"></th><th colspan="3">2D model</th><th colspan="3">3D model</th></tr><tr><th>q<sup>r</sup><sub>rz</sub></th><th>q<sup>ra</sup><sub>el</sub> &amp; q<sup>la</sup><sub>el</sub></th><th>q<sup>ra</sup><sub>cp</sub> &amp; q<sup>la</sup><sub>cp</sub></th><th>q<sup>r</sup><sub>rz</sub></th><th>q<sup>ra</sup><sub>el</sub> &amp; q<sup>la</sup><sub>el</sub></th><th>q<sup>ra</sup><sub>cp</sub> &amp; q<sup>la</sup><sub>cp</sub></th></tr><tr><th>position (°)</th><th>[-45, 45]</th><th>[0, 180]</th><th>[-135, 45]</th><th>position (°)</th><th>[-45, 45]</th><th>[0, 180]</th></tr></thead><tbody><tr><th>velocity (°/s)</th><td>—</td><td>[-5730, 5730]</td><td>[-5730, 5730]</td><th>velocity (°/s)</th><td>—</td><td>[-5730, 5730]</td></tr><tr><th>torque (Nm)</th><td>—</td><td>[-100, 100]</td><td>[-100, 100]</td><th>torque (Nm)</th><td>—</td><td>[-100, 100]</td></tr></tbody></table>

<table><caption>TABLE II: Weighting coefficients of the Lagrange and Mayer terms used in the cost function of the optimal control problem for PHP techniques. For the 8-dofs model, $\beta_{q_{cp}^{ra}}$ and $\beta_{q_{el}^{la}}$ are left blank because these DoFs are not included. For UHP techniques the coefficients are the same except for the $\delta$ column which becomes null.</caption><thead><tr><th></th><th>α</th><th>β<sub>q<sub>el</sub><sup>ra</sup></sub></th><th>β<sub>q<sub>el</sub><sup>la</sup></sub></th><th>β<sub>q<sub>cp</sub><sup>ra</sup></sub></th><th>β<sub>q<sub>el</sub><sup>la</sup></sub></th><th>γ</th><th>δ</th></tr></thead><tbody><tr><th>8-dofs model</th><td>-1e<sup>4</sup></td><td>1e<sup>-4</sup></td><td>1e<sup>-4</sup></td><td>-</td><td>-</td><td>1</td><td>1e<sup>-2</sup></td></tr><tr><th>10-dof model</th><td>-1e<sup>4</sup></td><td>1e<sup>-4</sup></td><td>1e<sup>-4</sup></td><td>1e<sup>-2</sup></td><td>1e<sup>-2</sup></td><td>1</td><td>1e<sup>-2</sup></td></tr></tbody></table>

### C. Quaternion modeling of the pelvis orientation

Unlike Euler angles, the description of 3D orientations by means of quaternions ($q$) is free of singularities, which is why

they are ideal mathematical tools for the computer graphics (visualization) or the robotics (control) communities. In order to avoid such singularities in the Euler angles representation presented above, the second angle of the pelvis rotation sequence was bounded to a domain excluding the possibility of gimbal lock ($q_{ry}^r(T)$ in $360 \pm 10^\circ$). However, this constraint is only met at the shooting nodes and might still lead to singularity during integration between two nodes. In addition to jeopardizing convergence this extra constraint might lead to less relevant solutions, although it is only a numerical artifice which has nothing to do with the kinematics or the dynamics of the problem.

In order to investigate to what extent this constraint could affect the Euler-angle-based solutions, we chose to implement a quaternion-based version of the optimal control problem introduced in Sec. II-B. Using such a representation requires to deal with numerical integration of quaternions which is a problem widely addressed in the literature [29], [30], [31]. Indeed, when representing orientations, quaternions must be unitary and thus belong to a constrained manifold (namely, the unit 3-sphere $S^3$). However, classical numerical integration schemes such as Runge-Kutta methods treat unit quaternions as if they were arbitrarily defined in $\mathbb{R}^4$. There exist at least two ways of addressing this problem: either perform a normalization step after each Runge-Kutta iteration to project the non-unitary quaternion onto $S^3$ (RKN) or using Crouch-Grossman methods which essentially use the exponential map of the Lie algebra instead of regular addition and multiplications. As shown in [32], for sufficiently small time steps, the performances of Crouch-Grossman and RKN methods are similar. We compared the RKN methods in our framework (i.e. time step imposed by the number of nodes of the OCP, and a constant rotation velocity ($\omega \approx 10$ rad/s) against an exact integration solving the first order differential equation $\dot{q} = \frac{1}{2}q \otimes \omega$. The RKN integration error was of the order of $1e^{-8}$ for 1 second of integration which is approximately the duration of our problem. In view of this acceptable error, we used RKN because it yields simpler relationships between the state variables and the control (integration of Eq. 3b) and facilitates the convergence of the optimization problem.

Since quaternions represent orientation and not rotation, the twist Mayer term of the objective function ($q_{rz}^r(T)$, in the Euler version of the OCP), eventually greater than one full rotation, was replaced by a discretized Lagrange term on the twisting velocity:

$$C_{QUAT} = \gamma \dot{q}_{ry}^r(0) + \sum_{i=1}^{nb_{\tau}} \beta_{q_i} \sum_{j=0}^{N} \tau_{ij}^2 + \alpha \sum_{j=0}^{N} \dot{q}_{rz}^r(T).$$

Both the quaternion and the Euler angles version of the problem were solved for the 8-DoF model, and did not show any change in the exploitation of the tilt angle of the root and the overall performance. Even if not constrained, the quaternion solution did not exceed $\pm \pi/6$ rad (same as for the Euler solution), with gimbal lock being at $\pm \pi/2$ rad.
<figure><img src="image_5.png" /><figcaption>Fig. 4: Mean (± 95% confidence interval) shoulder torques integrated over the full motion for 400 optimizations of a backward twisting somersault for the 2D model. Both the active and passive contributions of the shoulder joint are displayed for PHP and UHP techniques.</figcaption></figure>

The strategies differed slightly, which can be explained by the numerical discrepancies in the formulation of the two problems. One drawback of the quaternion approach is that it makes it impossible to directly express the cumulative number of rotations of the root (which is one of the Mayer terms of the cost function), therefore, to compute this quantity it was necessary to keep track of the instantaneous twisting velocity of the root and integrate it to yield its final orientation in Euler angles. In the rest of the study, the Euler approach was preferred merely because the resulting rotations are easier to interpret in terms of twisting performances.

### D. Multi-start approach

The first varying parameter of the multi-start approach was the number of shooting nodes ($N \in [295, 305] \subset \mathbb{N}$). This range of seemingly high shooting nodes numbers was chosen experimentally to obtain the best compromise between convergence rate and optimization time. The other parameters were related to the initial guess provided to the gradient-based solver (see Fig. 11 in appendix), which can yield different (sub-)optimal solutions:

* The twist rotation ($q_{rz}^r$) linearly increasing from 0 to $q_{rz}^r(T) \in [2, 5] \subset \mathbb{N}$,

* Uniform random arm elevation trajectories ($q_{el}^{ra}, q_{el}^{la} \in \mathcal{X}$),

* Uniformly distributed random arm torques $\underline{u}$ on $\mathcal{U}$

Each parameter was modified while the others were kept constant and the resulting grid of solutions was kept for results analysis (400 optimized solutions per model and per technique).

### E. Robustness analysis of the optimal solutions

Robustness of the optimized techniques was tested by adding 10° gaussian noise to the positions of the arms at the shooting nodes. The resulting noised state trajectory were then interpolated with quartic splines ([33]) in order to be differentiated and to compute the torques resulting from these "kinematics errors", using inverse dynamics. After applying the forward dynamics using these noised torques, the resulting number of twists were compared to the original ones,

<figure><img src="image_6.png" /><figcaption>Fig. 5: Mean (± 95% confidence interval) final number of twists generated by the optimal techniques according the model (2D vs 3D) and the objective function used in the optimization (PHP vs UHP). These statistics are computed on the best optimized solutions (leading to the highest number of twists) for each number of nodes. Therefore, each point represents the best-performing solution for a given model, a given objective type and a given number of nodes.</figcaption></figure>

yielding a level of robustness for each noised solution. Keeping track of these twist deviations allows to gain insight in the relevance of the optimal solutions. Indeed, ideal solutions are highly twisting and can be executed by athletes with different deviations from the optimal trajectory with approximately the same outcome in terms of number of twists (see Fig 7).

### F. Uncontrolled torques analysis

In the following, we highlight the fact that neglecting the change of plane of elevation at the shoulder is not a simplification from a motor control point of view. Indeed, numerically locking a DoF implies that no motion can be made around it, but it does not imply that the torques needed to meet this constraint are zero. In the case of human movement, every torques would result from muscle activations. In simulation, these uncontrolled torques are actually passively produced by the joint structure of the model, as opposed the controlled ones actively produced by the joint actuator. This is particularly the case in motions with fast rotations, as the expressions of torques in local frames attached to the joints include non-linear effects (gravity, Coriolis, centrifugal). Consequently, a motion restricted to one axis of a 2-DoF joint may still require joint torques on the two axes of the joint. Therefore, applying an active motor strategy from a simplified model to a model with a higher number of DoFs might lead to completely different whole-body behaviors. For this reason, we quantified the amount of passive joint torque at the shoulder during a backward twisting somersault on a 2D model.

Joints acceleration of the 2D model $\ddot{q}^{2D}$ were computed using the forward dynamics equation (Eq. 2). Then, the kinematics of the 2D model ($q^{2D}, \dot{q}^{2D}, \ddot{q}^{2D}$) was mapped to the 3D model ($q^{3D}, \dot{q}^{3D}, \ddot{q}^{3D}$) by imposing a null kinematics to the extra DoFs:

$$ \{ \dot{q}^{3D}[1..8], \dot{q}^{3D}[1..8], \ddot{q}^{3D}[1..8] \} = \{ q^{2D}, \dot{q}^{2D}, \ddot{q}^{2D} \}, \quad (6a) $$

$$ \{ \dot{q}^{3D}[9, 10], \dot{q}^{3D}[9, 10], \ddot{q}^{3D}[9, 10] \} = \{ \mathbf{0}, \mathbf{0}, \mathbf{0} \}, \quad (6b) $$
<figure><img src="image_7.png" /><figcaption>Fig. 6: Kernel density estimation representation of the distribution of twists generated by the optimal techniques in function of the model (2D vs 3D) and of the objective function used in the optimization (PHP vs UHP). These statistics are computed on every optimized solutions which yielded at least one twist. Note: the scales of the y-axis are different for the sake of readability.</figcaption></figure>

where $\dot{\boldsymbol{q}}^{3D} = [q_{tx}^r, q_{ty}^r, q_{tz}^r, q_{rx}^r, q_{ry}^r, q_{rz}^r, q_{el}^{ra}, q_{el}^{la}, q_{cp}^{ra}, q_{cp}^{la}]^T$. The passive torque absorbed by the structure of the joint in the 2D case can then be computed using the inverse dynamics function (Eq. 1), on the kinematics of the 3D model.

$$ \boldsymbol{\tau}^{3D} = M(\boldsymbol{q}^{3D})\ddot{\boldsymbol{q}}^{3D} + N(\boldsymbol{q}^{3D}, \dot{\boldsymbol{q}}^{3D}) + G(\boldsymbol{q}^{3D}) \quad (7a) $$

$$ = [\tau_{\text{root}}^{2D}, \tau_{\text{active}}^{2D}, \tau_{\text{passive}}^{2D}]^T \quad (7b) $$

Fig. 4 depicts the integral of active and passive torques produced/absorbed by the joint for the 2D model during a backward twisting somersault for 400 optimizations of both PHP and UHP techniques. Interestingly, the mean passive torque integrated during this motion is of the same order of magnitude as the active one. It is therefore not simpler, from a motor control point of view, to consider 2D kinetic strategies requiring torque control on 2 DoFs, rather than taking advantage of 3D kinetic strategies for the same level of control complexity. It is worth noticing that one consequence of penalizing the hand path (PHP vs UHP) is to limit the amount of both passive and active torques needed to execute the motion, leading to more relevant solutions (from a sports perspective), but less twists.

### III. RESULTS

Fig. 2 (resp. Fig. 3) is an animation of the most twisting solution obtained with the 2D (resp. 3D) model. This first result highlights that the best 3D solution (4.90 twists) is better than the best 2D solution (2.94 twists). The corresponding arms trajectories are displayed in Appendix (Figs. 8 and 9). The arms trajectories of another 3D solution yielding about the same number of twists (3.10) than the best 2D one is presented in Fig. 10 for comparison purposes. For a more systematic analysis, solutions generating the highest number of twist for each number of shooting nodes ($N \in [295, 305]$) were kept for further analysis (Fig. 5). Overall, good 3D techniques generated more twists than good 2D techniques whether penalized for their complexity or not, validating our hypothesis

<figure><img src="image_8.png" /><figcaption>Fig. 7: Logarithmic heatmaps of the twist deviation as a function of the optimal number of twist. Each optimal solution was randomly noised 500 times and the resulting performances (final number of twist) were computed. These performances were compared to the original one, giving a number of twist deviations (positive is more). These computations are performed on every optimized solutions which yielded at least one twist.</figcaption></figure>

about the improvement brought by the addition of an extra DoF at each shoulder. Optimized solutions convergence rate, mean optimization duration and mean number of twists are presented in Tab. III. These results show that the 3D problem is harder to solve than the 2D one (higher convergence times and lower convergence rates) but yields better performances in terms of twists. Although lower than in the 2D case, 3D convergence rates remain reasonably high to be used in a multi-start approach. Adding a penalization on the hand path (PHP vs UHP) also complicated the problem (see Tab. III) and yielded lower number of twist but more realistic solutions (less arm motion). This result supports the relevance of our custom cost function Eq. (5) which was designed to find shorter yet more efficient arm trajectories for twisting the body.

<table><caption>TABLE III: Optimization statistics for each model and each type of cost function. The convergence rate and the mean time of convergence are given for all the optimized solution (400 per line). The number of twist is given for the best twisting solutions per number of nodes, as in Fig. 5.</caption><thead><tr><th>Model</th><th>Cost</th><th>Convergence</th><th>Time (s)</th><th>Twists ± (95% ci)</th></tr></thead><tbody><tr><td>2D</td><td>UHP</td><td>99.1%</td><td>294</td><td>1.9 ± 0.5</td></tr><tr><td>3D</td><td>UHP</td><td>73.9%</td><td>1214</td><td>4.5 ± 0.3</td></tr><tr><td>2D</td><td>PHP</td><td>97.7%</td><td>387</td><td>1.4 ± 0.1</td></tr><tr><td>3D</td><td>PHP</td><td>69.5%</td><td>1226</td><td>3.6 ± 0.4</td></tr></tbody></table>

Fig. 6 displays a kernel density estimation of the distribution of twists generated by the optimal techniques in function of the model (2D vs 3D) and of the objective function used in the optimization (PHP vs UHP). First, this representation allows to note that there are clusters of solutions, which correspond to local minima of the constrained objective function. These clusters are narrower with the 2D model regardless of the hand path penalization. Secondly, regardless of the model, the most performing UHP techniques led to higher number of twists (2D: 1.6, 3D: 4.7) than the most performing PHP ones (2D: 1.3, 3D: 4), with more clusters of near-optimal solutions.
Finally, as already revealed, with the 3D model, there exist clusters of solutions yielding higher number of twists than the best clusters of the 2D model (see also Fig. 5). It is worth noticing the really low density of 2D solution yielding more than 1.75 twists, illustrating the importance of multistart via the impact of the initial guess on the optimized solution.

Fig. 7 presents logarithmic heatmaps of the twist deviation due to the introduced kinematic noise, as a function of the number of twist generated by the initial solution. The higher the heatmaps, the higher the number of solutions in the underlying discretized area. The main result from this representation is the existence of robust clusters of solutions, corresponding to warm-colored cells on the 0-twist-deviation horizontal line. For instance, for the 2D model, Fig. 7 shows a high concentration of solutions initially producing around 1 twist which remain stable (around 0 twist deviation). For the same model, there also exist solutions with higher initial number of twist (2.4), of which a substantial part remains stable. For the 3D model, looking at the 0-twist-deviation horizontal line reveals high-twisting clusters of stable solutions (2.2, 3, 3.5 and 4 twists).

## IV. DISCUSSION

The main objective of this work was to study optimal 3D arm strategies for maximizing twist rotation during the single back layout somersault of a multibody model.

Indeed, in acrobatic sports, it is commonly admitted that twists in backward somersaults are generated mainly with 2D motion of the arm [34]. Overall, considering change of plane of abduction/adduction at the shoulder significantly increased the twisting performance, by about two twists compared to abduction/adduction only. Apart from pure performance, another interesting outcome of this investigation lies in the comparison between 2D and 3D strategies yielding similar number of twists (about 3; see Figs. 8 and 10). It reveals that exploiting 3D arms trajectories can result in completely different solutions (less arm abd/abduction, compensated for by elevation change of plane) for a similar final performance.

Thus, our results temper the consensus about 2D arms strategies for twisting somersaults. We suggest that athletes would gain in performance if they used 3D motion in a proper manner. Moreover, studies maximizing twists with 2D models disregard the dimensionality of the torques needed to execute these planar trajectories [5], [3]. In the present work, it is shown that these torques are not negligible, but of the same order as the controlled ones. For all this reasons, we believe that adding an extra DoF at each shoulder is a relevant modeling approach.

As shown in this study, the optimal control formulation (Eq. 3) gives a wide flexibility in the choice of the cost function, allowing one to maximize the performance of the system while penalizing too complex and thus unrealistic solutions. For instance, the minimization of segment trajectory is used in sports such as barbell lifting [35], where the underlying goal is to minimize the work done by the athlete. In the present study, we have shown its relevance in artistic sports to control movement complexity. Indeed, two cost

functions (PHP and UHP) were implemented and the analysis of the results reveals that, as expected, when penalized in complexity, the optimal solutions resulted in 35% (resp. 20%) less twists with the 2D (resp. 3D) model (Fig. 5), but were less demanding in terms of joints torques (Fig. 4).

The multiple shooting transcription of the problem combined with the non-linear solver IPOPT was a relevant choice as it enabled us to develop a multi-start approach with poor initial guesses and a lot of room for the solver to find innovative solutions. We believe that this strategy was particularly interesting as it resulted in different clusters of near-optimal solutions (see Fig. 6), which will be qualitatively studied and explained in a near-future study. The existence of very low-density clusters of solutions illustrates the relevance of the multistart approach because it highlights the impact of the initial guess on the optimized solution. When looking for innovative motor strategies, as in sports sciences for instance, this kind of approach should be preferred as it enables a minimal amount of a priori knowledge (from the coaches, the community, etc...) to be fed into the solver, leaving random starts and gradient-based descents the chance to discover new techniques [20]. Of course, a direct drawback of such an implementation is its computational cost which is constantly reducing thanks to the advances in numerical optimization. In particular, the number of shooting nodes used for discretizing the optimization problem was relatively important because it reduces the integration time between two nodes and thus limits the risk of reaching a singularity on the root orientation during integration.

As pointed out already, techniques needing infinite precision to achieve targeted performances are irrelevant to any human-related application scenario (particularly acrobatic sports) because of the variability inherent to human motion [22]. In this regard, robustness of the optimized solution was studied after optimization to determine whether errors during the execution of the solution would lead to noticeable changes in its performance. Our choice was to introduce the errors at a kinematic level, which could be discussed from a motor control point of view, but remains a relevant test for robustness. The results show that, if a part of the solutions are not robust to the noise introduced (see dispersion in Fig. 7), there exist clusters of stable solutions. Interestingly, some of the highly-twisting solutions (3 or more twists) optimized with the 3D model are robust to kinematic noise, demonstrating a clear advantage over 2D solutions whose stable cluster generates only a bit more than one twist. This is a promising result since an ideal solution would score a high performance while being robust to execution errors. Before concluding on the transferability of such techniques to the sports community however, a more in-depth analysis of the robustness should be conducted (e.g. violation of the problem constraints after adding noise). The fact remains that, according to the results obtained in the present work, considering 3D trajectories of the shoulder for twisting somersault is a better choice on all fronts. Moreover, it was numerically demonstrated that increasing the kinematic fidelity of the model does not increase the level of control complexity on the real system, which is another
argument in favor of using the 3D model.

Another point concerns the use of quaternions for modeling
joints rotations in a context of numerical optimal control.
This was shown to be feasible, with no implications on the
nature of the solution when applied at the root segment of
the model. Although our model permits 3D arm motions, it
remains relatively simple and improvable especially by adding
hip or thorax joints known to be used in optimal techniques
[6]. In future works, it would also be relevant to add one
flexion at each elbow, but this would require to add another
rotation at the shoulder to control the full 3D orientation of the
arm. Therefore, Gimbal lock problem could occur with Euler
angles representation, a situation which could be avoided using
a quaternion representation of the orientation at the shoulder
[36]. Also, although self-collision of the arms with the body
where avoided thanks to kinematic constraints, future work
could involve general non-crossing segment constraints, in
order to avoid arms collisions for instance.

In conclusion, 3D arms techniques during backward
somersault are relevant since they produce more twist
than 2D ones without increasing the complexity of the
underlying motor control. The contributions of this work
should be of interest to acrobatic sports coaches and to the
biomechanics community, as the presented approach (large
batch optimization, selection of the solutions, robustness
analysis) could be generalized to other motion simulations.

ACKNOWLEDGMENT

This work was supported by the Mitacs Acceleration grant
#IT11976 and by the Institut National du Sport du Québec.

REFERENCES

[1] M. R. Yeadon and M. J. Hiley, “The control of twisting somersaults,” *Journal of biomechanics*, vol. 47, no. 6, pp. 1340–1347, 2014.

[2] J. Koschorreck and K. Mombaur, “Modeling and optimal control of human platform diving with somersaults and twists,” *Optimization and Engineering*, vol. 13, no. 1, pp. 29–56, 2012.

[3] H. R. Dullin and W. Tong, “Twisting somersault,” *SIAM Journal on Applied Dynamical Systems*, vol. 15, no. 4, pp. 1806–1822, 2016.

[4] C. Frohlich, “The physics of somersaulting and twisting,” *Scientific American*, vol. 242, no. 3, pp. 154–165, 1980.

[5] S. Bharadwaj, N. Duignan, H. R. Dullin, K. Leung, and W. Tong, “The diver with a rotor,” *Indagationes Mathematicae*, vol. 27, no. 5, pp. 1147–1161, 2016.

[6] M. R. Yeadon and M. J. Hiley, “Twist limits for late twisting double somersaults on trampoline,” *Journal of biomechanics*, vol. 58, pp. 174–178, 2017.

[7] P. Pigeon, S. B. Bortolami, P. DiZio, and J. R. Lackner, “Coordinated turn-and-reach movements. i. anticipatory compensation for self-generated coriolis and interaction torques,” *Journal of neurophysiology*, vol. 89, no. 1, pp. 276–289, 2003.

[8] A. J. van den Bogert, D. Blana, and D. Heinrich, “Implicit methods for efficient musculoskeletal simulation and optimal control,” *Procedia Iutam*, vol. 2, pp. 297–316, 2011.

[9] F. De Groote, A. L. Kinney, A. V. Rao, and B. J. Fregly, “Evaluation of direct collocation optimal control problem formulations for solving the muscle redundancy problem,” *Annals of biomedical engineering*, vol. 44, no. 10, pp. 2922–2936, 2016.

[10] M. R. Yeadon and M. J. Hiley, “The mechanics of the backward giant circle on the high bar,” *Human Movement Science*, vol. 19, no. 2, pp. 153–173, 2000.

[11] F. Leboeuf, G. Bessonnet, P. Seguin, and P. Lacouture, “Energetic versus sthenic optimality criteria for gymnastic movement synthesis,” *Multibody System Dynamics*, vol. 16, no. 3, pp. 213–236, 2006.

[12] N. V. R. K. N. Murthy and S. S. Keerthi, “Optimal control of a somersaulting platform diver: A numerical approach,” in [*1993 Proceedings IEEE International Conference on Robotics and Automation*. IEEE, 1993, pp. 1013–1018.

[13] M. Koh, L. Jennings, B. Elliott, and D. Lloyd, “A predicted optimal performance of the yurchenko layout vault in women’s artistic gymnastics,” *Journal of applied biomechanics*, vol. 19, no. 3, pp. 187–204, 2003.

[14] M. R. Yeadon, “The simulation of aerial movement—ii. a mathematical inertia model of the human body,” *Journal of biomechanics*, vol. 23, no. 1, pp. 67–74, 1990.

[15] R. R. Neptune, D. J. Clark, and S. A. Kautz, “Modular control of human walking: a simulation study,” *Journal of biomechanics*, vol. 42, no. 9, pp. 1282–1287, 2009.

[16] J. Higginson, R. R. Neptune, and F. Anderson, “Simulated parallel annealing within a neighborhood for optimization of biomechanical systems,” *Journal of biomechanics*, vol. 38, no. 9, pp. 1938–1942, 2005.

[17] J. F. Schutte, B.-I. Koh, J. A. Reinbolt, R. T. Haftka, A. D. George, and B. J. Fregly, “Evaluation of a particle swarm algorithm for biomechanical optimization,” *Journal of Biomechanical Engineering*, vol. 127, no. 3, pp. 465–474, 2005.

[18] L.-F. Lee and B. R. Umberger, “Generating optimal control simulations of musculoskeletal movement using opensim and matlab,” *PeerJ*, vol. 4, p. e1638, 2016.

[19] S. J. Wright, “Interior point methods for optimal control of discrete time systems,” *Journal of Optimization Theory and Applications*, vol. 77, no. 1, pp. 161–187, 1993.

[20] A. Huchez, D. Haering, P. Holvoët, F. Barbier, and M. Begon, “Local versus global optimal sports techniques in a group of athletes,” *Computer methods in biomechanics and biomedical engineering*, vol. 18, no. 8, pp. 829–838, 2015.

[21] A. Huchez, D. Haering, P. Holvoet, F. Barbier, and M. Begon, “Differences between expert and novice gymnasts performance of a counter movement forward in flight on uneven bars.” *Science of Gymnastics Journal*, vol. 8, no. 1, 2016.

[22] M. J. Hiley and M. R. Yeadon, “Investigating optimal technique in a noisy environment: application to the upstart on uneven bars,” *Human Movement Science*, vol. 32, no. 1, pp. 181–191, 2013.

[23] B. Michaud and M. Begon, “Biorbd: Toolbox for biomechanical analyses,” Web page, 2018. [Online]. Available: https://github.com/pyomeca/biorbd

[24] M. L. Felis, “Rbdl: an efficient rigid-body dynamics library using recursive algorithms,” *Autonomous Robots*, pp. 1–17, 2016. [Online]. Available: http://dx.doi.org/10.1007/s10514-016-9574-0

[25] R. Featherstone, *Rigid body dynamics algorithms*. Springer, 2014.

[26] J. A. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl, “Casadi: a software framework for nonlinear optimization and optimal control,” *Mathematical Programming Computation*, vol. 11, no. 1, pp. 1–36, 2019.

[27] A. Wächter and L. T. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” *Mathematical programming*, vol. 106, no. 1, pp. 25–57, 2006.

[28] S. Namdari, G. Yagnik, D. D. Ebaugh, S. Nagda, M. L. Ramsey, G. R. Williams Jr, and S. Mehta, “Defining functional shoulder range of motion for activities of daily living,” *Journal of shoulder and elbow surgery*, vol. 21, no. 9, pp. 1177–1183, 2012.

[29] P. E. Crouch and R. Grossman, “Numerical integration of ordinary differential equations on manifolds,” *Journal of Nonlinear Science*, vol. 3, no. 1, pp. 1–33, 1993.

[30] M. Boyle, “The integration of angular velocity,” *Advances in Applied Clifford Algebras*, vol. 27, no. 3, pp. 2345–2374, 2017.

[31] L. Seelen, J. Padding, and J. Kuipers, “Improved quaternion-based integration scheme for rigid body motion,” *Acta Mechanica*, vol. 227, no. 12, pp. 3381–3389, 2016.

[32] M. S. Andrle and J. L. Crassidis, “Geometric integration of quaternions,” *Journal of Guidance, Control, and Dynamics*, vol. 36, no. 6, pp. 1762–1767, 2013.

[33] G. Bessonnet, P. Seguin, and P. Sardain, “A parametric optimization approach to walking pattern synthesis,” *The International Journal of Robotics Research*, vol. 24, no. 7, pp. 523–536, 2005.

[34] M. R. Yeadon, *Twisting somersaults*. SB & MC, 2015.

[35] S. L. Nejadian, M. Rostami, and F. Towhidkhah, “Optimization of barbell trajectory during the snatch lift technique by using optimal control theory,” *American Journal of Applied Sciences*, vol. 5, no. 5, pp. 524–531, 2008.
<footnote>[36] C. Kuhlmann, T. L. Milani, et al., "Investigation of shoulder kinematics in volleyball spikes," in *ISBS-Conference Proceedings Archive*, 2008.</footnote>
<table><caption>TABLE IV: Summary table of the kinematics and kinetics bounds enforced by the optimization problem at the start, during the skill and at the end of it. — stands for infinite bounds</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="6">2D model</th><th colspan="2">3D model</th></tr><tr><th>q<sup>r</sup><sub>tx</sub></th><th>q<sup>r</sup><sub>ty</sub></th><th>q<sup>r</sup><sub>tz</sub></th><th>q<sup>r</sup><sub>rx</sub></th><th>q<sup>r</sup><sub>ry</sub></th><th>q<sup>r</sup><sub>rz</sub></th><th>q<sup>ra</sup><sub>el</sub> &amp; q<sup>la</sup><sub>el</sub></th><th>q<sup>ra</sup><sub>cp</sub> &amp; q<sup>la</sup><sub>cp</sub></th></tr></thead><tbody><tr><td rowspan="3">start</td><td>position (° or m)</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>[160, 180]</td><td>[-135, 45]</td></tr><tr><td>velocity (°/s or m/s)</td><td>[-10, 10]</td><td>[-10, 10]</td><td>[5.89, 6.09]</td><td>0</td><td>—</td><td>0</td><td>[-5730, 5730]</td><td>[-5730, 5730]</td></tr><tr><td>torque (Nm)</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>[-100, 100]</td><td>[-100, 100]</td></tr><tr><td rowspan="3">skill</td><td>position (° or m)</td><td>—</td><td>—</td><td>—</td><td>[-45, 45]</td><td>—</td><td>—</td><td>[0, 180]</td><td>[-135, 45]</td></tr><tr><td>velocity (°/s or m/s)</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>[-5730, 5730]</td><td>[-5730, 5730]</td></tr><tr><td>torque (Nm)</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>[-100, 100]</td><td>[-100, 100]</td></tr><tr><td rowspan="3">end</td><td>position (° or m)</td><td>[-0.1, 0.1]</td><td>[-0.1, 0.1]</td><td>[-0.1, 0.1]</td><td>[-15, 15]</td><td>[350, 370]</td><td>—</td><td>[155, 165]</td><td>[-135, 45]</td></tr><tr><td>velocity (°/s or m/s)</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>[-5730, 5730]</td><td>[-5730, 5730]</td></tr><tr><td>torque (Nm)</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>[-100, 100]</td><td>[-100, 100]</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Fig. 8: Right and left arm abd/adduction (elevation) of the most-twisting 2D solution (2.94 twits), without penalization on the hand trajectory.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9: Right and left arm abd/adduction (elevation) and elevation change of plane of the most-twisting 3D solution (4.90 twits), without penalization on the hand trajectory.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10: Right and left arm abd/adduction (elevation) and elevation change of plane of another 3D solution (3.10 twits), with penalization on the hand trajectory.</figcaption></figure>

<table><caption>Fig. 11: Initial guesses of each optimization variable, for the 3D and 2D models. Empty field means initialization to 0.</caption><thead><tr><th>variable</th><th>initial guess</th></tr></thead><tbody><tr><td><i>T</i></td><td>1.2s</td></tr><tr><td><i>u</i></td><td>uniformly distributed random arm torques in [-10, 10] Nm</td></tr><tr><td><i>q</i><sup>r</sup><sub>tx</sub>, <i>q</i><sup>r</sup><sub>ty</sub></td><td></td></tr><tr><td><i>q</i><sup>r</sup><sub>tz</sub></td><td>free fall parabolic trajectory for initial velocity of 6.0m/s</td></tr><tr><td><i>q</i><sup>r</sup><sub>rz</sub></td><td>free fall affine trajectory for initial velocity of 6.0m/s</td></tr><tr><td><i>q</i><sup>r</sup><sub>rx</sub></td><td></td></tr><tr><td><i>q</i><sup>r</sup><sub>ry</sub></td><td>linear increase from 0° to 360°</td></tr><tr><td><i>q</i><sup>r</sup><sub>ry</sub></td><td>300°/s</td></tr><tr><td><i>q</i><sup>r</sup><sub>rz</sub></td><td>linear increase from 0° to <i>q</i><sup>r</sup><sub>rz</sub>(T)</td></tr><tr><td><i>q</i><sup>r</sup><sub>rz</sub></td><td><i>q</i><sup>r</sup><sub>rz</sub>(T)/T</td></tr><tr><td><i>q</i><sup>ra</sup><sub>el</sub>, <i>q</i><sup>la</sup><sub>el</sub></td><td>uniformly distributed random arm elevation in [0, 180] °</td></tr><tr><td><i>q</i><sup>ra</sup><sub>cp</sub>, <i>q</i><sup>la</sup><sub>cp</sub></td><td></td></tr><tr><td><i>q</i><sup>ra</sup><sub>cp</sub>, <i>q</i><sup>la</sup><sub>cp</sub></td><td></td></tr></tbody></table>