À 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$.
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. $$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:
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...
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. $$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:
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()
TODO
Quelques livres sur le sujet:
Quelques liens intéressants: