1 La brachistochrone

Contexte: vous êtes ingénieur mandaté pour la conception de nouvelles montagnes russes pour un parc d’attraction. Afin de maximiser la vitesse du train, vous devez concevoir la rampe qui va partir du point le plus haut, arriver au point le plus bas (à une hauteur \(H\) sous le point le plus haut et à une distance horizontale \(L\) du point haut) le plus vite possible, sous l’influence de la gravité uniquement.

Afin de simplifier le raisonnement initial, on considère que le train est une masse ponctuelle qui glisse sans frottements le long de la rampe.

Note

Ce TD est inspiré de Mittal (2008) (lien PDF) qui développe la discrétisation du problème.

On va chercher une solution numérique de ce problème. Pour cela, formalisons:

On cherche une courbe \(y: [0, L] \rightarrow \mathbb R\) avec \(y(0) = 0\) et \(y(L) = H\) telle qu’une particule de masse \(m\) subissant l’accélération de la gravité \(\vec g\) mette le moins de temps pour aller de \(A = (0, 0)\) à \(B = (L, H)\) en suivant la courbe \(y(x)\).

Avertissement

Attention à la convention de signe pour \(y\): l’axe est dans le même sens que \(\vec g\).

Note

Ce problème mathématique est le problème de la brachistochrone, qui a une importance particulière dans l’histoire des mathématiques car c’est un problème d’optimisation où la solution est une fonction et non un nombre ou un vecteur. Il a d’ailleurs fait l’objet d’une compétition entre les frères Bernoulli, Newton et Leibniz, entre autres.

Le problème admet une solution analytique que l’on ne va pas détailler dans cet exercice : nous nous serviront de la solution uniquement pour valider l’approche numérique.

1.1 Fonction à minimiser

La formulation du problème nous donne la fonction à minimiser : le temps de trajet \(t(y)\) qui est fonction de la forme de la rampe choisie. On appelle \(s\) la position de la particule le long de la rampe (\(s \neq x\) !). La distance totale parcourue par la particule est :

\[ D = \int_A^B \mathrm ds, \]

soit la somme de petits morceaux de courbe \(\mathrm ds\). On note \(v(s)\) la vitesse de la particule à la position \(s\). Pour parcourir la distance \(\mathrm ds\), le temps nécessaire est:

\[ \mathrm dt = \frac{\mathrm ds}{v(s)}, \]

ce qui donne le temps de trajet total (que l’on doit minimiser) :

\[ T = \int_A^B \mathrm dt = \int_A^B \frac{1}{v(s)}\,\mathrm ds \]

On peut exprimer la longueur de l’élément de courbe \(\mathrm ds\) avec Pythagore:

\[ \mathrm ds = \sqrt{\mathrm dx^2 + \mathrm dy^2} = \sqrt{1 + \left(\frac{\mathrm dy}{\mathrm dx}\right)^2}\,\mathrm dx = \sqrt{1 + y'(x)^2}\,\mathrm dx\]

On peut donc, par changement de variable, exprimer \(T\) en fonction de \(y(x)\):

\[T(y) = \int_0^L \frac{\sqrt{1 + y'^2}}{v(x)}\,\mathrm dx\]

Consigne

En utilisant la conservation de l’énergie mécanique, montrer que le temps de trajet s’exprime :

\[ \sqrt{2g}\cdot T(y) = \int_0^L \sqrt{\frac{1 + y'^2}{y}}\,\mathrm dx \tag{1.1}\]

Puisque minimiser \(\sqrt{2g}\cdot T\) est équivalent à minimiser \(T\), on pose par la suite :

\[T(y) = \int_0^L \sqrt{\frac{1 + y'^2}{y}}\,\mathrm dx \tag{1.2}\]

Le problème d’optimisation est donc le suivant:

\[\begin{aligned} \min_{y}\quad\ & T(y) = \int_0^L \sqrt{\frac{1 + y'^2}{y}}\,\mathrm dx \\\text{sous contraintes } \quad & y(0) = 0\text{ et }y(L) = H\end{aligned}\]

1.2 Discrétisation

Pour résoudre numériquement le problème, on procède à une discrétisation par éléments finis de la courbe \(y\).

On discrétise l’intervalle \([0, L]\) en \(N_e\) éléments de même taille \(\Delta x\). On note \(x_1,\ldots x_N\) les positions des nœuds, \(N = N_e+1\). On fait une discrétisation de \(y(x)\) avec des fonctions de formes linéaires:

\[ y_h(x) = \sum_{i = 1}^N y_i \phi_i(x). \]

Les valeurs \(y_i\) sont les valeurs de \(y_h\) aux nœuds \(x_i\) et \(\phi_i\) sont les fonctions d’interpolations P1 classiques.

On peut donc discrétiser l’intégrale de \(T\):

\[ T(y_h) = \int_0^L \sqrt{\frac{1 + y_h'^2}{y_h}}\,\mathrm dx = \sum_{i=1}^{N_e} \int_{x_i}^{x_{i+1}} \sqrt{\frac{1 + y_h'^2}{y_h}}\,\mathrm dx = \sum_{i=1}^{N_e} T_i \tag{1.3}\]

Consigne

Monter que

\[ T_i = \int_{x_i}^{x_{i+1}} \sqrt{\frac{1 + y_h'^2}{y}}\,\mathrm dx = \frac{\sqrt{\Delta x^2 + \Delta y_i^2}}{ \sqrt{y_{i+1}} + \sqrt{y_i}} \qquad\text{avec } \Delta y_i = y_{i+1} - y_i \tag{1.4}\]

On pourra se servir du fait que dans l’intervale \(]x_i, x_{i+1}[\) on a \(y_h' = \Delta y_i / \Delta x\) qui est constant.

1.3 Optimisation à deux degrés de liberté (avec contrainte)

On va en premier lieu explorer graphiquement ce qui se passe avec deux éléments. On a donc trois valeurs inconnues \(y_1\), \(y_2\) et \(y_3\). On applique directement la contrainte \(y(0) = 0\) en posant \(y_1 = 0\), mais on va laisser la valeur de \(y_3\) libre et imposer la contrainte avec un multiplicateur de Lagrange.

Consigne

Monter que le temps de trajet en fonction de \(y_2\) et \(y_3\) s’exprime:

\[ T(y_2, y_3) = \frac{\sqrt{(L/2)^2 + y_2^2}}{\sqrt{y_2}} + \frac{\sqrt{(L/2)^2 + (y_3-y_2)^2}}{\sqrt{y_3}+\sqrt{y_2}}\]

On pose \(L\) = 10 et \(H\) = 3. On souhaite afficher sur un graphe les courbes de niveau de \(T(y_2, y_3)\) et la courbe de la contrainte \(y_3 = H\), et chercher graphiquement le point tangent où se trouve l’optimum. Pour cela, il faut implémenter la fonction

import numpy as np
import matplotlib.pyplot as plt

L = 10
H = 3

def temps_2ddl(y2, y3):
    return ...

Et construire l’espace de coordonnées \((y_2, y_3)\in [0, 2H]\times[0, 2H]\):

y2 = np.linspace(0, 2 * H, 100)
y3 = np.linspace(0, 2 * H, 100)
y2, y3 = np.meshgrid(y2, y3)

Consigne

En utilisant la fonction contour() de matplotlib, chercher graphiquement la valeur de \(y_2\) optimale. On pourra affiner les bornes de \(y_2\) et \(y_3\) dans le code ci-dessus (2 chiffres significatifs seront suiffsants).

Solution

Consigne

Afficher la forme de la rampe pour la solution trouvée.

Solution

1.4 Optimisation à \(N\) degrés de liberté

On cherche maintenant à généraliser l’approche d’optimisation à \(N\) degrés de liberté. Pour cela, on va se servir de deux bibliothèques Python:

Scipy est une bibliothèque généraliste de calcul scientifique. On va surtout se servir de la fonction minimize qui permet de minimiser des fonctions sous contraintes.
JAX qui est une bibliothèque similaire à Numpy, mais qui permet de calculer automatiquement le gradient de n’importe quelle fonction (cela nous évite de dériver l’expression du gradient \(\partial T_i/\partial y_j\) (Équation 1.4) à la main).

import jax.numpy as jnp
import scipy.optimize as opt

from jax import grad

Consigne

À l’aide des fonctions jnp.sqrt, jnp.sum, implémenter l’expression de l’Équation 1.3 dans la fonction suivante.

def temps_nddl(y, x):
    ...

Note

Pour faciliter l’implémentation, on note que les expressions suivantes sont équivalentes:

y[:] \(\Leftrightarrow \{y_1, y_2, \ldots, y_{N}\}\)
y[1:] \(\Leftrightarrow \{y_2, y_3, \ldots, y_{N}\}\)
y[:-1] \(\Leftrightarrow \{y_1, y_2, \ldots, y_{N-1}\}\)

Essayer d’écrire \(\Delta y_i = y_{i+1} - y_i\) à l’aide des expressions ci-dessus (il faut éviter d’écrire une boucle sur \(i\)).

Avertissement

Pendant la procédure d’optimisation il est possible que certaines valuers de \(y_i\) soient négatives. Pour éviter les soucis avec \(\sqrt{y_i}\), il faut calculer \(\sqrt{|y_i|}\) avec jnp.abs qui retourne la valeur absolue.

Vérifier que la nouvelle implémentation donne le même temps que l’implémentation à 2 DDL.

x_test2 = jnp.array([0, L/2, L])
y_test2 = jnp.array([0, H/2, H])

print(f"Temps NDDL = {temps_nddl(y_test2, x_test2)}\nTemps 2DDL = {temps_2ddl(y_test2[1], y_test2[2])}")
# ou encore
assert temps_nddl(y_test2, x_test2) == temps_2ddl(y_test2[1], y_test2[2])

Temps NDDL = 6.027713775634766
Temps 2DDL = 6.027713775634766

On note \(\mathbf y = [y_1, y_2, \ldots y_N]^\mathsf T \in \mathbb R^{N\times 1}\) le vecteur colonne contenant les inconnes. Le problème d’optimisation discrétisé que l’on cherche à résoudre est le suivant:

\[ \begin{aligned} \min_{\mathbf y}\quad & \sum_{i=1}^{N-1} \frac{\sqrt{\Delta x^2 + \Delta y_i^2}}{ \sqrt{y_{i+1}} + \sqrt{y_i}} \\ \text{sous contraintes} \quad & y_1 = 0 \\ & y_N = H \end{aligned} \]

Afin d’utiliser Scipy pour la résolution de ce problème, il nous faut exprimer les contraintes comme un système d’équations linéaires:

\[ \mathbf A\mathbf y = \mathbf b\]

où \(\mathbf A\in \mathbb R^{2\times N}\) et \(\mathbf b\in \mathbb{R}^{2\times 1}\). Le problème d’optimisation s’écrit donc:

\[ \begin{aligned} \min_{\mathbf y}\quad & \sum_{i=1}^{N-1} \frac{\sqrt{\Delta x^2 + \Delta y_i^2}}{ \sqrt{y_{i+1}} + \sqrt{y_i}} \\ \text{sous contraintes} \quad & \mathbf A\mathbf y = \mathbf b \end{aligned} \]

Consigne

Trouver les matrices \(\mathbf A\) et \(\mathbf b\) qui traduisent les deux contraintes \(y_1 = 0\) et \(y_2 = H\) et les implémenter.

N = 6 # par exemple

A = np.zeros((2, N))
b = np.zeros(2)

# Remplir les bonnes valeurs de A et b ici
# ...

Pour intégrer les contraintes d’égalité, on crée un objet spécial de type opt.LinearConstraint:

bords = opt.LinearConstraint(A, b, b)

Puisque l’expression du temps contient des racines carrées, on ajoute une contrainte de positivité pour \(\mathbf y\):

contrainte_positive = opt.LinearConstraint(np.eye(N), 0, np.inf)

On doit également créer le tableau des valeurs discrètes de \(x\) et les valeurs initiales de \(\mathbf y\).

x = jnp.linspace(0, L, N)
y_initial = jnp.ones_like(x)
temps_nddl(y_initial, x)

Array(5., dtype=float32)

On peut enfin appeler minimize:

solution = opt.minimize(temps_nddl, y_initial, args=(x,), jac=grad(temps_nddl), constraints=(bords, contrainte_positive), method='trust-constr')
print(f"Optimisation réussie ? {solution.success}")

Optimisation réussie ? True

Consigne

Le tableau solution.x contient les valeurs optimales de \(\mathbf y\). Afficher la solution obtenue pour \(N = 6\).

Solution

1.5 Convergence

On cherche à caractériser la solution quand \(N\) augmente.

Consigne

Utiliser le code précédent pour implémenter la fonction solution_optimale(N) qui donne la solution pour un nombre de segments donné.

def solution_optimale(N):
    x = jnp.linspace(0, L, N)
    
    ...
    
    return x, solution.x

Choisissez 5 valeurs de \(N\) entre 2 et 50 pour et afficher chaque solution.

Solution

1.6 Contrainte d’obstacle

La fonction opt.LinearConstraint(A, l, u) permet d’intégrer des contraintes de la forme :

\[ \mathbf {l} \leq \mathbf A \mathbf y \leq \mathbf u \]

On s’est servi ci-dessus du fait que quand l == u on a une contrainte d’égalité \(\mathbf {Ay} = \mathbf l = \mathbf u\). On a aussi défini la contrainte \(\mathbf y \geq 0\) en posant \(\mathbf A = \mathbf I_{N\times N}\), \(\mathbf l = 0\) et \(\mathbf u = +\infty\).

Consigne

On imagine que le rail de la montagne russe ne doit pas dépasser un obstactle défini par la fonction \(h(x)\), c’est-à-dire que le problème de minimisation continu est:

\[\begin{aligned} \min_{y}\quad\ & T(y) = \int_0^L \sqrt{\frac{1 + y'^2}{y}}\,\mathrm dx \\\text{sous contraintes } \quad & y(0) = 0\text{ et }y(L) = H \\ & y(x) \leq h(x) \end{aligned}\]

À l’aide de opt.LinearConstraint, implémenter la contrainte d’obstacle pour \(h(x) = 2H - x\left(1 - \frac{x}{L}\right)\). Tester avec \(N=20\) et afficher la solution optimale.

Solution

Références

Mittal, Aditya. 2008. « Numerical Solution to the Brachistochrone Problem ». Stanford University.