La méthode des multiplicateurs de Lagrange

À quoi ça sert ?

À trouver les extremums (minimums, maximums) d'une fonction $f$ d'une ou plusieurs variables $x_1, \dots, x_n$, sous réserve que l'ensemble solution respecte un contrainte d'égalité: $g(x_1, \dots, x_n) = 0$. Autrement dit la méthode des multiplicateurs de Lagrange va permettre de résoudre certains problèmes d'optimisation sous contraintes.

Exemple: maximiser $f(x_1,x_2)$ soumise aux contraintes $g(x_1,x_2)=0$.

Cas d'une fonction à une variable

Dans ce cas, on écrit simplement $x := x_1$.

On pose la fonction:

$$ \mathcal{L}(x, \lambda) = f(x) + \lambda g(x) $$

$\lambda$ est ce qu'on appelle un multiplicateur de Lagrange ; sa valeur n'est pas connue à priori.

Pour maximiser $\mathcal{L}$, on annule ses dérivées partielles (condition nécessaire du premier ordre). Le problème initial revient à résoudre le système d'équations à deux inconnues suivant:

$$ \left\{ \begin{array}{rcl} {\large \frac{\partial\mathcal{L}(x,\lambda)}{\partial x}} & = & 0 \\ {\large \frac{\partial\mathcal{L}(x,\lambda)}{\partial \lambda}} & = & 0 \end{array} \right. $$

Exemple

On cherche à résoudre le problème d'optimisation suivant:

$$ \begin{align} \max & \quad f(x) = a x^2 + b x + c \\ \text{s.t.} & \quad g(x) = a' x^2 + b' x + c' = 0 \end{align} $$

avec par exemple

$a = -\frac12$, $b = 0$, $c = 20$, $a' = 1$, $b' = 2$ et $c' = -10$

soit

$$ \begin{align} \max & \quad f(x) = -\frac12 x^2 + 20 \\ \text{s.t.} & \quad g(x) = x^2 + 2 x - 10 = 0 \end{align} $$

avec:

  • $f(x) = -\frac12 x^2 + 20$ la fonction à maximiser
  • $g(x) = x^2 + 2 x - 10 = 0$ la contrainte à respecter

La fonction de Lagrange correspondant à ce problème est:

$$ \begin{align} \mathcal{L}(x,\lambda) & = f(x) + \lambda g(x) \\ & = -\frac12 x^2 + 20 + \lambda (x^2 + 2 x - 10) \\ & = -\frac12 x^2 + 20 + \lambda x^2 + 2 \lambda x - 10 \lambda \\ & = (-\frac12 + \lambda) x^2 + 2 \lambda x + 20 - 10 \lambda \end{align} $$

Les conditions du premier ordre (annulation des dérivées premières) sont données par:

$$ \begin{align} \left\{ \begin{array}{rcl} {\large \frac{\partial\mathcal{L}(x,\lambda)}{\partial x}} & = & 0 \\ {\large \frac{\partial\mathcal{L}(x,\lambda)}{\partial \lambda}} & = & 0 \end{array} \right. & \Leftrightarrow \left\{ \begin{array}{rcl} 2 (-\frac12 + \lambda) x + 2 \lambda & = & 0 \\ x^2 + 2 x - 10 & = & 0 \end{array} \right. \\ & \\ & \Leftrightarrow \left\{ \begin{array}{rcl} -x + 2 \lambda x + 2 \lambda & = & 0 \\ x & = & -\sqrt{11} - 1 \quad \text{ou} \quad x = \sqrt{11} - 1 \end{array} \right. \\ & \\ & \Leftrightarrow \left\{ \begin{array}{rcl} 2 \lambda x + 2 \lambda & = & x \\ x & = & -\sqrt{11} - 1 \quad \text{ou} \quad x = \sqrt{11} - 1 \end{array} \right. \\ & \\ & \Leftrightarrow \left\{ \begin{array}{rcl} \lambda & = & \frac{x}{2x + 2} \\ x & = & -\sqrt{11} - 1 \quad \text{ou} \quad x = \sqrt{11} - 1 \end{array} \right. \\ & \\ & \Leftrightarrow \left\{ \begin{array}{rcl} \lambda & = & \frac{-\sqrt{11} - 1}{2(-\sqrt{11} - 1) + 2} \\ x & = & -\sqrt{11} - 1 \end{array} \right. \quad \text{ou} \quad \left\{ \begin{array}{rcl} \lambda & = & \frac{\sqrt{11} - 1}{2(\sqrt{11} - 1) + 2} \\ x & = & \sqrt{11} - 1 \end{array} \right. \\ & \\ & \Leftrightarrow \left\{ \begin{array}{rcl} \lambda & = & \frac{\sqrt{11} + 1}{2 \sqrt{11}} \\ x & = & -\sqrt{11} - 1 \end{array} \right. \quad \text{ou} \quad \left\{ \begin{array}{rcl} \lambda & = & \frac{\sqrt{11} - 1}{2 \sqrt{11}} \\ x & = & \sqrt{11} - 1 \end{array} \right. \end{align} $$

In [ ]:
# Check results

#def deriv_partielle_x(x, _lambda):
#    return 2. * (-0.5 + _lambda) * x + 2. * _lambda

#def deriv_partielle_lambda(x, _lambda):
#    return x**2. + 2. * x - 10.

#print((math.sqrt(11.) + 1.)/(2. * math.sqrt(11)))
#print((math.sqrt(11.) - 1.)/(2. * math.sqrt(11)))

#print(-math.sqrt(11.) -1.)
#print(math.sqrt(11.) - 1.)

#print(g_roots, lambda_roots)

#print(deriv_partielle_x(g_roots[0], lambda_roots[0]))      # ERR  TODO
#print(deriv_partielle_lambda(g_roots[0], lambda_roots[0])) # OK

#print(deriv_partielle_x(g_roots[1], lambda_roots[1]))      # ERR  TODO
#print(deriv_partielle_lambda(g_roots[1], lambda_roots[1])) # OK

# TODO: REVOIR LES CALCULS POUR DERIV_LAMBDA ET LAMBDA_ROOTS, LES RESULTATS NE VERIFIENT PAS LE SYSTEME...

TODO

Bon, ok, l'exemple 1D est un peu particulier... On se retrouve à devoir choisir entre 2 solutions, les 2 racines de $g$ et on pouvait très bien faire ça dés le début sans avoir recours aux multiplicateurs de Lagrange...

Cet exemple illustratif - dont la motivation initiale était de pouvoir représenter la fonction $\mathcal{L}$ et ses points stationnaires - n'est peut-être pas très pertinent du coup... Sauf à montrer en quoi les multiplicateurs de Lagrange deviennent réellement utiles en 2D et plus...

TODO: ajouter un exemple 2D : $\max f(x) = -x^2 + 10$ s.t. $g(x) = x_1^2 + 2 x_2^2 - 1 = 0$ : une ellipse centrée sur 0 (intersection entre f et un plan (incliné)).


In [ ]:
#%matplotlib notebook
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    a = -0.5
    b = 0.
    c = 20.
    return a * x**2 + b * x + c
    
class ConstraintFunction:
    def __init__(self):
        self.a = 1.
        self.b = 2.
        self.c = -10.
        
    def __call__(self, x):
        return self.a * x**2 + self.b * x + self.c
    
    def delta(self):
        return self.b ** 2. - 4 * self.a * self.c
    
    def roots(self):
        return np.array([
                         (-g.b + math.sqrt(g.delta()))/(2. * g.a),
                         (-g.b - math.sqrt(g.delta()))/(2. * g.a)
                        ])

g = ConstraintFunction()
lambda_roots = np.array([
                         (math.sqrt(11.) + 1.)/(2. * math.sqrt(11.)),
                         (math.sqrt(11.) - 1.)/(2. * math.sqrt(11.))
                        ])

g_roots = g.roots()

x_min = -10.
x_max = 10.

lambda_min = -3.       # TODO
lambda_max = 3.        # TODO

In [ ]:
# Build datas ###############

x = np.linspace(x_min, x_max, 200)

fx = f(x)
gx = g(x)

lambda_ = 1.
lx = fx + lambda_ * gx

# Plot data #################

fig = plt.figure(figsize=(8.0, 8.0))
ax = fig.add_subplot(111)

ax.plot(x, fx, "-b", label=r"$f(x)$")
ax.plot(x, gx, "-r", label=r"$g(x)$")
ax.plot(x, lx, "--k", label=r"$\mathcal{L}(x)$")

ax.plot(np.array([g_roots, g_roots, np.zeros(len(g_roots))]),
        np.array([np.zeros(len(g_roots)), f(g_roots), f(g_roots)]), 
        ".:r",
        lw=1)
        #label=r"$g(x) = 0$")

ax.hlines(y=0, xmin=x[0], xmax=x[-1], colors='k')
ax.vlines(x=0, ymin=f(x[0]), ymax=g(x[-1]), colors='k')

# Set title and labels ######

ax.set_title(r"Graphical representation of $f$, $g$ and $\mathcal{L}$", fontsize=18)
ax.set_xlabel(r"$x$", fontsize=18)
#ax.set_ylabel(r"$f(x)$", fontsize=18)

# Set legend ################

ax.legend(loc='upper right', fontsize=18)

# Plot ######################

plt.show()

In [ ]:
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
#from matplotlib.colors import LightSource

# Build datas ###############

N = 100

x = np.linspace(x_min, x_max, N)
lambda_ = np.linspace(lambda_min, lambda_max, N)

xx,ll = np.meshgrid(x, lambda_)

fx = f(x)
gx = g(x)

z = fx + ll * gx

# Plot data #################

fig = plt.figure(figsize=(10.0, 10.0))

ax = axes3d.Axes3D(fig)
#surf = ax.plot_surface(xx, ll, z, cmap=cm.jet, rstride=1, cstride=1, color='b', shade=False, alpha=1.)
#fig.colorbar(surf, shrink=0.5, aspect=5)

#ax = fig.gca(projection='3d')
#ax.plot_surface(xx, ll, z, rstride=5, cstride=5, alpha=0.3)
#cset = ax.contourf(xx, ll, z, zdir='z', offset=0, cmap=cm.coolwarm)

#ax = axes3d.Axes3D(fig)
#ax.plot_wireframe(xx, ll, z, color='k', alpha=0.3, cstride=0, rstride=1)
ax.plot_wireframe(xx, ll, z, color='k', alpha=0.3, cstride=0, rstride=2)
#ax.plot_wireframe(xx, ll, z, color='k', alpha=0.3, cstride=1, rstride=0)


lambda_0 = [0]
xx,ll0 = np.meshgrid(x, lambda_0)
z0 = fx + ll0 * gx
ax.plot_wireframe(xx, ll0, z0, color='r', alpha=1., cstride=0, rstride=1)


#x0 = np.array([root1])
#fx0 = f(x0)
#gx0 = g(x0)
#xx0,ll = np.meshgrid(x0, lambda_)
#z0 = fx0 + ll * gx0
#ax.plot_wireframe(xx0, ll, z0, color='g', alpha=1., cstride=1, rstride=0)


x0 = np.array(g_roots)
fx0 = f(x0)
gx0 = g(x0)
xx0,ll = np.meshgrid(x0, lambda_)
z0 = fx0 + ll * gx0
ax.plot_wireframe(xx0, ll, z0, color='g', alpha=1., cstride=1, rstride=0)

ax.scatter(g_roots,                                # x
           lambda_roots,                           # y
           f(g_roots) + lambda_roots * g(g_roots), #z
           label="points stationnaires")

ax.legend(loc='upper right', fontsize=14)

#fig, ax = plt.subplots(figsize=(10.0, 10.0), subplot_kw=dict(projection='3d'))
#ls = LightSource(270, 45)
## To use a custom hillshading mode, override the built-in shading and pass
## in the rgb colors of the shaded surface calculated from "shade".
#rgb = ls.shade(z, cmap=cm.gist_earth, vert_exag=0.1, blend_mode='soft')
#surf = ax.plot_surface(xx, ll, z, rstride=1, cstride=1, facecolors=rgb,
#                       linewidth=0, antialiased=False, shade=True)


ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$\lambda$')
ax.set_zlabel(r'$\mathcal{L}(x)$');

Les points stationnaires de $\mathcal{L}$ sont quelquepart sur les droites vertes...

Cas d'une fonction à deux variables

TODO: reécrire cette partie

On pose la fonction:

$$ \mathcal{L}(x_1,x_2,\lambda) = f(x_1,x_2) + \lambda g(x_1,x_2) $$

$\lambda$ est ce qu'on appelle un multiplicateur de Lagrange ; sa valeur n'est pas connue à priori.

Pour maximiser $\mathcal{L}$, on annule ses dérivées partielles (condition nécessaire du premier ordre). Le problème initial revient à résoudre le système d'équations à trois inconnues suivant:

$$ \left\{ \begin{array}{rcl} {\large \frac{\partial\mathcal{L}(x_1,x_2,\lambda)}{\partial x_1}} & = & 0 \\ {\large \frac{\partial\mathcal{L}(x_1,x_2,\lambda)}{\partial x_2}} & = & 0 \\ {\large \frac{\partial\mathcal{L}(x_1,x_2,\lambda)}{\partial \lambda}} & = & 0 \end{array} \right. $$

Exemple

On cherche le rectangle d'aire maximum et de périmètre constant $P$.

Plus formellement, on cherche à résoudre le problème d'optimisation suivant:

$$ \begin{align} \max & \quad f(x_1,x_2) = x_1 x_2 \\ \text{s.t.} & \quad g(x_1,x_2) = 2x_1 + 2x_2 - P = 0 \end{align} $$

avec:

  • $x_1$ et $x_2$ les dimensions du rectangle (respectivement la largeur et la hauteur)
  • $f(x_1,x_2) = x_1 x_2$ l'aire du rectangle (la fonction à maximiser)
  • $g(x_1,x_2) = 2x_1 + 2x_2 - P = 0$ la contrainte sur le périmètre du rectangle ($2x_1 + 2x_2 = P$)

La fonction de Lagrange correspondant à ce problème est:

$$ \mathcal{L}(x_1,x_2,\lambda) = x_1 x_2 + \lambda (2x_1 + 2x_2 - P) $$

Les conditions du premier ordre (annulation des dérivées premières) sont données par:

$$ \begin{align} \left\{ \begin{array}{rcl} {\large \frac{\partial\mathcal{L}(x_1,x_2,\lambda)}{\partial x_1}} & = & 0 \\ {\large \frac{\partial\mathcal{L}(x_1,x_2,\lambda)}{\partial x_2}} & = & 0 \\ {\large \frac{\partial\mathcal{L}(x_1,x_2,\lambda)}{\partial \lambda}} & = & 0 \end{array} \right. & \Leftrightarrow \left\{ \begin{array}{rcl} x_2 + 2 \lambda & = & 0 \\ x_1 + 2 \lambda & = & 0 \\ 2x_1 + 2x_2 - P & = & 0 \end{array} \right. \\ & \\ & \Leftrightarrow \left\{ \begin{array}{rcl} \color{red} \lambda & \color{red} = & {\color{red} {\large \frac{-x_2}{2}}} \\ \color{red} \lambda & \color{red} = & {\color{red} {\large \frac{-x_1}{2}}} \\ 2x_1 + 2x_2 - P & = & 0 \end{array} \right. \\ & \\ & \Leftrightarrow \left\{ \begin{array}{rcl} x_1 & = & x_2 \\ 2x_1 + 2x_2 - P & = & 0 \end{array} \right. \\ & \\ & \Leftrightarrow \left\{ \begin{array}{rcl} x_1 & = & x_2 \\ P & = & 2x_1 + 2x_1 = 4x_1 = 4x_2 \end{array} \right. \end{align} $$

On a donc ${\large \frac{-x_2}{2}} = {\large \frac{-x_1}{2}}$ c'est à dire $x_1 = x_2$. Ainsi, le carré ($x_1 = x_2$) est le rectangle d'aire maximum pour un périmètre donné $P$.

Remarque: il est préférable d'éliminer $\lambda$ dés le début des calculs car c'est une variable auxiliaire dont la valeur n'est pas utile.


In [ ]:
#%matplotlib notebook
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt

# Build datas ###############

P = 10
x = np.linspace(0., 5.0, 50)
y = np.linspace(0., 5.0, 50)
X, Y = np.meshgrid(x, y)

Z1 = X * Y
Z2 = 2 * X + 2 * Y

# Plot data #################

fig, ax = plt.subplots(figsize=(8,8))

# SURFACE ###################

ax.imshow(Z1,
          origin='lower',
          extent=(0,5,0,5),
          alpha=0.5,
          cmap='Blues_r')

max_value = np.max(Z1)
levels = np.array([0.1*max_value, 0.3*max_value, 0.6*max_value])

cs = plt.contour(x, y, Z1, levels,
                 linewidths=(2, 2, 3), linestyles=('dotted', 'dashed', 'solid'),
                 alpha=0.5, colors='blue')
ax.clabel(cs, inline=False, fontsize=12)

# Set legend ################

lines = [ cs.collections[0]]
labels = ['surface']

# PERIMETRE #################

levels = np.array([P])

cs = plt.contour(x, y, Z2, levels,
                 linewidths=(2, 2, 3), linestyles=('dotted', 'dashed', 'solid'),
                 alpha=0.5, colors='red')
ax.clabel(cs, inline=False, fontsize=12)

# Set title and labels ######

ax.axis('equal')              # <- SAME SCALE ON X AND Y
ax.set_title("Example for P = " + str(P), fontsize=20)
ax.set_xlabel(r"$X_1$", fontsize=20)
ax.set_ylabel(r"$X_2$", fontsize=20)

# Set legend ################

lines.append(cs.collections[0])
labels.append('périmètre P = ' + str(P))

ax.legend(lines, labels, prop={'size': 14}, loc='best', fancybox=True, framealpha=0.5)

# The optimal point #########

ax.plot([0, 5], [0, 5], ":k", alpha=0.5)
ax.plot([2.5], [2.5], "xk")
ax.text(2.4, 2.2, r"$x^*$", fontsize=14)

# Plot ######################

plt.grid()
plt.show()

Cas générale

TODO

Bibliographie

Quelques livres sur le sujet:

  • Optimisation et contrôle des systèmes linéaires (chapitre 3) de Maïtine Bergounioux aux editions Dunod
  • Toutes les mathématiques et les bases de l'informatique (p. 566-567) de Horst Stöcker aux editions Dunod