I. Une approche physique▲
Le problème considéré dans l’intégralité de cet article est d’empiler des billes dans un récipient. Le but est de savoir comment la position de chaque bille évolue au cours du temps. Pour cela, il faut connaître les forces qui s’appliquent entre les billes lorsque celles-ci sont en contact. Et chaque force au niveau d’un contact dépend des forces qui s’appliquent aux autres contacts. On utilise alors des méthodes numériques pour résoudre ce problème.
|
|
Pour simplifier le problème, on considère des billes complètement lisses, sans 
frottementWikipedia ni glissement entre billes, ce qui revient à faire abstraction de l’orientation et de la vitesse de rotation de la bille. En effet, les forces de contact sans frottement sur une bille sont perpendiculaires à sa surface, la force passe donc par le centre de la bille et ne génère pas de coupleWikipedia sur celle-ci.
Ainsi, lorsque deux billes se rencontrent, une impulsionWikipedia qui s’applique au point de contact entre les deux billes va annuler la vitesse relative entre elles (on suppose que les billes ne rebondissent pas). Par le principe de conservation de quantité de mouvementWikipedia, cette impulsion va s’appliquer avec la même intensité mais de façon opposée sur les deux billes.
Ainsi, en notant
et
les vitesses du centre de la bille 1 et du centre de la bille 2 et 
un vecteur unitaire indiquant la direction de contact entre les deux billes, nous devons résoudre l’annulation des vitesses après l’impulsion dans la direction de contact :
sachant que les vitesses après impulsion sont
et 
.
Alors, on trouve comme impulsion 
.
Nous appliquons donc cette impulsion si les billes se rapprochent et si elles sont en contact.
|
|
II. Impulsion séquentielle naïve▲
Nous avons trouvé les forces à appliquer lors d’un contact entre billes. Et si, pour résoudre notre système, nous parcourions l’ensemble des contraintes et appliquions l’impulsion de contact ci-dessus ? Nous résoudrions peut-être notre problème. En effet, nous pouvons voir déjà apparaître des propriétés intéressantes : par exemple, l’énergie décroît strictement lors d’une impulsion, cela nous assure que notre méthode de résolution ne va pas diverger.
Ainsi, une première simulation peut être réalisée. À chaque itération, l’accélération de la gravité est ajoutée à la vitesse pour toutes les billes 
.
Ensuite, on parcourt toutes les collisions possibles et, s’il y a une collision, on corrige les vitesses selon les équations déterminées au-dessus. Ce processus peut être itéré plusieurs fois dans un même pas de temps pour améliorer les résultats.
Ensuite, on incrémente la position en fonction des vitesses de toutes les billes
, et on reboucle l’algorithme.
Les résultats sont intéressants, mais les billes commencent à passer les unes à travers les autres. On peut déjà multiplier les passes de calcul d’impulsion. Passer à 20 itérations améliore grandement ce problème mais les billes finissent toujours par se traverser. Pour éviter cela, une force de rappel basée sur la position est nécessaire.
Impulsion naïve sans correction en position
II-A. Stabilisation en position▲
Pour corriger les erreurs en position qui s’accumulent, on peut ajouter une impulsion basée sur l’erreur de position. La façon simple et optimale est de réaliser l’impulsion qui annule l’erreur de position en un pas de temps.
Si l’erreur de position est
et que le pas de temps est
, la correction en vitesse à donner est donc
. Finalement on obtient l’impulsion :
.
En rajoutant la correction en position, on perd la propriété de décroissance stricte de l’énergie lors d’un choc.
On peut donc réaliser une simulation. On observe que les billes ne se passent plus à travers, mais le système ne se stabilise pas. Les billes continuent à trembler.
Impulsion naïve sans anticipation du choc
L’explication vient du fait que la force de rappel en position n’est pas physiquement réaliste et n’est qu’une manière de stabiliser le système ; mais notre méthode est incomplète. En effet, on attend que les billes soient en intersection pour appliquer l’impulsion, mais c’est déjà trop tard, la contrainte en position est violée, les billes se passent à travers et, à la prochaine itération, les billes vont se repousser trop fort.
La solution est d’anticiper la collision et d’appliquer l’impulsion avant que la contrainte soit violée.
II-B. Anticipation du choc▲
L’anticipation du choc est réalisée de manière très simple. Au lieu d’appliquer l’impulsion si les objets sont en intersection, l’impulsion est générée dès que l’impulsion totale (en prenant en compte la correction en position) est positive :
. On garde donc l’équation pour l’impulsion précédente.
L’impulsion se comporte donc différemment lorsque les boules ne se touchent pas encore. On peut distinguer deux cas :
- les boules se rapprochent (
), mais elles sont trop éloignées pour que l’impulsion totale soit positive (
). Cela signifie soit que l’erreur ne sera pas violée à la prochaine itération, soit que les boules n’entreront pas en collision au prochain pas de temps. Aucune impulsion n’est appliquée dans ce cas-là ; - les boules ne se touchent pas (
), mais elles se rapprochent assez pour que l’impulsion totale soit positive
. Cela signifie que l’erreur sera violée à la prochaine itération si aucune impulsion n’est appliquée. L’impulsion est positive, elle est donc appliquée. En analysant le phénomène, on remarque que, après impulsion, au prochain pas de temps, l’erreur de position sera nulle et il ne restera plus qu’une erreur de vitesse à corriger. Le système se stabilise donc très vite.
Si on résume le fonctionnement de l’algorithme lors d’un choc :
- si à l’itération suivante, l’erreur n’est pas violée, on n’applique pas d’impulsion ;
- si à l’itération suivante, l’erreur est violée, on applique une impulsion qui donne une erreur en position nulle à l’itération suivante ;
- passé l’étape 2, l’erreur en position est très faible et il ne reste plus qu’à corriger l’erreur en vitesse.
Toutes ces étapes sont automatiquement réalisées par l’équation totale de l’impulsion ramenée à 0 lorsque celle-ci est négative.
La simulation fonctionne effectivement très bien avec cette méthode !
Simulations réalisées avec stabilisation et anticipation.
Mais nous avons résolu le problème sans faire appel aux formalismes mathématiques usuels tels que la formulation de complémentarité linéaire. Est-ce grave ?
Eh bien on peut modéliser un système physique de différentes manières, ainsi la méthode des impulsions naïves en est une, le problème de complémentarité linéaire en est une autre…
III. Complémentarité linéaire▲
Le problème de complémentarité linéaire est la façon usuelle de décrire les problèmes de contact. De façon simple, le problème est de trouver les impulsions
positives ou nulles telles que l’erreur de vitesse à chaque contact soit :
- positive lorsque la force est nulle ;
- nulle lorsque la force est positive.
Voilà un problème de complémentarité linéaire.
On le pose de façon équivalente :
;
;
.
est le vecteur d’impulsion ;
est l’erreur de vitesse initiale ;
est la matrice qui, multipliée aux vecteurs des impulsions, calcule les erreurs en vitesses générées. est symétrique ;
est l’erreur de vitesse après application des impulsions.
On emploie en général des méthodes itératives pour résoudre le problème. On part d’un vecteur d’impulsion nul et, pendant toutes les itérations, on n’utilise que des impulsions positives. Jusqu’à convergence, les contraintes et 
ne sont pas respectées. Si, en augmentant le nombre d’itérations, les contraintes, et tendent vers 0, alors la méthode converge vers la solution du problème de complémentarité linéaire.
On peut se demander si la méthode des impulsions naïves, qui donne des résultats très intéressants, résout le problème de complémentarité linéaire. La méthode des impulsions naïves ne converge pas vers la solution du problème de complémentarité linéaire. Si les contraintes de positivité sont bien remplies (et), la contrainte de complémentarité ne l’est pas.
IV. Méthode de Gauss-Seidel projetée▲
Pour rendre notre méthode d’impulsion naïve capable de résoudre le problème de complémentarité linéaire, une seule petite modification est nécessaire. Cependant, il faut comprendre d’où vient le problème, pourquoi la méthode ne permet pas de respecter les contraintes de complémentarité.
On a vu que, lorsqu’on appliquait une impulsion pour corriger une erreur entre deux billes, l’erreur devient nulle : la contrainte de complémentarité est donc respectée localement, 
.
Mais il suffit qu’une autre collision se produise en même temps pour que ce ne soit pas le cas : par exemple, l’erreur de vitesse deviendrait positive
, donc
. Pour résoudre ce problème, il faut alors réduire à nouveau
de telle sorte qu‘une contrainte parmi
et
soit respectée et sans que ni l’un ni l’autre terme ne soit négatif. En résumé, on a donc :
- l’impulsion naïve :
; - la méthode de Gauss-Seidel projetée :
.
Cette modification autorise l’impulsion à diminuer lorsque l’erreur de vitesse est par la suite corrigée par un autre choc. Ici,
est la somme des impulsions qui ont été appliquées sur le contact entre les billes 1 et 2 au cours des différentes itérations.
On retombe d’ailleurs sur le solveur utilisé par BulletpyBullet. C’est une façon optimisée de faire une méthode de Gauss-Seidel projetée, car il n’est pas nécessaire de construire la matrice.
Bullet appelle un tel solveur « sequential impulse solver », car son implémentation ressemble beaucoup à la méthode des impulsions successives, mais il résout bien le problème de complémentarité linéaire et est strictement équivalent à la méthode de Gauss-Seidel projetée.
Finalement, à l’œil nu, on ne voit pas beaucoup de différences avec notre précédent solveur, mais, en calculant les contraintes, on voit bien que notre nouvel algorithme converge vers la solution du problème de complémentarité.
Quelques variantes de cette méthode existent aussi comme la méthode de Gauss-Jacobi projetée ou la méthode de sur-relaxation successive projetée.
IV-A. Méthode de sur-relaxation successive projetée▲
Puisque la méthode de Gauss-Seidel permet de corriger des impulsions générées en trop, pourquoi ne pas générer volontairement des impulsions trop fortes ? Par exemple, appliquer une impulsion 1,75 fois plus forte pourra améliorer la vitesse de convergence. Cette méthode est simplement une méthode d’overshooting qui permet de converger plus vite dans certains cas. En effet, pour tout coefficient d’overshooting
, la méthode converge. En général, un coefficient inférieur à 1 diminue la vitesse de convergence et un coefficient trop proche de 2 la diminuera aussi. Les meilleures performances sont atteintes entre 1 et 2. Dans la suite, nous utiliserons un coefficient de sur-relaxation de
.
L’impulsion devient donc :
Il ne faut pas oublier de diviser par
l’impulsion totale pour que l’impulsion finale soit toujours positive.
On peut d’ailleurs régler le paramètre de sur-relaxation dans les moteurs physiques usuels tels que Bullet. Par défaut, celui-ci est à 1.
IV-B. Méthode de Jacobi projetée▲
Pour réaliser la méthode de Jacobi, les impulsions sont d’abord toutes calculées d’un coup. Ensuite, elles sont appliquées aux vitesses, contrairement à la méthode de Gauss-Seidel ou de sur-relaxation, où les vitesses sont calculées à chaque impulsion. Cette méthode converge moins vite que la méthode de Gauss-Seidel, elle nécessite de réduire l’impulsion, sinon le système diverge — un coefficient de
sur l’impulsion évitera que la méthode diverge dans notre cas. Cette méthode est présentée à titre de comparaison, mais son utilisation n’est pas conseillée.
Cette méthode était utilisée dans le moteur physique TokamakWikipedia.
V. Gradient conjugué projeté▲
La méthode de Gauss-Seidel projetée et ses dérivées sont parfois trop lentes à trouver la solution. Pour des problèmes avec de gros ratios de masses, par exemple lorsqu’une bille de 1 kg est empilée sur une bille de 1 g, la méthode de Gauss-Seidel ne réussit pas à résoudre le système en un temps raisonnable. Ce problème se traduit mathématiquement par un mauvais conditionnementWikipedia de la matrice.
Pour résoudre de tels cas de figure, des méthodes avec un plus grand taux de convergence peuvent être utilisées. La méthode du gradient conjugué est proposée. Elle possède la séduisante propriété de résoudre le système au bout d’un nombre fini d’itération.
V-A. Gradient conjugué▲
La méthode du gradient conjugué est un algorithme pour résoudre des systèmes d’équations linéaires dont la matrice est symétrique et définie positiveWikipedia. C’est une méthode qui converge en un nombre fini d’itérations (au plus égal à la dimension du système linéaire). Le gradient conjugué est une méthode de choix pour résoudre les systèmes linéaires en mécanique, particulièrement avec la méthode des éléments finis où la matriceest creuseWikipedia et le temps de calcul diminue beaucoup par rapport au pivot de Gauss.
Voici l’algorithme :
répéter :
tant que
est le résidu, soit l’erreur que l’on cherche à réduire à 0. Il est aussi égal à
. Dans notre cas, c’est l’erreur de vitesse.
est la direction de recherche de la solution.
est la solution qu’on cherche. Dans notre cas, c’est le vecteur des impulsions.
L’un des gros intérêts du gradient conjugué est qu’il ne nécessite pas encore une fois de construire la matrice A. Il suffit d’avoir une fonction qui multiplie un vecteur par A.
V-B. Calcul du produit matriciel▲
Dans notre cas, la fonction « multiplier par A » est calculée en deux étapes. La première étape prend en paramètre un vecteur d’impulsion et applique ces impulsions au niveau des points de contact entre billes et calcule, à partir de ces impulsions, les vitesses générées sur les billes. La deuxième étape récupère les vitesses au niveau des billes et va calculer les erreurs de vitesse au niveau des billes.
|
|
Pour vérifier que notre fonction « multiplier par A » est correcte, on peut construire la matrice A en multipliant successivement les vecteurs de la base canonique par A, soit (1,0,0,0...) puis (0,1,0,0…) puis (0,0,1,0…).
Si la matrice formée des vecteurs ainsi créés est symétrique, c’est que nous sommes sur la bonne voie...
Un premier test peut être réalisé, dans un cas où le problème de contact se réduit à un problème statique non hyperstatique : par exemple, des billes empilées au repos. Le gradient conjugué résout en 5 itérations seulement le problème d’empilement de 5 billes. Cela est prévisible, vu qu’il n’y a que 5 impulsions à trouver.
V-C. Projection▲
Le gradient conjugué fonctionne très bien pour des problèmes linéaires, mais, pour les problèmes de complémentarité linéaire, une modification de l’algorithme est nécessaire. Il y a plusieurs approches possibles dans la littérature. Nous allons nous concentrer sur une formulation qui reste proche du gradient conjugué qu’on peut trouver dans 
Un algorithme de gradient conjugué projeté préconditionné pour la résolution de problèmes unilatéraux, Nicolas Tardieu, Fabien Youbissi, Éric Chamberland
En effet, dans cette formulation, si on applique notre algorithme à un problème linéaire symétrique, le gradient conjugué projeté se comportera exactement comme un gradient conjugué.
Dans la version projetée de l’algorithme du gradient conjugué, il y a quelques modifications à faire :
- lorsqu’un nouveau vecteur d’impulsion
est calculé, il ne faut prendre que les valeurs positives et mettre les valeurs négatives à 0. S‘il y avait des valeurs négatives, c’est que l’état de contrainte a changé ; - le résidu ne peut plus se calculer de la même manière qu’avant, il est nécessaire d’utiliser le calcul direct du résidu :
. Ensuite, il faut projeter le résidu : si l’erreur en vitesse au contact est positive et que l’impulsion est nulle, c’est que la contrainte est respectée, on pose alors le résidu à 0 ; - enfin, la conjugaison n’a d’intérêt que si l’état de contrainte n’a pas changé. Si c’est le cas, la direction de recherche
est réinitialisée au résidu courant
.
Voici l’algorithme :
si
et
alors on pose
répéter :
si
et 
alors on pose 
.
si l’état de contrainte est le même, sinon.
tant que
.
Finalement, l’algorithme est approximativement deux fois plus lent que le gradient conjugué simple, car il nécessite une multiplication par A en plus pour calculer le résidu.
V-D. Pré-conditionnement de Jacobi▲
Des méthodes pour accélérer les performances de convergences du gradient conjugué existent. Il s’agit du pré-conditionnement. Ici un pré-conditionnement de Jacobi a été utilisé, c’est le plus simple et le plus basique. L’amélioration des performances du pré-conditionnement est minime mais non négligeable dans certains cas. Le pré-conditionnement de Jacobi est plutôt simple. Le lecteur peut se référer à La méthode du gradient conjugué pré-conditionné
VI. Benchmark▲
Il est à noter que les méthodes ont quasiment le même temps de calcul par itération, sauf pour les méthodes de gradient conjugué projeté, où il y a deux multiplications matricielles au lieu d’une ce qui la rend deux fois plus lente que les autres.
VI-A. Problème 3D▲
|
|
|
Dans cette simulation, 350 billes de 20 cm de diamètre sont empilées dans une sphère de 2 m de diamètre.
On retrouve bien le classement des performances obtenues. La méthode du gradient conjugué arrive à la précision limite (
) en 1300 itérations. La méthode de sur-relaxation (PSOR) y arrive en 4000 itérations environ.
VI-B. Problème 2D▲
|
|
|
Dans cette simulation, 50 billes de 20 cm de diamètre sont empilées dans un disque de 2 m de diamètre. Le problème est plan.
Il est difficile d’interpréter ces résultats. En effet, toutes les méthodes stagnent à un résidu de
. On pourrait considérer que le problème est résolu, pourtant, au bout de 500 itérations voire 1000 itérations, le gradient conjugué arrive à descendre à un résidu de. Cependant, par sécurité, il est peut-être plus intéressant de stopper le gradient conjugué à un résidu de.
VI-C. Problème 1D▲
|
|
|
Dans cette simulation, 7 billes de 20 cm de diamètre sont empilées. Il y a donc 7 impulsions à trouver. Le gradient conjugué converge alors en 7 itérations. L’intérêt du gradient conjugué par rapport aux autres méthodes est ici très net. Ici les courbes NAIV et PGS sont superposées.
VI-D. Problème mal conditionné▲
|
|
 |
Dans ce problème, la bille du haut est 1000 fois plus lourde que celle du bas. Le problème est donc mal conditionné. On voit clairement que ce problème pose beaucoup de difficultés aux méthodes basées sur la méthode de Gauss-Seidel, alors que le gradient conjugué converge en 2 itérations seulement (vu qu’il y a 2 degrés de liberté). Ici les courbes NAIV et PGS sont superposées.
VI-E. Conclusion▲
Globalement, dans tous les cas, on peut classer les vitesses de convergence des méthodes comme suit :
Jacobi projeté (PJ) < impulsions Naives (NAIV) < Gauss-Seidel projeté (PGS) < sur-relaxation projetée (PSOR) < gradient conjugué projeté (PCG) < gradient conjugué pré-conditionné projeté (PPCG).
Nous avons vu que les méthodes du gradient conjugué et de sur-relaxation convergent plus vite que la méthode de Gauss-Seidel. Pourtant, la méthode de Gauss-Seidel est la plus utilisée dans la simulation interactive.
En effet, par exemple pour le premier benchmark, la méthode Gauss-Seidel est dépassée seulement à partir de l’itération 100. Mais, pour une simulation interactive, 100 itérations, c’est déjà trop long. Dans notre cas c’est le nombre d’itérations limite de notre solveur pour le premier benchmark avec un fonctionnement à 50hz. Pourtant, dans notre cas, la rotation et les frottements des billes ont été négligés.
Ainsi, typiquement, un solveur interactif effectue seulement une vingtaine d’itérations. Au bout de 20 itérations, la méthode de Gauss-Seidel fournira un meilleur résultat que les deux autres. Voilà pourquoi la méthode de Gauss-Seidel est la plus utilisée dans un contexte interactif.
Ainsi, faire une simulation du premier ou deuxième benchmark avec la méthode du gradient conjugué et 20 itérations seulement donne un résultat ridicule : les billes remuent de manière visible, alors qu’avec les autres méthodes le système se stabilise.
PPCG avec trop peu d’itérations.
Le dernier benchmark permet aussi de voir la limite de la méthode de sur-relaxation en temps réel. Avec 20 itérations, la bille du dessous acquiert beaucoup de vitesse avec les chocs successifs ; au bout de 20 itérations, la bille du dessous aura une telle vitesse qu’elle passera à travers l’autre bille.
VII. Sources▲
Le programme a été réalisé en C++ avec Qt 4.8 pour linux et Qt 5 pour windows ainsi que Eigen3. Python 3.6 a été utilisé pour l’affichage des courbes de performances.
Voici l’architecture du projet :
une classeWindowhéritée deQGLWidgetgère l’affichage OpenGL et les entrées sorties clavier. Pour cela un raycasting sur les billes est réalisé en GLSL (calcul principalement réalisé par le fragment shader qui est fourni). Avec les touches du clavier, il est aussi possible de modifier le scénario, le solveur et le nombre d’itérations du solveur.
Une classeSolveurhéritée deQThreadcalcule en parallèle l’évolution des billes avec le type de solveur sélectionné.
La librairieEigenest utilisée pour les opérations vectorielles simples en double précision (produit scalaire, norme de vecteur…).
Lorsque le programme principal est fermé, il génère avant de quitter des courbes de performances qui dépendent de l’état du système au moment où celui-ci est fermé. Pour cela, la classe Solveur enregistre les courbes de performances dans un fichier numpy (.npy) puis lance le programme affiche.py qui doit se trouver dans le même répertoire que le fichier benchmark.npy créé et qui affiche les courbes de performances.
VIII. Remerciements▲
Merci à Dourouc5 pour m’avoir proposé d’écrire cet article et pour l’avoir relu. Merci aussi à LittleWhite, Escartefigue et Bousk pour leur relecture.