3  Optimisation d’une poutre en treillis

\(\newcommand{\pvec}{{\{\mathbf p\}}}\) \(\newcommand{\Fvec}{\{\mathbf F\}}\) \(\newcommand{\Kmat}{[\mathbf K(\pvec)]}\) \(\newcommand{\uvec}{{\{\mathbf u\}}}\) \(\newcommand{\upvec}{\{\mathbf u(\pvec)\}}\) \(\newcommand{\dp}{\mathrm d_{\mathbf p}}\) \(\newcommand{\drp}{\partial_{\mathbf p}}\) \(\newcommand{\du}{\mathrm d_{\mathbf u}}\) \(\newcommand{\dru}{\partial_{\mathbf u}}\) \(\newcommand{\Lcal}{\mathcal L}\) \(\newcommand{\lvec}{\{\mathbf \lambda\}}\)

Dans cet exercice, on cherche à optimiser des structures en treillis (structures composées de bielles et tirants) pour minimiser la déformation sous un chargement donné.

Figure 3.1: Illustration de treillis optimal pour un porte-à-faux

3.1 Rappels éléments finis

On considère un treillis composé de \(N\) nœuds dont les positions sont donnés par le vecteur \(\pvec\). Le vecteur \(\Fvec\) défini les forces appliquées à chaque nœud, et \(\uvec\) est le déplacement de chaque nœud. On suppose que toutes les barres sont du même matériau au module d’Young \(E\) et de section \(A\) identique pour chaque barre. On note \(\Kmat\) la matrice de raideur de la structure.

Équilibre et minimisation

L’équation d’équilibre de la structure est, pour un \(\pvec\) fixé :

\[ \Kmat\uvec = \Fvec \tag{3.1}\]

Cette équation correspond au problème d’optimisation suivant :

\[\min_{\uvec} \quad \frac{1}{2}\uvec^T\Kmat\uvec - \Fvec^T\uvec\]

Beaucoup de problèmes de mécanique peuvent s’écrire comme des problèmes de minimisation (grâce au principe de moindre action) sous contraintes. Les méthodes d’optimisation des structures sont donc utiles au-delà du génie mécanique !

Avant de traiter le problème d’optimisation, on va procéder à un calcul de la déformation du treillis initial de la Figure 3.1. Pour cela, on défini quelques fonctions pour le calcul de la matrice de raideur et l’affichage du treillis.

import jax.numpy as jnp
import matplotlib.pyplot as plt


def longueurs(positions, connectivite):
    "Calcule la longueur des barres"
    dx = positions[connectivite[:, 1]] - positions[connectivite[:, 0]]
    return jnp.linalg.norm(dx, axis=1)

def mat_rotations(positions, connectivite):
    "Calcule les matrices de rotation de chaque barre"
    l = longueurs(positions, connectivite)
    dx = positions[connectivite[:, 1]] - positions[connectivite[:, 0]]
    c, s = dx.T / l
    R = jnp.zeros([connectivite.shape[0], 2, 4])
    R = R.at[:, (0, 1), (0, 2)].set(c[:, jnp.newaxis])
    R = R.at[:, (0, 1), (1, 3)].set(s[:, jnp.newaxis])
    return R

def mat_raideur_locale(aires, E, longueurs):
    "Retourne les matrices de raideurs dans leur repère local"
    ke = jnp.array([[1, -1], [-1, 1]], dtype=float)
    return E * (aires / longueurs)[:, jnp.newaxis, jnp.newaxis] * ke[jnp.newaxis, :, :]

def mat_raideur_element(klocal, rotations):
    "Calcule T^t * k_e * T, les matrices de raideur dans le repère global"
    return jnp.einsum('ekl,eki,eij->elj', rotations, klocal, rotations)

def assemble_raideur(kelem, connectivite, ndofs):
    "Assemble les raideurs par élément dans la matrice globale"
    K = jnp.zeros([ndofs, ndofs])
    for conn, ke in zip(connectivite, kelem):
        ddl = jnp.array([2 * conn[0], 2 * conn[0] + 1,
                         2 * conn[1], 2 * conn[1] + 1])
        K = K.at[jnp.ix_(ddl, ddl)].add(ke)
    return K

def mat_raideur(positions, connectivite, E, aires):
    "Calcule la matrice de raideur de la structure libre"
    L = longueurs(positions, connectivite)
    kloc = mat_raideur_locale(aires, E, L)
    T = mat_rotations(positions, connectivite)
    kel = mat_raideur_element(kloc, T)
    return assemble_raideur(kel, connectivite, positions.size)

def blocages(K, ddl_bloques):
    "Applique les déplacements bloqués à la matrice de raideur"
    K = K.at[ddl_bloques, :].set(0).at[:, ddl_bloques].set(0)
    K = K.at[(ddl_bloques, ddl_bloques)].set(1)
    return K

def mat_raideur_bc(positions, connectivite, E, aires, ddl_bloques):
    "Retourne la matrice de raideur avec les déplacements bloqués appliqués"
    return blocages(mat_raideur(positions, connectivite, E, aires), ddl_bloques)

def plot_treillis(ax, positions, connectivite, couleur=None, **kwargs):
    "Affiche un treillis en position initiale"
    if couleur is None:
        couleur = ax.plot([], [])[0].get_color()

    for conn in connectivite:
        noeuds = positions[conn].T
        ax.plot(*noeuds, marker='o', color=couleur, **kwargs)

def plot_deformee(ax, positions, connectivite, deplacement, facteur=1.):
    "Affiche un treillis déformé (avec facteur d'échelle)"
    plot_treillis(ax, positions + facteur * deplacement.reshape(positions.shape), connectivite)

def plot_fleches(ax, positions, deplacement, **kwargs):
    "Affiche un vecteur sur chaque nœud (par exemple déplacement ou force appliquée)"
    ax.quiver(*positions.T, *deplacement.reshape(positions.shape).T, **kwargs)

La Figure 3.2 montre les coordonnées et numérotation des nœuds du treillis.

Figure 3.2: Position et numérotation des nœuds

Consigne

Remplir les tableaux positions et connectivite et vérifier le treillis avec la fonction plot_treillis.

positions = jnp.array([
    [0, 0], # position nœud 0
    [2, 0], # position nœud 1
    # ... à remplir
], dtype=float)

connectivite = jnp.array([
    [0, 1], # connectivité d'une barre (à modifier)
    # ... à remplir
], dtype=int)

fig, ax = plt.subplots()

plot_treillis(ax, positions, connectivite)
ax.set_aspect('equal')

On va d’abord faire un calcul d’équilibre avec le vecteur \(\Fvec \in \mathbb R^{10}\) et les appuis donnés à la Figure 3.2 et afficher la déformée.

Consigne

On prend \(E = 10 \text{ GPa}\), \(A = 200\text{ mm}^2\), \(F_x = 1\text{ kN}\). Construire le vecteur des forces \(\Fvec\) et des degrés de libertés bloqués:

F = np.zeros(positions.size)

F[...] = ... # composantes non nulles de F

ddl_bloques = np.array([
    # numéro des DDL bloqués
], dtype=int)

Assembler ensuite la matrice de raideur du système avec mat_raideur_bc.

K = mat_raideur_bc(...)

Puis résoudre Équation 3.1 à l’aide de jnp.linalg.solve et afficher la déformée.

fig, ax = plt.subplots()

plot_treillis(ax, positions, connectivite, couleur='k')
plot_deformee(ax, positions, connectivite, u)

Solution

Déformée du treilis

3.2 Problème d’optimisation

Nous allons chercher à minimiser l’énergie potentielle de déformation du treillis. Définissons la fonction coût en fonction de \(\uvec\) et \(\pvec\):

\[ f(\uvec, \pvec) = \uvec^T\Kmat\uvec\]

Si le vecteur \(\uvec\) est solution de Équation 3.1 (\(\uvec\) est en équilibre), alors on a:

\[ f(\uvec, \pvec) = \Fvec^T\uvec = f(\uvec)\]

On cherche à trouver la position des nœuds \(\pvec\) qui minimise \(f\), on a donc le problème de minimisation suivant:

\[\begin{aligned} \min_{\pvec} \quad & f(\uvec) = \Fvec^T\uvec \\ \text{sous contrainte}\quad & \Kmat \uvec = \Fvec \\ & [\mathbf A]\pvec = \{\mathbf b\}\end{aligned} \tag{3.2}\]

La deuxième contrainte représente le fait qu’on ne veut pas changer la position des nœuds 0, 1 et 5.

Avertissement

À première vue \(f\) ne dépend pas de \(\pvec\) : la dépendance est en réalité implicite, puisque \(\uvec\) dépend de \(\pvec\) via l’équation d’équilibre.

Note

La structure du problème 3.2 est caractéristique des problèmes d’optimisation rencontrés par les ingénieurs dans un large éventail de domaines : on cherche à minimiser un coût qui dépend directement de la solution à un problème de mécanique. Par exemple, on peut chercher à minimiser le coefficient de traînée d’une aile d’avion, qui dépend de la vitesse d’écoulement du fluide, qui elle-même est solution des équations de Navier–Stokes.

3.3 Minimisation

Afin de minimiser \(f\) par rapport à \(\pvec\), nous avons besoin de la dérivée par rapport à \(\pvec\) de \(f(\upvec)\), que l’on note \(\dp f\). En développant la dérivée, on obtient:

\[ \dp f = \du f \cdot \dp\uvec \]

La dérivée \(\du f\) est facile à calculer (c’est \(\Fvec\)), en revanche, la dérivée \(\dp \uvec\) fait intervenir la dérivée de l’inverse de la matrice de raideur \(\Kmat^{-1}\) qui est beaucoup trop complexe à calculer dans tous les cas. On pourrait faire une évaluation numérique de par différences finies \(\dp f\sim f(\{\mathbf u(\pvec + \{\delta \mathbf p\})\}) - f(\upvec)\), mais pour calculer toutes les composantes du gradient il faudrait résoudre \(2N\) problèmes aux éléments finis, ce qui est très coûteux. On se tourne donc vers une autre méthode.

3.3.1 Méthode de l’adjoint

La méthode de l’adjoint permet d’évaluer \(\dp f\) avec uniquement un calcul supplémentaire d’équilibre. La démonstration est donnée à titre indicatif pour les curieuses et curieux.

Démonstration de la méthode de l’adjoint

La dérivation ci-dessous est issue de Bradley (2024) (lien PDF).

On écrit la contrainte de 3.2 sous la forme:

\[\Kmat\upvec = \Fvec \quad\Leftrightarrow\quad g(\uvec, \pvec) = \Kmat\uvec - \Fvec = 0\]

Le problème de minimisation est donc:

\[\begin{aligned} \min_{\pvec} \quad & f(\uvec) \\ \text{sous contrainte}\quad & g(\uvec, \pvec) = 0\end{aligned} \tag{3.3}\]

Le lagrangien du problème est:

\[\Lcal(\uvec, \pvec, \lvec) = f(\uvec) - \lvec^Tg(\uvec, \pvec)\]

Le vecteur \(\lvec\) est le multiplicateur de Lagrange. Quand \(\uvec\) est solution de \(g(\uvec, \pvec) = 0\), on a \(\forall \lvec\)

\[\begin{aligned} \dp f = \dp \Lcal & = \du f\cdot \dp\uvec - \cancel{\dp \lvec^T g} - \lvec^T(\dru g\cdot \dp\uvec + \drp g)\\ & = (\du f - \lvec^T\dru g)\cdot \dp\uvec - \lvec^T \drp g \end{aligned}\]

Puisque le multiplicateur de Lagrange est arbitraire, on choisit \(\lvec\) tel que

\[ \du f - \lvec^T\dru g = 0\quad\Leftrightarrow\quad [\dru g]^T \lvec = \du f^T \tag{3.4}\]

Cela annule le terme en facteur de \(\dp\uvec\) (la dérivée compliquée), et il reste:

\[ \dp f = -\lvec^T \drp g \]

Le système d’équations 3.4 est appelé “système adjoint”, et la solution \(\lvec\) est l’adjoint (discret). En général, l’adjoint est aussi facile (voire plus facile si \(g\) est une fonction non-linéaire) à calculer que \(\uvec\). La dérivée \(\dp g\) est parfois difficile mais bien plus facile que \(\dp \uvec\).

En appliquant à notre cas, le système adjoint est:

\[ \Kmat \lvec = \Fvec \]

La dérivée \(\drp g\) est obtenue en dérivant \(\Kmat\):

\[\drp g = (\dp \Kmat)\uvec\]

On cherche un vecteur \(\lvec\), l’adjoint, solution du système d’équations suivant:

\[ \Kmat \lvec = \Fvec\]

Puis on peut obtenir le gradient de \(f\) par rapport à \(\pvec\):

\[\dp f = -\lvec^T (\dp \Kmat)\uvec \]

Astuce

L’évaluation de \(\dp f\) est utile pour l’optimisation des structures mais aussi pour l’analyse de sensibilité, par exemple pour savoir si la présence d’un défaut affecte la performance d’une pièce.

Consigne

Calculer \(\dp f\) à l’aide de la méthode de l’adjoint et l’afficher à l’aide de plot_fleches.

On pourra se servir de la fonction jacobian de JAX pour calculer \(\dp \Kmat\):

from jax import jacobian

def mat_derivee_raideur(positions):
    return jacobian(lambda p: mat_raideur_bc(p.reshape(-1, 2), connectivite, E, A, ddl_bloques))(positions.ravel()).T

Solution

Gradient \(\dp f\) avec plot_fleches(ax, positions, gradient, scale=.005)

3.3.2 Descente de gradient

Muni du gradient \(\dp f\), nous pouvons nous servir de la propriété que \(-\dp f\) est la direction de plus grande décroissance de \(f\) par rapport à \(\pvec\) pour minimiser \(f\). Une façon naturelle de procéder est de mettre à jour de manière itérative les positions des nœuds à l’aide du gradient. Soit \(\pvec^n\) la position des nœuds actuelle, si \(\dp f \neq 0\) on peut supposer que:

\(\pvec^{n+1} = \pvec^n - \alpha \dp f(\pvec^n)\)

est une meilleure position des nœuds du treillis. Le scalaire \(\alpha\) est la taille du “pas” que l’on fait dans la direction de la plus forte pente. Cet algorithme est appelé descente de gradient, et l’on peut arrêter les itérations lorsque \(\|\dp f\| < \epsilon\) où \(\epsilon\) est la tolérance (par exemple \(10^{-3}\)).

import numpy as np
import jax.numpy as jnp

import matplotlib.pyplot as plt
from PIL import Image
from jax import jacobian

# Module d'Young
E = 1.0

positions = jnp.array([
    [0, 0],
    [2, 0],
    [0, 2.5],
    [1, 2.5],
    [2, 2.5],
    [1, 5],
], dtype=float)

connectivite = jnp.array([
    [0, 2],
    [0, 3],
    [2, 3],
    [3, 1],
    [3, 4],
    [1, 4],
    [2, 5],
    [5, 4],
], dtype=int)

aires = jnp.ones(connectivite.shape[0], dtype=float)

ddl_bloques = jnp.array([
    0, 1, 2, 3
], dtype=int)

contraintes = jnp.array([
    0, 1, 2, 3, 11,
], dtype=int)

def longueurs(positions, connectivite):
    dx = positions[connectivite[:, 1]] - positions[connectivite[:, 0]]
    return jnp.linalg.norm(dx, axis=1)

def mat_rotations(positions, connectivite):
    l = longueurs(positions, connectivite)
    dx = positions[connectivite[:, 1]] - positions[connectivite[:, 0]]
    c, s = dx.T / l
    R = jnp.zeros([connectivite.shape[0], 2, 4])
    R = R.at[:, (0, 1), (0, 2)].set(c[:, jnp.newaxis])
    R = R.at[:, (0, 1), (1, 3)].set(s[:, jnp.newaxis])
    return R

def mat_raideur_locale(aires, E, longueurs):
    ke = jnp.array([[1, -1], [-1, 1]], dtype=float)
    return E * (aires / longueurs)[:, jnp.newaxis, jnp.newaxis] * ke[jnp.newaxis, :, :]

def mat_raideur_element(klocal, rotations):
    return jnp.einsum('ekl,eki,eij->elj', rotations, klocal, rotations)

def assemble_raideur(kelem, connectivite, ndofs):
    K = jnp.zeros([ndofs, ndofs])
    for conn, ke in zip(connectivite, kelem):
        ddl = jnp.array([2 * conn[0], 2 * conn[0] + 1,
                         2 * conn[1], 2 * conn[1] + 1])
        K = K.at[jnp.ix_(ddl, ddl)].add(ke)
    return K

def mat_raideur(positions, connectivite, E, aires):
    L = longueurs(positions, connectivite)
    kloc = mat_raideur_locale(aires, E, L)
    T = mat_rotations(positions, connectivite)
    kel = mat_raideur_element(kloc, T)
    return assemble_raideur(kel, connectivite, positions.size)

def blocages(K, ddl_bloques):
    K = K.at[ddl_bloques, :].set(0).at[:, ddl_bloques].set(0)
    K = K.at[(ddl_bloques, ddl_bloques)].set(1)
    return K

def mat_raideur_bc(positions, connectivite, E, aires, ddl_bloques):
    return blocages(mat_raideur(positions, connectivite, E, aires), ddl_bloques)

K = mat_raideur_bc(positions, connectivite, E, aires, ddl_bloques)

F = jnp.zeros(positions.size)
F = F.at[10].set(1)

u = jnp.linalg.solve(K, F)
u

Array([ 0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
        1.0260627e+01,  3.1249998e+00,  9.7606268e+00,  0.0000000e+00,
        1.0260627e+01, -3.1249998e+00,  2.7833736e+01, -1.4991825e-07],      dtype=float32)

def plot_treillis(ax, positions, connectivite, couleur='k', **kwargs):
    for conn in connectivite:
        noeuds = positions[conn].T
        ax.plot(*noeuds, marker='o', color=couleur, **kwargs)

def plot_deformee(ax, positions, connectivite, deplacement, facteur=1.):
    plot_treillis(ax, positions + facteur * deplacement.reshape(positions.shape), connectivite, couleur='C0')

def plot_fleches(ax, positions, deplacement, **kwargs):
    ax.quiver(*positions.T, *deplacement.reshape(positions.shape).T, **kwargs)

fig, ax = plt.subplots()
plot_treillis(ax, positions, connectivite)
plot_deformee(ax, positions, connectivite, u, facteur=0.1)
#plot_fleches(ax, positions, u)

F.T @ u

Array(27.833736, dtype=float32)

K_prime = jacobian(mat_raideur_bc)

gradf = - u.T @ K_prime(positions, connectivite, E, aires, ddl_bloques).T @ u
gradf = gradf.T
gradf = gradf.flatten().at[contraintes].set(0)

fig, ax = plt.subplots()
plot_treillis(ax, positions, connectivite)
#plot_deformee(ax, positions, connectivite, -gradf, facteur=0.05)
ax.set_aspect('equal')
plot_fleches(ax, positions, -gradf.T)

def gradient_f(positions, u):
    gradf = - u.T @ K_prime(positions, connectivite, E, aires, ddl_bloques).T @ u
    gradf = gradf.T
    return gradf.flatten().at[contraintes].set(0)

fig, axs = plt.subplots(1, 2)

step = 0.05

pos = positions.copy()
plot_treillis(axs[0], pos, connectivite, couleur=f"C0")

for i in range(100):
    K = mat_raideur_bc(pos, connectivite, E, aires, ddl_bloques)
    u = jnp.linalg.solve(K, F)
    gradf = gradient_f(pos, u)
    pos -= step * gradf.reshape(pos.shape)

plot_treillis(axs[0], pos, connectivite, couleur=f"C{i}")



im = np.asarray(Image.open("figs/treillis_prager.png"))
axs[1].imshow(im)
axs[1].spines[:].set_visible(False)
axs[1].set_xticks([])
axs[1].set_yticks([])
axs[0].set_aspect('equal')

#fig.savefig("treillis.png")

Bradley, Andrew M. 2024. « PDE-constrained Optimization and the Adjoint Method ». Stanford University.