# Solveurs de systèmes de complémentarité et application à la tenségrité

Pierre Alart, Frédéric Dubois, David Dureisseix, Shéhérazade Nineb

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

## To cite this version:

Pierre Alart, Frédéric Dubois, David Dureisseix, Shéhérazade Nineb. Solveurs de systèmes de complémentarité et application à la tenségrité. 17e Congrès Français de Mécanique - CFM 2005, Sep 2005, Troyes, France. pp.1-6. hal-00325950

HAL Id: hal-00325950

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

Submitted on 21 Jan 2022

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Solveurs de systèmes de complémentarité et application à la tenségrité

P. Alart, F. Dubois, D. Dureisseix, S. Nineb

Laboratoire de Mécanique et Génie Civil, UM2-CNRS, CC048, Place Eugène Bataillon F-34095 Montpellier
CEDEX 5
nineb@lmgc.univ-montp2.fr

## Résumé :

Un Problème de Complémentarité Linéaire (Linear Complementary Problem-LCP) relie linéairement deux variables positives et orthogonales entre elles. On propose de comparer des méthodes de résolution de LCP, appliquées à un système de tenségrité, système réticulé à barres et câbles (ces derniers devant rester en tension) dont la rigidité et la stabilité sont conditionnées par l'existence d'un état d'autocontrainte.

## Abstract :

A Linear Complementary Problem (LCP) deals with two dual variables that have to be positive and orthogonal, and that are linearly dependent on each other. Solving algorithms for LCP problems are compared when applied to a tensegrity system. Such a system is built with bars and cables (that should remain in tension). Its stiffness and stability is subjected to a self balanced internal stress field.

## Mots-clefs :

## LCP, non régularité, système discret

### 1 Motivations

De nombreux problèmes non réguliers peuvent se ramener après discrétisation, à un système algébrique contenant des contraintes d'inégalités. Le problème générique est un problème de complémentarité linéaire (LCP-Linear Complementary Problem). Celui-ci est constitué d'un système linéaire et de contraintes de positivité et d'orthogonalité sur un couple de variables duales. La forme canonique d'un tel problème est la suivante :

$$ \begin{cases} w - Mz = q \\ 0 \le w \perp z \ge 0 \end{cases} \tag{1} $$

Un problème typique en mécanique est le cas du contact unilatéral. Il existe plusieurs extensions du LCP, dues à différentes problématiques ; par exemple le couplage contact/frottement induit une généralisation du LCP en NLCP (Non Linear Complementary Problem). Pour traiter ces problèmes, il n'existe que peu de librairies et d'outils génériques par comparaison aux résolutions de systèmes linéaires. Pourtant les problèmes de complémentarité intéressent plusieurs communautés, au delà de la mécanique, au point qu'ils sont au cœur d'un projet européen, nommé SICONOS (SImulation and COntrol of NOnSmooth dynamical Systems), regroupant mathématiciens, mécaniciens, roboticiens, électrotechniciens, automaticiens et économistes. Une de nos contributions à ce projet consiste en l'élaboration d'une librairie de solveurs dans laquelle on trouve des algorithmes dédiés à la résolution de LCP entre autres. Les méthodes de résolutions numériques que nous proposons de comparer sont les suivantes : la méthode de Lemke, méthode directe de pivotage [1], la méthode de Gauss-Seidel Non linéaire (GSNL), méthode itérative non linéaire utilisée pour les simulations des matériaux granulaires [2], une méthode de type Gradient Projeté Conjugué (GPC), algorithme que l'on peut qualifier de quasi-optimisation [3] et la méthode LA-TIN, stratégie itérative non incrémentale qui a été initialement proposée pour résoudre les problèmes non linéaires d'évolution [4]. Ces dernières ont été simplifiées pour le traitement des problèmes de
complémentarité. La comparaison de ces différents algorithmes se fera sur un système réticulé spatial particulier que l'on nomme système de tenségrité.

## 2 Formulation LCP des problèmes de tenségrité

### 2.1 Système de tenségrité

Les systèmes de tenségrité [5] peuvent-être rencontrés dans le domaine architectural, artistique voire même plus récemment en biomécanique cellulaire. Ce sont des systèmes réticulés spatiaux qui possèdent des éléments rectilignes soumis à des sollicitations de compression ou de traction. Les liaisons entre les éléments se font par l'intermédiaire de rotules parfaites nommées nœuds. Les éléments rigides en compression portent le nom de barres ceux rigides en traction portent le nom de câbles.

On s'intéresse à la mise en place d'une formulation de type LCP pour ces systèmes. Pour ce faire, décrivons tout d'abord le problème associé.

La configuration de référence choisie est celle donnée par l'état du système avant assemblage. On suppose que le système réticulé est constitué de *m* éléments (barres et câbles) liés par *n* nœuds. Le nombre de degrés de libertés imposés aux nœuds du système par l'extérieur (c'est à dire par le bâti) est noté *k*. On en déduit ainsi le nombre *N* de ddl du système dans l'espace, soit $N = 3n - k$.

Les *N* équations d'équilibre du système de tenségrité, qui traduisent l'équilibre des nœuds sous l'action des forces extérieures et des efforts normaux, conduisent à un système d'équations linéaires pouvant se mettre sous la forme matricielle suivante :

$$A q = f \qquad (2)$$

où *A* est la matrice d'équilibre à *N* lignes et *m* colonnes, *q* est le vecteur à *m* composantes des efforts dans les éléments et *f* est le vecteur à *N* composantes des forces extérieures appliquées aux nœuds. En séparant les valeurs relatives aux câbles (indicées par *c*) et celles relatives aux barres (indicées par *b*), on a $A = [A_b, A_c]$ et l'équilibre devient $A_b q_b + A_c q_c = f$.

La condition de compatibilité qui lie le vecteur à *N* composantes des déplacements des nœuds *u* depuis l'état de référence au vecteur à *m* composantes de variations de longueur *e* des éléments s'écrit, en notant la transposée avec un exposant *t*,

$$e = A^t u \quad \text{ou} \quad e_b = A_b^t u \quad \text{et} \quad e_c = A_c^t u \qquad (3)$$

La loi de comportement linéaire vérifiée par les barres s'écrit :

$$q_b = k_b e_b \qquad (4)$$

où $k_b$ est une matrice diagonale constituée des raideurs des barres.

La loi de comportement des câbles est non régulière : l'effort interne est non nul sous traction et nul sous compression. Après changement de variable, on retrouve le graphe de Signorini avec pour variables $q_c$ et $R_c$. L'expression de la variable introduite par le changement de variable est alors donnée par :

$$R_c = k_c^{-1} q_c - e_c - e_c^0. \qquad (5)$$

où $e_c^0$ est la variation de longueur des câbles après assemblage (précharge).

On cherche dans la suite à formuler le problème précédent sous la forme canonique d'un système LCP.

### 2.2 Formulation de type décomposition de domaine

La non régularité étant essentiellement contenue dans le comportement des câbles, il paraît naturel d'utiliser une approche du type de celles adoptées en décomposition de domaine, après avoir séparé
les éléments constitutifs d'un système de tenségrité en un sous domaine contenant les câbles, et en un sous domaine contenant les barres. La loi de comportement des barres (4) associée à l'équilibre mène à l'expression suivante $K_b u + A_c q_c = f$ où $K_b = A_b k_b A_b^t$ est la matrice de rigidité du système de barres. Celle-ci est potentiellement singulière (mais symétrique et positive), tout comme celle du système de câbles, $K_c = A_c k_c A_c^t$, la somme des deux $K = K_b + K_c$ étant supposée régulière si le problème est initialement bien posé. La condensation des inconnues $u$ peut être réalisée en utilisant une inverse généralisée $K_b^+$ de $K_b$ [6] : en notant $R$ la matrice constituée des modes singuliers de $K_b$ ($K_b R = 0$) et $\alpha$ les coordonnées possibles du déplacement sur ces modes. L'emploi de l'inverse généralisée nous permet alors d'obtenir le système (6) :

$$ \left\{ \begin{array}{l} u = K_b^+(f - A_c q_c) + R\alpha \\ R^t(f - A_c q_c) = 0 \end{array} \right. \qquad (6) $$

La deuxième équation traduit la nécessaire orthogonalité aux modes singuliers. Une combinaison adéquate de (6) et (5) conduit au système qui suit :

$$ \left\{ \begin{array}{l} W q_c - N R_c = g \\ 0 \le q_c \perp R_c \ge 0 \end{array} \right. \qquad (7) $$

$$ \text{où } N = k_c A_c^t R M^{-1} R^t A_c k_c - k_c \text{ avec } M = R^t A_c k_c A_c^t R, \\ W = 1 + N A_c^t K_b^+ A_c $$

$$ \text{et } g = k_c A_c^t K_b^+ f + k_c A_c^t R M^{-1} R^t (f - A_c k_c e_c^0 - K_c K_b^+ f) + k_c e_c^0. $$

La forme trouvée n'est pas celle d'un LCP standard. En outre W n'est pas symétrique et ceci restreint le champ des méthodes de résolution applicables, notamment exclut celles définies sur un problème d'optimisation (GPC par exemple). De plus, la résolution de ce système peut s'avérer coûteuse car les modes singuliers sont nombreux.

Une formulation plus adaptée est nécessaire. Elle peut cependant être encore construite à partir de la séparation initiale en barres et câbles.

## 2.3 Formulation algébrique

En utilisant simultanément l'équilibre, la loi de comportement des barres et le changement de variable proposé pour les câbles, on arrive à l'expression suivante :

$$ Ku + A_c k_c R_c = f - A_c k_c e_c^0 \qquad (8) $$

Une formulation symétrique pour $R_c$ peut être obtenue en multipliant à gauche par $k_c A_c^t K^{-1}$, cela donne :

$$ k_c A_c^t u + k_c A_c^t K^{-1} A_c k_c R_c = k_c A_c^t K^{-1} (f - A_c k_c e_c^0) \qquad (9) $$

$e_c = A_c^t u$ peut finalement être remplacé en réutilisant le changement de variable (5), et on obtient :

$$ q_c - M R_c = h \qquad (10) $$

$$ \text{où } M = k_c - k_c A_c^t K^{-1} A_c k_c \text{ et } h = k_c A_c^t K^{-1} (f - A_c k_c e_c^0) + k_c e_c^0. $$

Finalement, nous parvenons à la forme souhaitée d'un LCP standard :

$$ \left\{ \begin{array}{l} q_c - M R_c = h \\ 0 \le q_c \perp R_c \ge 0 \end{array} \right. \qquad (11) $$

On peut montrer qu'outre la symétrie, M possède la propriété de positivité et que la dimension de son noyau est inférieure ou égale à celle du noyau de $K_b$ (hors du champ d'investigation de cet article). Cependant, M est généralement singulière, ce qui complique l'analyse mathématique du problème et la maîtrise du comportement des algorithmes de résolution.
## 2.4 Illustrations

Nous décrivons dans cette partie brièvement deux cas tests bidimensionnels, qui bien qu'étant des assemblages de barres et de câbles, ne sont pas à proprement parler des systèmes de tenségrité : celui d'un "mât" et celui d'un "cerf-volant".

<figure><img src="image_3.png" /><figcaption>FIG. 1 – Mât à gauche et cerf-volant à droite</figcaption></figure>

- Le mât. Nous présentons le cas non symétrique où l'angle $α_1$ est différent de l'angle $α_2$. Les câbles sont supposés avoir la même raideur $k_c$. L'écriture de l'équilibre dans la base $(n, t)$, liée à la barre, donne les expressions suivantes :

$$A_c = \begin{pmatrix} -\cos \alpha_1 & -\cos \alpha_2 \\ \sin \alpha_1 & -\sin \alpha_2 \end{pmatrix}, A_b = \begin{pmatrix} -1 \\ 0 \end{pmatrix}$$

et

$$K = k_b \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} + k_c \begin{pmatrix} \cos^2 \alpha_1 + \cos^2 \alpha_2 & \frac{1}{2}[-\sin(2\alpha_1) + \sin(2\alpha_2)] \\ \frac{1}{2}[-\sin(2\alpha_1) + \sin(2\alpha_2)] & (\sin^2 \alpha_1 + \sin^2 \alpha_2) \end{pmatrix}$$

Dans ce cas, la taille du noyau de $K_b$, comme celle de $M$, est de 1. Comme $M$ est copositive plus, l'algorithme de Lemke est applicable comme méthode directe. On peut vérifier que les différentes

<figure><img src="image_4.png" /><figcaption>FIG. 2 – Effort interne en fonction de la charge imposée verticalement.</figcaption></figure>

algorithmes redonnent la même solution tant qu'il y a unicité (c'est à dire non détente simultanée des câbles) et qu'elles ne diffèrent que de R quand il y a plusieurs solutions différentes (détente simultanée des câbles au delà d'une charge critique).

- Le cerf-volant. Pour traiter ce cas, il faut écrire l'équilibre des 4 nœuds et écarter ensuite les différentes variables sur lesquelles les déplacements ont été imposés. A ce stade, on obtient le système matriciel (2). Dans ce cas, la dimension du noyau de $M$ est 3.
# 3 Exemple d'une nappe de tenségrité

## 3.1 Description du problème

Les structures étudiées sont des grilles planes [7], obtenues par assemblage de 9 puis de 256 demi-cuboctaèdre (cf. (3)) (quadruplex particulier ayant 8 nœuds, 12 câbles et 4 barres). L'assemblage de ces modules se fait au niveau de leurs nœuds. Connaissant un état d'autocontrainte du module c'est à dire un état de tension non nul vérifiant l'équilibre sans chargement extérieur, on peut déduire un état d'autocontrainte de l'ensemble. Un déplacement nul est imposé sur les bords gauche et droit de la nappe, et des efforts verticaux constants sont appliqués sur tous les nœuds de celle-ci.

On dénombre respectivement 108 et 3072 contraintes de complémentarité sur ces deux cas. La matrice M étant copositive plus, la méthode de Lemke est applicable.

## 3.2 Comparaison de solveurs LCP

Pour chaque méthode itérative, on choisit comme solution de référence celle obtenue soit par la méthode directe dans le cas de la petite nappe, soit celle obtenue après un grand nombre d'itérations par Gauss-Seidel que l'on notera $R_c^{ref}$. Pour comparer la qualité des solutions, on effectue le calcul de la différence entre les solutions $R_c^i$ obtenues par les différentes méthodes numériques et $R_c^{ref}$, qui représente l'erreur donnée par la formule ci-dessous :

$$ \eta = \sqrt{(R_c^i - R_c^{ref})^t M (R_c^i - R_c^{ref}) / (R_c^i + R_c^{ref})^t M (R_c^i + R_c^{ref})} $$

Les résultats sont donnés par la figure 3. Différents paramètres de direction de recherche, k, ont été pris pour la LATIN.

<figure><img src="image_5.png" /><figcaption>FIG. 3 – Résultats de convergence sur la structure à 9 modules à gauche, sur la structure à 256 modules au centre et déformée à droite.</figcaption></figure>

La méthode de Lemke, montre ses limites en terme de taille de matrice assez rapidement, ne convergeant plus au delà d'une centaine de contraintes. A cette étape de l'étude, il est difficile de séparer ce qui relève du comportement standard d'une méthode directe (qui reste performante dans le cas linéaire pour plusieurs milliers de variables), de ce qui est dû à la singularité de la matrice M (dont la taille du noyau peut être assez conséquent).

La comparaison des méthodes itératives en nombre d'itérations montre une dépendance vis à vis de la taille du problème. Pour un problème de taille réduite, la méthode de Gauss-Seidel converge très rapidement; pour la méthode LATIN, il existe une valeur optimale du paramètre k de la direction de recherche, qui dépend de la taille du problème.
## 4 Conclusions

Pour conclure sur la performance des méthodes, il est encore nécessaire de faire une étude plus précise du coût relatif, en terme de complexité. Une comparaison plus large d’algorithmes va exiger de fiabiliser les méthodes itératives, en particulier vis à vis du test d’arrêt, car il semble que la singularité du problème diminue la robustesse de ces méthodes à l’exception de Gauss-Seidel qui confirme une robustesse déjà constatée sur d’autres types de problèmes.

## Références

[1] Murty Katta. *Principal pivoting methods for LCP*. Department of Industrial and Operations Engineering, University of Michigan, 1997.

[2] M. Jean. The non-smooth contact dynamics method. *Comp. Meth. Appl. Mech. Engrg*, 177:235–257, 1999.

[3] P. Alart and M. Renouf. Un nouvel algorithme de quasi-optimisation pour la résolution des problèmes multicontacts et application aux milieux granulaires. In *Actes du 16e Congrès Français de Mécanique*, pages 1–4, Nice, 2003.

[4] P. Ladevèze. *Nonlinear Computational Structural Mechanics — New Approaches and Non-Incremental Methods of Calculation*. Springer Verlag, 1999.

[5] R. Motro. *Tensegrity*. Hermes Science Publishing, London, 2003.

[6] Charbel Farhat and Daniel Rixen. *Encyclopedia of Vibration*. Academic Press, 2002. ISBN 0-12-227085-1.

[7] Kazi-Aoual M.N. Quirant, J. and Motro R. Designing tensegrity systems : the case of a double layer grid. *Engineering Structures*, 25(9):1121–1130, 2003.