import jax.numpy as jnp
import scipy.optimize as opt
from jax import grad2 La caténaire
caténaire1 .te.nɛʁ féminin
(Géométrie) Forme que prend un câble (ou une chaîne) lorsqu’il est suspendu par ses extrémités et soumis à une force gravitationnelle uniforme (son propre poids).
Les câbles sont des structures mécaniques particulières en ce sens qu’ils sont souvent modélisé comme des éléments capables de résister à des efforts de traction uniquement. Un câble soumis à sont poids propre décrit donc une courbe dont les efforts intérieurs sont uniquement en traction. On appelle cette courbe la caténaire, et on parle de structure funiculaire, qui ne possède aucun moment de flexion sous poids propre.
En renversant la caténaire, on obtient une arche qui ne subit que des efforts de compressions sous son poids propre. Un exemple spectaculaire est la Gateway Arch à St Louis dans le Missouri.
Les structures funiculaires sont particulièrement intéressantes pour optimiser la quantité de matière des structures en béton. Ce matériau a une excellente performance en compression mais est très fragile et peu résistant en traction. Pour construire une poutre en béton, il faut donc renforcer les parties tendues par des armatures en acier, ce qui a des conséquences sur la durée de vie des structures. La conception d’une structure intégralement comprimée permet donc d’éviter les efforts de tractions et de minimiser ou d’éliminer le besoin d’armature, permettant une meilleure durabilité des structures (p. ex. le Panthéon de Rome est en béton non-armé, mais ce n’est pas une structure funiculaire). L’ingénieur suisse Heinz Isler était connu pour la conception de structures optimales en béton à partir de modèles en filets fins gelés pour former une coque (Chilton s. d.; Shepherd et Chilton s. d.).
On cherche donc dans cet exercice à trouver des structures funiculaires en traction en minimisant l’énergie potentielle d’un câble.
2.1 Fonction à minimiser
Un câble au repos a une énergie potentielle minimale. On suppose que la forme du câble est donnée par la fonction \(y: [0, L] \rightarrow \mathbb R\), et que le câble est attaché aux points \(A = (0, 0)\) et \(B = (L, H)\) à ses extrémités. On a deux contributions à l’énergie potentielle: l’énergie potentielle de pesanteur et l’énergie élastique de déformation du câble (on suppose le câble extensible). La première contribution est caractérisée par \(\mu\), la masse linéique constante du câble, et \(g\) l’accélération de pesanteur. La seconde contribution est caractérisée par \(k\), la raideur totale du câble, et \(L_0\) la longueur initiale du câble.
On prend \(x \in [0, L]\) la coordonnée horizontale et \(s \neq x\) la coordonnée le long du câble. La masse totale du câble et la longueur du câble sont données par:
\[ m = \int_A^B \mu \,\mathrm ds,\quad \ell = \int_A^B \mathrm ds\]
L’énergie potentielle de pesanteur de chaque morceau de câble \(\mathrm ds\) est donnée par \(g\mu y\,\mathrm ds\). L’énergie potentielle totale est donnée par:
\[ U = \int_A^B g\mu y\,\mathrm ds + \frac{1}{2}k(\ell - L_0)^2 \]
Consigne
En utilisant l’identité \(\mathrm ds = \sqrt{\mathrm dx^2 + \mathrm dy^2}\), montrer que \(U\) est donné par:
\[ \frac{U}{g\mu} = \int_0^L y\sqrt{1 + y'^2}\,\mathrm dx + \frac{1}{2}\alpha\left(\int_0^L \sqrt{1 + y'^2}\,\mathrm dx - L_0\right)^2,\quad \text{avec } \alpha = \frac{k}{g\mu} \]
Comme minimiser \(U/(g\mu)\) est identique à minimiser \(U\), on a le problème de minimisation suivant:
\[ \begin{aligned} \min_y \quad& \int_0^L y\sqrt{1 + y'^2}\,\mathrm dx + \frac{1}{2}\alpha\left(\int_0^L \sqrt{1 + y'^2}\,\mathrm dx - L_0\right)^2 \\ \text{sous contraintes} \quad& y(0) = 0\text{ et }y(L) = H \end{aligned} \]
2.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.
Consigne
On note les deux intégrales dans l’énergie:
\[U_p(y) = \int_0^L y\sqrt{1 + y'^2}\,\mathrm dx\quad\text{et}\quad \ell(y) = \int_0^L \sqrt{1 + y'^2}\,\mathrm dx \]
Montrer que :
\[U_p(y_h) = \sum_{i=1}^{N-1} \Delta x \frac{y_i + y_{i+1}}{2}\sqrt{1 + \left(\frac{\Delta y_i}{\Delta x}\right)^2}\] et \[\ell(y_h) = \sum_{i=1}^{N-1} \Delta x \sqrt{1 + \left(\frac{\Delta y_i}{\Delta x}\right)^2}\]
Implémentation à \(N\) degrés de liberté
On va directement implémenter la fonction à minimiser pour \(N\) degrés de liberté en utilisant les expressions ci-dessus:
Pour cela, on va fixer la valeur des paramètres:
| Paramètre | Valeur |
|---|---|
| \(L\) | 10 |
| \(H\) | 3 |
| \(L_0\) | 15 |
| \(\alpha\) | 10 |
Consigne
Avec les fonctions jnp.sqrt, jnp.sum, implémenter la fonction qui calcule l’énergie potentielle du câble:
L0 = 15
L = 10
alpha = 10
H = 3
def longueur_nddl(y, x):
return ...
def energie_nddl(y, x):
return ...2.3 Optimisation à deux degrés de liberté
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.
@np.vectorize
def longueur_2ddl(y2, y3):
x = jnp.array([0, L/2, L])
y = jnp.array([0, y2, y3])
return longueur_nddl(y, x)
@np.vectorize
def energie_2ddl(y2, y3):
x = jnp.array([0, L/2, L])
y = jnp.array([0, y2, y3])
return energie_nddl(y, x)
Consigne
On souhaite afficher sur un graphe les courbes de niveau de \(U(y_2, y_3)\) et la courbe de la contrainte \(y_3 = H\), et chercher graphiquement le point tangent où se trouve l’optimum.
En utilisant le code ci-dessous, afficher la courbe des valeurs de \(y_2\) et \(y_3\) pour lesquelles \(\ell = L_0\). Changer la valeur de \(\alpha\) à 100, commenter la position de l’optimum.
ax.contour(y2, y3, longueur_2ddl(y2, y3), levels=[L0], colors='k', ls='--')
Solution
Consigne
Afficher la ou les forme(s) optimale(s) du câble pour \(\alpha = 10\).
Solution
2.4 Optimisation à \(N\) degrés de liberté
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 & U_p(\mathbf y) + \frac{1}{2}\alpha(\ell(\mathbf y) - L_0)^2 = \sum_{i=1}^{N-1} \Delta x \frac{y_i + y_{i+1}}{2}\sqrt{1 + \left(\frac{\Delta y_i}{\Delta x}\right)^2} + \frac{1}{2}\alpha\left(\sum_{i=1}^{N-1} \Delta x \sqrt{1 + \left(\frac{\Delta y_i}{\Delta x}\right)^2} - L_0\right)^2 \\ \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 & U_p(\mathbf y) + \frac{1}{2}\alpha(\ell(\mathbf y) - L_0)^2 \\ \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)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)
energie_nddl(y_initial, x)Array(135., dtype=float32)On peut enfin appeler minimize:
solution = opt.minimize(energie_nddl, y_initial, args=(x,), jac=grad(energie_nddl), constraints=(bords,), 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
2.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.xChoisissez 5 valeurs de \(N\) entre 2 et 50 pour et afficher chaque solution.
Solution
2.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\).
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\ & U(y) = U_p(y) +\frac{1}{2}\alpha (\ell(y) - L_0)^2 \\\text{sous contraintes } \quad & y(0) = 0\text{ et }y(L) = H \\ & h(x) \leq y(x) \end{aligned}\]
À l’aide de opt.LinearConstraint, implémenter la contrainte d’obstacle pour \(h(x) = -H\). Tester avec \(N=20\) et afficher la solution optimale.
Solution
Références
Chilton, John. s. d. « Heinz Isler: Natural Hills on Different Edge Lines ». https://drawingmatter.org/heinz-isler-natural-hills-on-different-edge-lines?page&name=heinz-isler-natural-hills-on-different-edge-lines. Consulté le 14 août 2026.
Shepherd, Paul, et John Chilton. s. d. « Heinz Isler Model ». https://drawingmatter.org/heinz-isler-model?page&name=heinz-isler-model. Consulté le 24 août 2026.