# Derivative-Free Optimization Approaches for Force Polytopes Prediction

Gautier Laisné, Jean-Marc Salotti, Nasser Rezzoug

► To cite this version:

Gautier Laisné, Jean-Marc Salotti, Nasser Rezzoug. Derivative-Free Optimization Approaches for Force Polytopes Prediction. ESANN 2023 - European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning, Oct 2023, Bruges, Belgium. pp.339-344, 10.14428/esann/2023.ES2023-122. hal-04330766

HAL Id: hal-04330766

https://hal.science/hal-04330766v1

Submitted on 8 Dec 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.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Derivative-Free Optimization Approaches for Force Polytopes Prediction

G. Laisné¹, J-M. Salotti¹,² and N. Rezzoug¹,³ *

1 - INRIA Center at University of Bordeaux
Talence - France

2 - University of Bordeaux, CNRS, Bordeaux INP, IMS, UMR 5218, F-33400
Talence - France

3 - PPrime Institute, CNRS, University of Poitiers, ENSMA, UPR 3346
Poitiers - France

**Abstract.** Hand force capacities reflect an individual's ability to generate forces in all directions, considering a given upper-limb posture. These capacities are described as polytopes by means of an upper-limb musculoskeletal model. However, such a model needs to be adapted to an individual for more accuracy. The model parameter space is investigated using derivative-free algorithms which do not require the optimization function to be differentiable: genetic algorithms and SRACOS, a classification-based algorithm. Results demonstrate that employing a genetic algorithm with a polytope representation in 26 vertices yields the most accurate prediction of force capacities in a validation posture.

## 1 Introduction

In biomechanics, human hand force capacities refer to the set of feasible forces exertable at the hand considering arm posture and muscle tensions. It sees application in physical Human-Robot Interaction (pHRI), in which an interaction between the operator's hand and the robot's end-effector is considered. Knowing the operator's force capacities allows the robot to adjust its assistance to avoid exceeding force limits [1]. However, these limits depend on the operator's posture, muscle geometry and biomechanical properties.

The human arm can be described by means of a musculoskeletal (MSK) model formalism comprised of a set of bodies linked through joints and muscles producing joint torques. Force capacities are then represented as a 3D convex bounded polytope, called the *force polytope* (sections 2.1 to 2.3).

Supervised training of force polytopes on chosen postures can be employed for prediction across different postures. In [4], a neural network is used, although MSK model parameters are assumed, limiting predictions to a specific subject.

In this work, the arm and muscles geometry are known, unlike the muscle biomechanical parameters. Only 3 postures are considered to accommodate to the experimental difficulties of gathering force capacities. A non-differentiable function to optimize is described in section 2.4 to fit MSK model parameters to

<footnote>*The simulations in this paper were carried out using the PlaFRIM experimental testbed (Inria, CNRS, LABRI, IMB, Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine) (see https://www.plafrim.fr).</footnote>
reference force polytopes. Two solvers are considered: a genetic algorithm and SRACOS, a classification-based algorithm [3] (section 2.5). A hyperparameter space is then defined to explore 4 fitting postures set and 3 polytope representations. Each instance is repeated using 5 different initial populations (section 2.6). The results section evaluates the prediction quality of the best MSK model parameters. Finally, results are discussed in the conclusion.

# 2 Methods

## 2.1 Polytope

A *polytope* is the generalization of a polygon in higher dimensions. It is assumed to be convex and bounded, so that it can be described by its vertices [2]. In this work, a polytope *P* is constructed from the intersection of a *zonotope* *Z* (the projection of an hypercube) with an affine subspace *A* [2] (Fig. 1).

<figure><img src="image_2.png" /><figcaption>Fig. 1: A 2D polytope *P* created from a 3D zonotope *Z* and a plane *A*.</figcaption></figure>

Enumerating the vertices of *P* is impractical, due to the high combinatorics involved in the zonotope computation [2]. The Iterative Convex Hull (ICH) method allows a reasonably fast approximation of the vertices within a chosen tolerance [1]. Using a 2x8-cores (2.61GHz) Intel i9-11950H processor, the ICH (with a tolerance of 1N) mean time for 100 randomly generated 3D affine subspaces and 7D zonotopes built from 50D hypercubes is $1.03 \pm 0.05$ seconds.

## 2.2 Force polytope

In biomechanics, force capacities correspond to all possible forces that can be exerted in every direction by an individual at a specific position (for example, at the hand). Experimentally, measuring force capacities is long and demanding, due to the challenge of achieving maximal force in specific directions. However, using a musculoskeletal model representing a human upper-limb, force capacities can be described as a 3D polytope called the *force polytope* [4]. The shape and volume vary according to the individual's posture.

The upper-limb is considered as a *n* degree-of-freedoms kinematic chain actuated by *m* muscles. In static conditions, the force polytope localized at the center of the hand for a posture **x** ∈ ℝⁿ is the convex 3D polytope defined as $P(\mathbf{x}) = \{\mathbf{f} \in \mathbb{R}^3 \mid J^T(\mathbf{x})\mathbf{f} = -L^T(\mathbf{x})\mathbf{t}, \mathbf{t} \in [\mathbf{t}_{\min}(\mathbf{x}), \mathbf{t}_{\max}(\mathbf{x})]\},$ where $J(\mathbf{x}) = \frac{\partial X}{\partial \mathbf{x}} \in \mathbb{R}^{3 \times n}$ is the cartesian end-effector jacobian positionned at X, $L(\mathbf{x}) = \frac{\partial L(\mathbf{x})}{\partial \mathbf{x}} \in \mathbb{R}^{m \times n}$ is the muscle-length jacobian of the system and $\mathbf{t} \in \mathbb{R}^m$ are the possible muscle tensions with $\mathbf{t} \in [\mathbf{t}_{\min}, \mathbf{t}_{\max}]$.
## 2.3 Polytope proximity assessment

Assessing proximity of polytopes is non-trivial due to their geometric nature [2]. Three representations of P are considered and noted $R_1^P$, $R_2^P$ and $R_3^P$ (Fig. 2).

<figure><img src="image_3.png" /><figcaption>Fig. 2: The polytope representations in 2D. For $R_1$, all the vertices are considered. For $R_2$, the vertices result from intersecting P with lines oriented of angles $0, \frac{\pi}{4}, \frac{\pi}{2}$ and $\frac{3\pi}{4}$. For $R_3$, the eigenvectors $u_1$ and $u_2$ of P are computed and vertices are created by intersecting them with P.</figcaption></figure>

* $R_1$ corresponds to the vertices computed via the ICH algorithm.

* $R_2$ consists of the norms of 26 vertices created by intersecting P with 13 lines defined to span the 3D space in diverse directions. Each line induces 2 vertices: one in the same direction, noted $v^+$, and the other one noted $v^-$. Thus, $R_2 := (\|v_{0,-\frac{\pi}{4}}^+\|, \|v_{0,-\frac{\pi}{4}}^-\|, \dots, \|v_{0,\frac{\pi}{2}}^+\|, \|v_{0,\frac{\pi}{2}}^-\|) \in \mathbb{R}^{26}$ such that for all spherical azimuth and polar angles $(\theta, \varphi) \in (\{0, \frac{\pi}{4}, \frac{\pi}{2}, \frac{3\pi}{4}\} \times \{-\frac{\pi}{4}, 0, \frac{\pi}{4}\}) \cup (0, \frac{\pi}{2}), v_{\theta,\varphi}^+$ and $v_{\theta,\varphi}^-$ are the vertices obtained from $P \cap L(\theta, \varphi)$ where $L(\theta, \varphi)$ is the line in $\mathbb{R}^3$ oriented by angles $(\theta, \varphi)$.

* $R_3(u_1, u_2, u_3)$ is similar to $R_2$, except with 3 lines corresponding to the covariance matrix' eigenvectors generated by another polytope $P'$, leading to 6 vertices. These lines define the ellipsoid axes fitting the vertices of $P'$.

To compare polytopes, functions $f_1, f_2$ and $f_3$ are defined as follows:

* $f_1$ is the Hausdorff distance between two polytopes $P_1$ and $P_2$. For vertices $v_1 \in R_1^{P_1}$ and $v_2 \in R_1^{P_2}$, it is defined as:

$$f_1(P_1, P_2) = \max\{\sup_{v_1} \inf_{v_2} \sqrt{\|v_1 - v_2\|}, \sup_{v_2} \inf_{v_1} \sqrt{\|v_1 - v_2\|}\}$$

* $f_2$ returns the mean of all distances between each pair of vertices being on the same line and with the same direction: $f_2(P_1, P_2) = \text{MSE}(R_2^{P_1}, R_2^{P_2})$ with MSE being the mean squared error.

* $f_3$ is similar to $f_2$ except for the polytopes representation: $f_3(P_1, P_2) = \text{MSE}(R_3^{P_1}(u_1, u_2, u_3), R_3^{P_2}(u_1, u_2, u_3))$ with $u_1, u_2$ and $u_3$ the covariance matrix' eigenvectors of $P_1$.
## 2.4 Optimization problem

Let $\theta \in \mathbb{R}^p$ be a vector of $p$ parameters of a MSK model and consider $S$ to be a finite set of postures in $\mathbb{R}^n$. Our goal is to find a parameter vector $\theta$ such that its generated force polytopes are close enough to the force polytopes generated by the parameters $\theta_0$, according to a comparison function $f$.

$$ \theta^* = \arg \min_{\theta \in \mathbb{R}^p} \frac{1}{|S|} \sum_{\mathbf{x} \in S} f(P_{\theta_0}(\mathbf{x}), P_{\theta}(\mathbf{x})) $$

The prediction quality of $\theta^*$ is measured by $f(P_{\theta^*}(\mathbf{x}_{val}))$, for $\mathbf{x}_{val}$ a posture not in $S$. The smaller is the value, the better is the quality.

## 2.5 Derivative-free optimization

The defined optimization involves intersections, which are not linearizable nor differentiable so that non-gradient methods should be used as solvers.

A *derivative-free optimization* learns a model distribution in the search space and uses sampling to understand the function topology, allowing local researches to give better solutions, without assuming the function to be differentiable.

In genetic algorithms, the model is a set of good solutions and the sampling is done through variations on solutions. Our implementation uses the following strategy: the 5 best solutions minimizing the cost function are kept for the next generation and mutations are derived from them. The mutation process generates 5 new solutions randomly selected in the neighbourhood of the parent parameter set. The neighbourhood size varies depending on the number of times a solution has been selected as a parent: below 4 times, the neighbourhood includes up to 10% around the parent solution. Between 5 and 9 times, up to 1%. Between 10 and 29, up to 0.5% and above 30 times, up to 0.1%.

In SRACOS, the model is a hypercube and the sampling is from the uniform distribution in the hypercube. It learns to classify solutions as either *positive* or *negative*, using reinforcement learning described as a Markov decision process [3]. The solutions are then sampled from the *positive* areas. SRACOS is well-suited for optimization on non-differentiable and non-convex functions.

## 2.6 Hyperparameters

The hyperparameter space $\mathcal{H}$ includes 12 elements ($S, R$) such that:

* $S \in \{S_1, S_2, S_3, S_4\}$ is a set of postures such that $S_1 = \{\mathbf{x}_1, \mathbf{x}_2, \mathbf{x}_3\}$, $S_2 = \{\mathbf{x}_1, \mathbf{x}_2, \mathbf{x}_4\}$, $S_3 = \{\mathbf{x}_2, \mathbf{x}_3, \mathbf{x}_4\}$ and $S_4 = \{\mathbf{x}_1, \mathbf{x}_3, \mathbf{x}_4\}$. The posture $\mathbf{x}_i \in \mathbb{R}^7$ is defined using Stanford's model in the global reference frame (the $x$-axis is normal to the coronal plane, the $y$-axis normal to the transverse plane and the $z$-axis normal to the sagittal plane). They are described by 7 Euler angles in degrees using a $y-z'-y''$ sequence for the shoulder, a $z-y'$ sequence for the elbow and a $x-z'$ sequence for the wrist: $\mathbf{x}_1 = (31^\circ, 12^\circ, -34^\circ, 74^\circ, 17^\circ, 15^\circ, 15^\circ)$, $\mathbf{x}_2 = (79^\circ, 54^\circ, -73^\circ, 53^\circ, 90^\circ, 0^\circ, 0^\circ)$, $\mathbf{x}_3 = (124^\circ, 56^\circ, -75^\circ, 32^\circ, 88^\circ, 0^\circ, 0^\circ)$ and $\mathbf{x}_4 = (27^\circ, 66^\circ, 59^\circ, 34^\circ, 88^\circ, 0^\circ, 0^\circ)$;
• $R \in \{R_1, R_2, R_3\}$ is the polytope representation to use (cf. section 2.3).

The quality of the best solution found for given hyperparameters ($S, R$) is computed for the posture not used in the postures set $S$. We denote this posture $\mathbf{x}_{val}^S = \{\mathbf{x}_1, \mathbf{x}_2, \mathbf{x}_3, \mathbf{x}_4\} \setminus S$.

## 2.7 Datasets

Polytopes are computed based on Stanford's upper-limb MSK model, which uses a Hill-based muscle model [5].

A solution $\theta \in \mathbb{R}^{150}$ parametrizes Stanford's MSK model via 3 parameters per muscle (the maximum isometric force, the optimal length at which the muscle creates its maximal isometric force and the tendon slack length). $\theta_0$ denotes the model parameter set initially parametrizing the Stanford's model.

# 3 Experimental results

In the series of tests, each algorithm employs a specific solving method (either the genetic algorithm or SRACOS) and is initialized with a set of hyperparameters ($S, R$) across 5 different initial generations, leading to 120 instances distributed over 2x64-core AMD Zen3 machines. All scripts are written in Python 3.8 and the packages `biorbd`, `pycapacity` 1.2.19 and `ZOOpt` 0.3.0 have been used to describe the musculoskeletal model, the ICH method and the SRACOS implementation. The ICH method is set to a tolerance of 1N [1]. All algorithms stop after exploring 32000 solutions.

In average for all initial generations, the best polytope prediction comes from the solution found by the genetic algorithm for postures set $S_3$ and the polytope representation $R_2$ (Fig. 3 and Table 1). SRACOS tends to overfit more than genetic algorithms.

<figure><img src="image_4.png" /><figcaption>Fig. 3: (a) Polytope to estimate, generated by $\theta_0$ at validation posture $\mathbf{x}_{val}^{S_3}$. (b) Polytope predicted by the solution $\theta^*$ found by the genetic algorithm with hyperparameters ($S_3, R_2$) at posture $\mathbf{x}_{val}^{S_3}$ during the first initial population.</figcaption></figure>
<table><caption>Table 1: White lines are the results for genetic algorithms and gray's are for SRACOS. For postures set $S_3$, genetic algorithms using polytope representation $R_2$ lead, in average, to the best prediction quality on validation posture $x_1$.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Evaluation of f<br>on fitting postures in S<sub>3</sub></th><th colspan="3">Prediction quality on<br>validation posture x<sup>S<sub>3</sub></sup><sub>val</sub> = x<sub>1</sub></th></tr><tr><th>x<sub>2</sub></th><th>x<sub>3</sub></th><th>x<sub>4</sub></th><th>f<sub>1</sub></th><th>f<sub>2</sub></th><th>f<sub>3</sub></th></tr></thead><tbody><tr><td>R<sub>1</sub></td><td>271±21</td><td>*</td><td>*</td><td>115±13</td><td>*</td><td>*</td></tr><tr><td>R<sub>2</sub></td><td>398±30</td><td>81±45</td><td>875±938</td><td><b>125±34</b></td><td><b>252±170</b></td><td><b>157±83</b></td></tr><tr><td>R<sub>3</sub></td><td>403±25</td><td>456±288</td><td>100±36</td><td>150±18</td><td>*</td><td>*</td></tr><tr><td>R<sub>1</sub></td><td>130±119</td><td>*</td><td>*</td><td>81±75</td><td>*</td><td>*</td></tr><tr><td>R<sub>2</sub></td><td>346±50</td><td>13±8</td><td>935±968</td><td>*</td><td>*</td><td>*</td></tr><tr><td>R<sub>3</sub></td><td>403±75</td><td>432±270</td><td>13±5</td><td>129±25</td><td>181±117</td><td>569±735</td></tr></tbody></table>

* means that the mean value and the standard deviation are greater than 1000.

## 4 Conclusion

Force capacities were represented as 3D polytopes, using a musculoskeletal model parametrized by $\theta \in \mathbb{R}^{150}$. Model-fitting was described as a non-differentiable optimization problem to find a solution $\theta^*$ generating force polytopes similar to those created by a reference parameter set $\theta_0$. Two derivative-free algorithms, a genetic algorithm and SRACOS, were used as solvers. 4 sets of postures and 3 polytope representations were considered. Each algorithm was instantiated through 5 different initial populations. For each instance, the prediction quality was evaluated using a posture not used in the optimization process.

Results showed that using a genetic algorithm with a representation of polytopes in 26 vertices using the postures set $S_3$ leads to the smallest prediction errors. In perspective of this work, a possible option could be to take into account the vertices neighbourhood using a graph representation. The method should also be extended to other derivative-free algorithms, such as CMA-ES, Bayesian and optimistic optimizations.

## References

[1] A. Skuric, V. Padois, N. Rezzoug and D. Daney, On-Line Feasible Wrench Polytope Evaluation Based on Human Musculoskeletal Models: An Iterative Convex Hull Method, *IEEE Robotics and Automation Letters*, 7(2), 5206-5213, 2022.

[2] B. Grünbaum, *Convex Polytopes*, Springer Science and Business Media, 2013.

[3] Y.-Q. Hu, H. Qian, and Y. Yu, Sequential classification-based optimization for direct policy search, *proceedings of the 31st AAAI Conference on Artificial Intelligence*, pages 2029-2035, San Francisco, CA, 2017

[4] V. Hernandez, N. Rezzoug, P. Gorce, and G. Venture, Force feasible set prediction with artificial neural network and musculoskeletal model, *Computer Methods in Biomechanics and Biomedical Engineering*, 21(14), 740-749, 2018.

[5] K. Holzbaur, W. Murray, S. Delp, A Model of the Upper Extremity for Simulating Musculoskeletal Surgery and Analyzing Neuromuscular Control, *Annals of Biomedical Engineering*, 33(6), 829-840, 2005.