Correction sur la trace de l'algo
[Projet_Recherche_Operationnelle.git] / rapport / ProjetOptimRO.tex
index 6ffc1f64bc64cce5f108f9648e5acedcf70703f0..ccce4d2d1ae6ad6280595a3c715deff9bdb60f71 100644 (file)
@@ -21,6 +21,8 @@
 \usepackage{tocbibind}
 \usepackage{lmodern}
 \usepackage{enumitem}
+\usepackage{algorithm2e}
+\usepackage{algorithmic}
 
 
 %%%%%Marges & en-t\^etes
@@ -233,7 +235,11 @@ Définissons quelques notions supplémentaires de base nécessaires à la suite
  $ A \subset \mathbb{R}^n $ est un fermé $ \iff A = \overline{A} $.
 \end{Rmq}
 \begin{Def}
- Soient une fonction $ f: \mathbb{R}^n \longrightarrow \mathbb{R} $ et $ x^\ast \in \mathbb{R}^n $.
+ Soient $ f : \mathbb{R}^n \longrightarrow \mathbb{R} $ et $ S \subset \mathbb{R}^n $. On définit $ \mathrm{argmin} $ de $ f $ sur $ S $ par :
+ $$ \underset{x \in S}{\mathrm{argmin}} f(x) = \{ x \in \mathbb{R}^n \ | \ x \in S \land \forall y \in S \ f(y) \geq f(x) \} $$
+\end{Def}
+\begin{Def}
+ Soient une fonction $ f : \mathbb{R}^n \longrightarrow \mathbb{R} $ et $ x^\ast \in \mathbb{R}^n $.
  \newline
  On dit que $ f $ est continue en $ x^\ast $ si
  $$ \forall \varepsilon \in \mathbb{R}_{+}^{*} \ \exists \alpha \in \mathbb{R}_{+}^{*} \ \forall x \in \mathbb{R}^n \ \norme{x - x^\ast} \leq \alpha \implies |f(x) - f(x^\ast)| \leq \varepsilon $$
@@ -250,7 +256,7 @@ Définissons quelques notions supplémentaires de base nécessaires à la suite
  cette dérivée.
 \end{Def}
 \begin{Def}
- Soient une fonction $ f: \mathbb{R}^n \longrightarrow \mathbb{R} $
+ Soient une fonction $ f : \mathbb{R}^n \longrightarrow \mathbb{R} $
  et $ x^\ast, h \in \mathbb{R}^n $.
  \newline
  On dit que $ f $ est différentiable en $ x^\ast $ si il existe une application linéraire $ d_{x^\ast}f $ de $ \mathbb{R}^n $ dans $ \mathbb{R} $ telle que
@@ -279,7 +285,7 @@ Définissons quelques notions supplémentaires de base nécessaires à la suite
  $ \forall h \in \mathbb{R}^n \ d_{x^\ast}f(h) = \langle \nabla f(x^\ast),h \rangle = \nabla f(x^\ast)^\top h $
 \end{Rmq}
 \begin{Def}
- Soit $ f: \mathbb{R}^n \longrightarrow \mathbb{R} $ un fonction de classe $ \mathcal{C}^2 $.
+ Soit $ f: \mathbb{R}^n \longrightarrow \mathbb{R} $ une fonction de classe $ \mathcal{C}^2 $.
  On définit la matrice hessienne de $ f $ en $ x^\ast $ par :
  $$ H[f](x^\ast) =
   \begin{pmatrix}
@@ -304,8 +310,8 @@ Définissons quelques notions supplémentaires de base nécessaires à la suite
 
 \subsection{Conditions d'existence d'un extremum}
 
-On peut démontrer que $ \mathcal{C }$ est un ensemble fermé de $ \mathbb{R}^n $ si $ g $ et $ h $ sont continues.
-On peut en déduire que si $ J $ est continue, $ \mathcal{C } $ est un ensemble fermé et borné de $ \mathbb{R}^n $.
+On peut démontrer que $ \mathcal{C}$ est un ensemble fermé de $ \mathbb{R}^n $ si $ g $ et $ h $ sont continues.
+On peut en déduire $ \mathcal{C} $ est un ensemble fermé et borné de $ \mathbb{R}^n $.
 \begin{Th}[Théorème de Weierstrass]
  Soient $ \mathcal{C} \neq \emptyset \subset \mathbb{R}^n $ un fermé borné et $ f : \mathcal{C} \longrightarrow \mathbb{R} $ une fonction continue.
  \newline
@@ -314,7 +320,7 @@ On peut en déduire que si $ J $ est continue, $ \mathcal{C } $ est un ensemble
  \newline
  De la même façon, il existe un maximum global de $ J $ sur $ \mathcal{C} $.
 \end{Th}
-On en déduit que $ \mathcal{P} $ admet au moins une solution dans le cas où $ J, g ,h $ sont continues \cite{LJK,RON}. L'étude de la convexité de $ J $ sur $ \mathcal{C} $ permet d'explorer l'unicité de la solution \cite{LJK,RON}.
+Si $ J $ est continue, on en déduit que $ \mathcal{P} $ admet au moins une solution dans le cas où $ J, g ,h $ sont continues \cite{LJK,RON}. L'étude de la convexité de $ J $ sur $ \mathcal{C} $ permet d'explorer l'unicité de la solution \cite{LJK,RON}.
 
 \subsection{Conditions de caractérisation d'un extremum}
 
@@ -322,7 +328,7 @@ Dans le cas où $ J, g, h $ sont continûment différentiable et ses dérivées
 \newline
 On peut en déduire que une condition nécessaire et suffisante pour que $ x^\ast \in \mathring{\mathcal{C}} $ soit un des extremums locaux de $ J $ est que $ \nabla J(x^\ast) = 0 $. Mais si $ x^\ast \in \overline{\mathcal{C}}\setminus\mathring{\mathcal{C}} $ (la frontière de $ \mathcal{C} $) alors $ \nabla J(x^\ast) $ n'est pas nécessairement nul. Il sera par conséquent nécessaire de trouver d'autres caratérisations d'un extremum local \cite{FEA,WAL}.
 
-\subsubsection{Conditions de Karuch-Kuhn-Tucker}\label{KKT}
+\subsubsection{Conditions nécessaires de Karuch-Kuhn-Tucker ou \textit{KKT}}\label{KKT}
 
 \begin{Th}
  Soient $ x^\ast \in \mathbb{R}^n $, $ I = \{ 1,\ldots,p \} $ et $ J = \{ 1,\ldots,q \} $.
@@ -336,12 +342,26 @@ On peut en déduire que une condition nécessaire et suffisante pour que $ x^\as
  et
  $$ \forall i \in I \ \exists \mu_i \in \mathbb{R}_{+} \land \forall j \in J \ \exists \lambda_j \in \mathbb{R} \ \nabla J(x^\ast) + \sum_{i \in I}\mu_i{\nabla g_i(x^\ast)} + \sum_{j \in J}\lambda_j{\nabla h_j(x^\ast)} = 0 \land \forall i \in I \ \mu_i \nabla g_i(x^\ast) = 0 $$
  On appelle $ (\mu_i)_{i \in I}$ les multiplicateurs de Kuhn-Tucker et $ (\lambda_j)_{j \in J}$ les multiplicateurs de Lagrange.
+ \newline
+ On nomme également les conditions \textit{KTT} conditions nécessaires d'optimalité du premier ordre.
 \end{Th}
 \begin{proof}
  Elle repose sur le lemme de Farkas \cite{FEA,RON}.
 \end{proof}
 Il est à noter que une condition d'égalité peut se répresenter par deux conditions d'inégalité : $ \forall x \in \mathbb{R}^n \ \forall i \in \{ 1,\ldots,q \} \ h_i(x) = 0 \iff h_i(x) \leq 0 \land h_i(x) \geq 0 $ \cite{FEA}, ce qui peut permettre de réécrire le problème $ \mathcal{P} $ en éliminant les contraintes d'égalités et change la forme des conditions \textit{KKT} à vérifier mais rajoute $ 2q $ conditions d'inégalités et donc $ 2q $ multiplicateurs de Kuhn-Tucker.
+\begin{Def}
+ On appelle un point admissible $ x^\ast \in \mathcal{C} $ un point critique de $ \mathcal{P} $ si il statisfait les conditions \textit{KKT}.
+\end{Def}
+
+\subsubsection{Conditions suffisantes du deuxième ordre}
 
+\begin{Th}
+ Les conditions suffisantes en plus de celles \textit{KKT} pour que $ x^\ast \in \mathcal{C} $ soit un minimum local de $ J $ sont :
+ \begin{enumerate}[label=(\roman*)]
+  \item relâchement complémentaire dual\footnote{La définition de cette notion ne sera pas donnée car elle n'est pas nécessaire pour l'étude de la méthode PQS.} strict en $ x^\ast $.
+  \item $ \forall v \in \mathbb{R}^n \land v \neq 0 \ \langle H_x[L](x^\ast,\lambda,\mu)v,v \rangle > 0 $.
+ \end{enumerate}
+\end{Th}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -459,7 +479,7 @@ Illustrées par les méthodes de descente de gradient, aucune de ces deux strat
 \end{Def}
 \begin{Def}
  Dans le cas où $ J $ est différentiable sur $ \mathcal{C} $, on dit que un algorithme de descente converge ssi
- $$ \lim\limits_{k \rightarrow +\infty} \norme{\nabla J(x_k)} = 0 $$
+ $$ \forall x_0 \in \mathbb{R}^n \lim\limits_{k \rightarrow +\infty} \norme{\nabla J(x_k)} = 0 $$
 \end{Def}
 
 \subsubsection{Principe de démonstration de convergence}
@@ -472,7 +492,19 @@ En sommant ces inégalités pour $ k $ variant de $ 0 $ à $ N - 1 $, on obtient
 $$ \forall N \in \mathbb{N} \ J(x_0) - J(x_N) \geq c \sum_{i=0}^{N-1}\norme{\nabla J(x_i)}^2 $$
 Si $ J $ est bornée inférieurement, alors nécessairement $ J(x_0 ) - J(x_N) $ est majorée et donc la somme partielle est majorée, et donc la série $ (\sum\limits_{i=0}^{N-1}\norme{\nabla J(x_i)}^2)_{N \in \mathbb{N}} $ converge, ce qui implique :
 $$ \lim\limits_{k \rightarrow +\infty} \norme{\nabla J(x_k)} = 0 $$
-L'étude plus détaillée de différents algorithmes de descente qui utilisent différentes méthodes de recherche linéaire pour optimiser $ \varphi $ et le choix d'une direction ainsi que leurs convergences sort du cadre de ce projet.
+\begin{Def}
+ On considère $ (x_n)_{n \in \mathbb{N}} $, la suite des itérés donnés par un algorithme convergent. On note $ x^\ast $ la limite de la suite $ (x_n)_{n \in \mathbb{N}} $ et on suppose que $ x_k \neq x^\ast $, pour tout $ k \in \mathbb{N} $. La convergence de l’algorithme est alors dite :
+ \begin{itemize}
+  \item linéaire si l'erreur décroît linéairement i.e. :
+        $$ \exists \tau \in ]0,1[ \ \lim_{k \rightarrow +\infty} \frac{\norme{x_{k+1} - x^\ast}}{\norme{x_k - x^\ast}} = \tau $$
+  \item superlinéaire si :
+        $$ \lim_{k \rightarrow +\infty} \frac{\norme{x_{k+1} - x^\ast}}{\norme{x_k - x^\ast}} = 0 $$
+  \item d'ordre $ p $ si :
+        $$ \exists \tau \geq 0 \ \lim_{k \rightarrow +\infty} \frac{\norme{x_{k+1} - x^\ast}}{\norme{x_k - x^\ast}^p} = \tau $$
+        En particulier, si $ p = 2 $, la convergence est dite quadratique.
+ \end{itemize}
+\end{Def}
+L'étude plus détaillée de différents algorithmes de descente qui utilisent différentes méthodes de recherche linéaire pour optimiser $ \varphi $ ainsi que leurs convergences sort du cadre de ce projet.
 
 \section{Méthode Newtonienne}
 
@@ -484,7 +516,7 @@ En supposant $ J $ de classe $ \mathcal{C}^2 $ et la matrice hessienne $ H[J](x_
 $$ x_{k+1} = x_k - H[J](x_k)^{-1} \nabla J(x_k), $$
 où $ d_k = -H[J](x_k)^{-1} \nabla J(x_k) $ est appelée direction de Newton. La direction $ d_k $ est également l’unique solution du problème :
 $$ \underset{d \in \mathbb{R}^n}{\mathrm{argmin}} \ J(x_k) + \langle \nabla J(x_k),d \rangle + \frac{1}{2}\langle H[J](x_k)d,d \rangle $$
-Autrement dit, $ d_k $ est le point de minimum global de l’approximation de second ordre de $ J $ au voisinage du point courant $ x_k $.
+Autrement dit, $ d_k $ est le point de minimum global de l’approximation quadratique de $ J $ au voisinage du point courant $ x_k $.
 À condition que la matrice $ H[J](x_k) $ soit définie positive à chaque itération, la méthode de Newton est bien une méthode de descente à pas fixe égal à $ 1 $.
 \newline
 Les propriétés remarquables de cet algorithme sont :
@@ -521,11 +553,34 @@ Nous ne répondrons pas à ces questions qui sont hors du cadre de ce projet. Ce
 
 \section{Méthode PQS (ou SQP)}
 
-Nous supposons les fonctions $ J,g,h $ à valeurs réelles et de classe $ \mathcal{C}^1 $. Trouver une solution d’un problème d’optimisation sous contraintes fonctionnelles consiste à déterminer un point optimal $ x^\ast $ et des multiplicateurs associés $ (\lambda^\ast,\mu^\ast) $. Deux grandes familles de méthodes peuvent être définies pour la résolution des problèmes d’optimisation sous contraintes : les méthodes primales et les méthodes duales. Les approches primales se concentrent sur la détermination du point $ x^\ast $, les multiplicateurs $ (\lambda,\mu) $ ne servant souvent qu’à vérifier l’optimalité de $ x^\ast $. Les méthodes duales quant à elles mettent l’accent sur la recherche d’un multiplicateur en travaillant sur un problème d’optimisation déduit du problème initial par \textit{dualité}.
+Nous supposons les fonctions $ J,g,h $ à valeurs réelles et de classe $ \mathcal{C}^1 $. Trouver une solution d’un problème d’optimisation sous contraintes fonctionnelles consiste à déterminer un point optimal $ x^\ast $ et des multiplicateurs associés $ (\lambda^\ast,\mu^\ast) $. Deux grandes familles de méthodes peuvent être définies pour la résolution des problèmes d’optimisation sous contraintes : les méthodes primales et les méthodes duales. Les approches primales se concentrent sur la détermination du point $ x^\ast $, les multiplicateurs $ (\lambda,\mu) $ ne servant souvent qu’à vérifier l’optimalité de $ x^\ast $. Les méthodes duales quant à elles mettent l’accent sur la recherche des multiplicateurs en travaillant sur un problème d’optimisation déduit du problème initial par \textit{dualité}.
+
+\subsection{Problème quadratique sous contraintes linéaires}
+
+Nous introduisons les différentes approches développées pour la résolution des problèmes de programmation quadratique avec contraintes d'égalités et d’inégalités linéaires.
+\newline
+Ce type de problème quadratique se pose sous la forme :
+$$
+ \mathcal{PQ} \left \{
+ \begin{array}{l}
+  \displaystyle\min_{x \in \mathbb{R}^n} c^\top x + \frac{1}{2} x^\top \mathcal{Q} x \\
+  A^\top x + b \leq 0                                                                \\
+  A^{\prime^\top} x + b^\prime = 0
+ \end{array}
+ \right .
+$$
+où $$ \mathcal{Q} \in \mathcal{M}_n(\mathbb{R}) \ symétrique, c \in \mathbb{R}^n, A \in  \mathcal{M}_{n,p}(\mathbb{R}), b \in \mathbb{R}^p, A^\prime \in \mathcal{M}_{n,q}(\mathbb{R}), b^\prime \in \mathbb{R}^q $$
+Or
+$$  A^{\prime^\top} x + b^\prime = 0 \iff A^{\prime^\top} x + b^\prime \leq 0 \land   -A^{\prime^\top} x - b^\prime \leq 0 $$
+Donc le problème se ramène à :
+
+\subsubsection{Algorithme 1}
+
+\subsubsection{Algorithme 2}
 
 \subsection{Algorithmes Newtoniens}
 
-Les algorithmes newtoniens sont basés sur la linéarisation d’équations caractérisant les solutions que l’on cherche, fournies par les conditions d’optimalité d’ordre $ 1 $. Ces algorithmes sont \textit{primaux-duaux} dans le sens où ils génèrent à la fois une suite primale $ (x_k )_{k \in \mathbb{N}} $ convergeant vers une solution $ \overline{x} $ du problème considéré, et une suite duale $ (\lambda_k)_{k \in \mathbb{N}} $ (resp. $ ((\lambda_k, \mu_k))_{k \in \mathbb{N}} $) de multiplicateurs convergeant vers un multiplicateur optimal $ \overline{\lambda} $ (resp. $ (\overline{\lambda},\overline{\mu}) $) associé à $ \overline{x} $.
+Les algorithmes newtoniens sont basés sur la linéarisation d’équations caractérisant les solutions que l’on cherche, fournies par les conditions d’optimalité d’ordre $ 1 $. Ces algorithmes sont \textit{primaux-duaux} dans le sens où ils génèrent à la fois une suite primale $ (x_k )_{k \in \mathbb{N}} $ convergeant vers une solution $ \overline{x} $ du problème considéré, et une suite duale $ (\lambda_k)_{k \in \mathbb{N}} $ (resp. $ ((\lambda_k, \mu_k))_{k \in \mathbb{N}} $) de multiplicateurs convergeant vers un multiplicateur optimal $ \overline{\lambda} $ (resp. $ (\overline{\lambda},\overline{\mu})) $ associé à $ \overline{x} $.
 
 \subsection{Algorithme PQS}
 
@@ -566,7 +621,7 @@ $$ \begin{pmatrix}
   h(x_k)
  \end{pmatrix}  $$
 où $ D_h(x) $ désigne la matrice jacobienne de l’application $ h : \mathbb{R}^n \longrightarrow \mathbb{R}^q $ définie par :
-$$ D_h(x)^\top = \begin{bmatrix} \nabla h_1(x)\ldots\nabla h_q(x) \end{bmatrix} $$
+$$ D_h(x)^\top = \begin{bmatrix} \nabla h_1(x)^\top\ldots\nabla h_q(x)^\top \end{bmatrix} $$
 Posons : $ H_k = H_x[L](x_k,\lambda_k), \ d = x_{k+1} - x_k $ et $ \mu = \lambda_{k+1} $. L'itération s'écrit donc :
 $$ \begin{pmatrix}
   H_k      & D_h(x_k)^\top \\
@@ -580,7 +635,7 @@ $$ \begin{pmatrix}
  \end{pmatrix} $$
 et est bien définie à condition que la matrice $ H_x[L](x_k,\lambda_k) $ soit inversible. Ce sera le cas si :
 \begin{enumerate}[label=(\roman*)]
- \item Les colonnes $ \nabla h_1(x_k),\ldots,\nabla h_q(x_k) $ de $ D_h(x_k)^\top $ sont linéairement indépendants : c’est l’hypothèse de qualification des contraintes.
+ \item Les colonnes $ \nabla h_1(x_k)^\top,\ldots,\nabla h_q(x_k)^\top $ de $ D_h(x_k)^\top $ sont linéairement indépendants : c’est la condition première de \textit{KTT} ou condition de qualification des contraintes.
  \item Quel que soit $ d \neq 0 $ tel que $ D_h(x_k)d = 0, \ d^\top H_k d > 0 $ : c’est la condition suffisante d’optimalité du second ordre dans le cas de contraintes d’égalité.
 \end{enumerate}
 Revenons à l’itération. Elle s’écrit encore :
@@ -656,8 +711,14 @@ $$
 $$
 où $ J: \mathbb{R}^n \longrightarrow \mathbb{R} $, $g: \mathbb{R}^n \longrightarrow \mathbb{R}^p$ et $h: \mathbb{R}^n \longrightarrow \mathbb{R}^q$ sont supposées au moins différentiables.
 \newline
-Selon le même principe qu’avec contraintes d’égalité seules, on linéarise les contraintes et on utilise une approximation quadratique du Lagrangien :
+Selon le même principe qu’avec contraintes d’égalité seules, on linéarise les contraintes et on utilise une approximation quadratique du Lagrangien à l'aide de développements de Taylor-Young en $ x_k $ et $ (x_k,\lambda_k,\mu_k) $ respectivement :
 $$ L(x,\lambda,\mu) = J(x) + \lambda^\top g(x) + \mu^\top h(x), \ \lambda \in \mathbb{R}_+^p \land \mu \in \mathbb{R}^q $$
+Soit à l'ordre 2 pour le Lagrangien :
+$$ L(x,\lambda,\mu) \approx L(x_k,\lambda_k,\mu_k) + \nabla L(x_k,\lambda_k,\mu_k)^\top (x - x_k) + \frac{1}{2} (x - x_k)^\top H[L](x_k,\lambda_k,\mu_k) (x - x_k) $$
+et à l'ordre 1 pour les contraintes :
+$$ g(x) \approx g(x_k) + \nabla g(x_k)^\top(x - x_k) $$
+$$ h(x) \approx h(x_k) + \nabla h(x_k)^\top(x - x_k) $$
+En posant $ d = x - x_k $ et $ H_k = H[L](x_k,\lambda_k,\mu_k) $, on obtient le sous problème quadratique $ \mathcal{PQ}_k $ :
 
 \hrulefill
 \newline
@@ -675,7 +736,7 @@ ALGORITHME PQS AVEC CONSTRAINTES D'ÉGALITÉ ET D'INEGALITÉ.
                \mathcal{PQ}_k \left \{
                \begin{array}{l}
                 \displaystyle\min_{d \in \mathbb{R}^n} \nabla J(x_k)^\top d + \frac{1}{2}d^\top H_k d \\
-                g_j(x_k) + \nabla g_j(x_k)^\top d = 0, \ \forall j \in \{1,\ldots,p\}                 \\
+                g_j(x_k) + \nabla g_j(x_k)^\top d \leq 0 \\, \ \forall j \in \{1,\ldots,p\}                 \\
                 h_i(x_k) + \nabla h_i(x_k)^\top d = 0, \ \forall i \in \{1,\ldots,q\}
                \end{array}
                \right .
@@ -692,6 +753,299 @@ Afin que le sous-programme quadratique $ \mathcal{PQ}_k $ admette une unique sol
 \newline
 Etant une méthode newtonienne, l’algorithme PQS converge localement quadratiquement pourvu que les points initiaux  $ (x_0,\lambda_0 ) $ (resp. $ (x_0,\lambda_0,\mu_0) $) soient dans un voisinage d’un point stationnaire $ \overline{x} $ et de ses multiplicateurs associés $ \overline{\lambda} $ (resp. $ (\overline{\lambda},\overline{\mu}) $). Bien entendu, il est possible de globaliser l’algorithme en ajoutant une étape de recherche linéaire.
 
+\subsection{Stratégie d'approximation de la hessienne}
+
+\subsubsection{Équation de sécante et approximation}
+
+L'approximation $ H_k $  de la hessienne du Lagrangien peut être obtenu par la relation :
+$$ \nabla L(x_{k+1},\lambda_{k+1},\mu_{k+1}) - \nabla L(x_{k},\lambda_{k+1},\mu_{k+1}) \approx H[L](x_{k+1},\lambda_{k+1},\mu_{k+1})(x_{k+1} - x_k) $$
+On construit une approximation $ H_{k+1} $ de $ H[L](x_{k+1},\lambda_{k+1},\mu_{k+1}) $ comme solution de l’équation :
+$$ H_{k+1}(x_{k+1} - x_k) = \nabla L(x_{k+1},\lambda_{k+1},\mu_{k+1}) - \nabla L(x_{k},\lambda_{k+1},\mu_{k+1}) $$
+appelée équation de sécante ou équation de quasi-Newton.
+\newline
+De façon similaire, on peut construire une approximation $ B_{k+1} $ de $ H[L](x_{k+1},\lambda_{k+1},\mu_{k+1})^{-1} $ comme solution de l’équation :
+$$ B_{k+1}(\nabla L(x_{k+1},\lambda_{k+1},\mu_{k+1}) - \nabla L(x_{k},\lambda_{k+1},\mu_{k+1})) = x_{k+1} - x_k $$
+Dans les deux cas, les équations de quasi-Newton forment un système sous-déterminé à $ n $ équations et $ n^2 $ inconnues. Il existe donc une infinité de matrices $ H_{k+1} $ pouvant convenir.
+\newline
+Une stratégie commune est de calculer $ (x_{k+1},\lambda_{k+1},\mu_{k+1}) $ pour une matrice $ H_k $ donnée et faire une mise à jour de $ H_k $ de rang 1 ou 2 :
+$$ H_{k+1} = H_k + U_k $$
+
+\subsubsection{Mises à jour DFP et BFGS}
+
+\subsection{Exemple d'utilisation de PQS}
+
+Considérons le problème $ \mathcal{P} $ suivant :
+$$
+ \mathcal{P} \left \{
+ \begin{array}{l}
+  \displaystyle\min_{(x,y,z) \in \mathbb{R}^3} J(x,y,z) = x^2 + y^2 + z^2 -r^2       \\
+  g(x,y,z) = (g_1(x,y,z), g_2(x,y,z)) = (x^2 + y^2 - r_1^2, x^2 + z^2 -r_2^2) \leq 0 \\
+ \end{array}
+ \right .
+$$
+où $$ (r,r_1,r_2) \in \mathbb{R}_+^3. $$
+\textit{Entrées} : $ J $ et $ g $ de classe $ \mathcal{C}^2 $, $ \varepsilon = 0.01 $ la précision, $ (x_0,y_0,z_0) = $ point initial et $ (\lambda_{0_1},\lambda_{0_2}) = $ multiplicateur initial.
+\newline
+Le Lagrangien $ L $ de $ \mathcal{P} $ : $$ L((x,y,z),(\lambda_1,\lambda_2)) = x^2 + y^2 + z^2 -r^2 + \lambda_1(x^2 + y^2 - r_1^2) + \lambda_2(x^2 + z^2 -r_2^2). $$
+\newline
+Le gradient de $ J $ : $$ \nabla J(x,y,z) = (\frac{\partial J}{\partial x}(x,y,z),\frac{\partial J}{\partial y}(x,y,z),\frac{\partial J}{\partial z}(x,y,z)) = (2x,2y,2z). $$
+\newline
+Le gradient de $ g $ : $$ \nabla g(x,y,z) = (\nabla g_1(x,y,z),\nabla g_2(x,y,z)) $$
+$$ = ((\frac{\partial g_1}{\partial x}(x,y,z),\frac{\partial g_1}{\partial y}(x,y,z),\frac{\partial g_1}{\partial z}(x,y,z)),(\frac{\partial g_2}{\partial x}(x,y,z),\frac{\partial g_2}{\partial y}(x,y,z),\frac{\partial g_2}{\partial z}(x,y,z)) $$
+$$ = ((2x,2y,0),(2x,0,2z)). $$
+\newline
+Le gradient du Lagrangien $ L $ :
+$$ \nabla L((x,y,z),(\lambda_1,\lambda_2)) = \nabla J(x,y,z) + \lambda_1 \nabla g_1(x,y,z) + \lambda_2 \nabla g_2(x,y,z)) $$
+\newline
+La matrice hessienne de $ J $ : $$ H[J](x,y,z) =
+ \begin{pmatrix}
+  \frac{\partial^2 J}{\partial^2 x}(x,y,z)         & \frac{\partial^2 J}{\partial x\partial y}(x,y,z) & \frac{\partial^2 J}{\partial x\partial z}(x,y,z) \\
+  \frac{\partial^2 J}{\partial y\partial x}(x,y,z) & \frac{\partial^2 J}{\partial^2 y}(x,y,z)         & \frac{\partial^2 J}{\partial y\partial z}(x,y,z) \\
+  \frac{\partial^2 J}{\partial z\partial x}(x,y,z) & \frac{\partial^2 J}{\partial z\partial y}(x,y,z) & \frac{\partial^2 J}{\partial^2 z}(x,y,z)         \\
+ \end{pmatrix} =
+ \begin{pmatrix}
+  2 & 0 & 0 \\
+  0 & 2 & 0 \\
+  0 & 0 & 2 \\
+ \end{pmatrix} = 2Id_{\mathbb{R}^3} $$
+On en déduit que $ H[J](x,y,z) $ est inversible et que $ H[J](x,y,z)^{-1} = \frac{1}{2}Id_{\mathbb{R}^3} $.
+
+\hrulefill
+
+\subsection{Trace d'éxécution de PQS}
+
+Utilisons le problème $ \mathcal{P} $ précédent :
+
+$$
+ \mathcal{P} \left \{
+ \begin{array}{l}
+  \displaystyle\min_{(x,y,z) \in \mathbb{R}^3} J(x,y,z) = x^2 + y^2 + z^2 -r^2       \\
+  g(x,y,z) = (g_1(x,y,z), g_2(x,y,z)) = (x^2 + y^2 - r_1^2, x^2 + z^2 -r_2^2) \leq 0 \\
+ \end{array}
+ \right .
+$$
+où $$ (r,r_1,r_2) \in \mathbb{R}_+^3. $$
+\textit{Entrées} : $ J $ et $ g $ de classe $ \mathcal{C}^2 $, $ \varepsilon = 0.01 $, $ (x_0,y_0,z_0) = (100, 100 ,0)$  et $(\lambda_{0_1},\lambda_{0_2}) = (1 , 1)$, les rayons : $r= 100$  et $r1 = r2 = 10$.
+\newline
+Le Lagrangien $ L $ de $ \mathcal{P} $ : $$ L((x,y,z),(\lambda_1,\lambda_2)) = x^2 + y^2 + z^2 -r^2 + \lambda_1(x^2 + y^2 - r_1^2) + \lambda_2(x^2 + z^2 -r_2^2). $$
+\newline
+Le Lagrangien $ L $ de $ \mathcal{P} $ avec les valeurs :
+ $ L((100,100,0),(1,1)) = 100^2 + 100^2 + 0^2 -100^2 + 1 * (100^2 +100^2 - 10^2) + \lambda_2(100^2 + 100^2 -10^2). $
+ $ L((100,100,0),(1,1)) = 1000 + 1000 - 1000 + (1000 + 1000 - 100) + (1000 + 1000 -100). $
+ $ L((100,100,0),(1,1)) = 4800. $
+
+\newpage
+\begin{algorithmfloat}[#Algo 1]
+ \caption {Trace d'éxécution du PQS du problème $ \mathcal{P} $}
+ \begin{algorithmic}
+ \REQUIRE $g(x_0,y_0,z_0)\leq 0$, $(x_0,y_0,z_0) = (10, 10 ,10)$
+ \ENSURE $\min_{(x,y,z) \in \mathbb{R}^3} J(x,y,z) = x^2 + y^2 + z^2 -r^2$ and \newline $g(x,y,z) = (g_1(x,y,z), g_2(x,y,z)) = (x^2 + y^2 - r_1^2, x^2 + z^2 -r_2^2) \leq 0 $
+ \STATE \textbf{Data :}
+ \STATE $k \leftarrow 0, (x_k, y_k, z_k)  \leftarrow (100, 100, 0), r \leftarrow 100$
+ \STATE $r_1 = r2 \leftarrow 10, \varepsilon \leftarrow 0.01$
+ \STATE $\lambda_1 = \lambda_2 = 1$
+ \STATE $ H[J](x,y,z)^{-1}\leftarrow  \begin{pmatrix}
+  0.5 & 0 & 0 \\
+  0 & 0.5 & 0 \\
+  0 & 0 & 0.5 \\ \end{pmatrix} $
+\newline
+
+ \STATE{//Calcule du gradient de $ J $ :}
+ \STATE $ \nabla J(x,y,z) = (2x_k,2y_k,2z_k)$  \hfill $  // résultat : (100,100,0) $
+\newline
+ \STATE {//calcule des deux sous partie de du gradient de $ g $: }
+ \STATE $ \nabla g_1(x_a,y_a,z_a) = ((2x_a,2y_a,0)$  \hfill $ //résultat : (20, 20, 0)$
+ \STATE $ \nabla g_2(x_a,y_a,z_a) = (2x_a,0,2z_a))$  \hfill $ //résultat : (20, 0, 20)$
+ \STATE $ \nabla g(x_k,y_k,z_k) = (\nabla g_1(x_k,y_k,z_k), \nabla g_2(x_k,y_k,z_k))$
+\newline
+ \WHILE{$ (\norme{\nabla L(x_k,\lambda_k,\mu_k)} > \varepsilon $ or k $ \leq 10)$}
+
+ \STATE { //première itération :}
+
+\STATE {//Calcule du gradient de $ L $ : }
+\STATE $\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2)) = \nabla J(x_k,y_k,z_k) + \lambda_1 \nabla g_1(x_a,y_a,z_a) + \lambda_2 \nabla g_2(x_a,y_a,z_a)) $ \hfill $// résultat : (220, 220, 40)$
+ \STATE $  \nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2)) = (x_L , y_L, z_L) $
+\newline
+ \STATE {//Calcule de la direction de la pente dk (méthode de Newton) : }
+ \STATE $ d_k = -H[J](x,y,z)^{-1}* J(x,y,z)$ \hfill $ //résultat : (-(50,50,0))$
+ \newline
+ \STATE {//Calcul nouvelles valeurs des coordonnées}
+ \STATE $ (x_{k+1},y_{k+1},z_{k+1}) = (x_k,y_k,z_k)+ d_k $ \hfill $ //résultat : (50,50,0)$
+ \newline
+ \STATE {//Incrémentation de k}
+ \STATE $ k \leftarrow k+1$\hfill $ //k = 1$
+\newline
+
+ \STATE {//Deuxième itération :}
+ \STATE{//Calcule du gradient de $ J $ :}
+ \STATE $ \nabla J(x,y,z) = (2x_k,2y_k,2z_k)$  \hfill $  // résultat : (100,100,0) $
+\newline
+\STATE {//Calcule du gradient de $ L $ : }
+\STATE $\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2)) = \nabla J(x_k,y_k,z_k) + \lambda_1 \nabla g_1(x_a,y_a,z_a) + \lambda_2 \nabla g_2(x_a,y_a,z_a)) $ \hfill $// résultat : (120, 120, 0)$
+ \STATE $  \varepsilon _1 = \norme{\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2))}$
+\newline
+ \STATE {//Calcule de la direction de la pente dk (méthode de Newton) : }
+ \STATE $ d_k = -H[J](x,y,z)^{-1}* J(x,y,z)$ \hfill $ //résultat : (-(25,25,0))$
+ \STATE {//Calcul nouvelles valeurs des coordonnées}
+ \STATE $ (x_{k+1},y_{k+1},z_{k+1}) = (x_k,y_k,z_k)+ d_k $ \hfill $ //résultat : (25,25,0)$
+ \newline
+ \STATE {//Incrémentation de k}
+\STATE $ k \leftarrow k+1$\hfill $ //k = 2$
+\newline
+
+\STATE {//Troisième itération :}
+\STATE{//Calcule du gradient de $ J $ :}
+\STATE $ \nabla J(x,y,z) = (2x_k,2y_k,2z_k)$  \hfill $  // résultat : (50,50,0) $
+\newline
+\STATE {//Calcule du gradient de $ L $ : }
+\STATE $\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2)) = \nabla J(x_k,y_k,z_k) + \lambda_1 \nabla g_1(x_a,y_a,z_a) + \lambda_2 \nabla g_2(x_a,y_a,z_a)) $ \hfill $// résultat : (70, 70, 0)$
+\STATE $  \varepsilon _1 = \norme{\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2))}$
+\newline
+\STATE {//Calcule de la direction de la pente dk (méthode de Newton) : }
+\STATE $ d_k = -H[J](x,y,z)^{-1}* J(x,y,z)$ \hfill $ //résultat : (-(12.5,12.5,0))$
+\STATE {//Calcul nouvelles valeurs des coordonnées}
+\newline
+\STATE $ (x_{k+1},y_{k+1},z_{k+1}) = (x_k,y_k,z_k)+ d_k $ \hfill $ //résultat : (12.5,12.5,0)$
+\STATE {//Incrémentation de k}
+\STATE $ k \leftarrow k+1$\hfill $ //k = 3$
+\newline
+
+\STATE {//Quatrième itération :}
+\STATE{//Calcule du gradient de $ J $ :}
+\STATE $ \nabla J(x,y,z) = (2x_k,2y_k,2z_k)$  \hfill $  // résultat : (25,25,0) $
+\newline
+\STATE {//Calcule du gradient de $ L $ : }
+\STATE $\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2)) = \nabla J(x_k,y_k,z_k) + \lambda_1 \nabla g_1(x_a,y_a,z_a) + \lambda_2 \nabla g_2(x_a,y_a,z_a)) $ \hfill $// résultat : (45, 45, 0)$
+\STATE $  \varepsilon _1 = \norme{\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2))}$
+\newline
+\STATE {//Calcule de la direction de la pente dk (méthode de Newton) : }
+\STATE $ d_k = -H[J](x,y,z)^{-1}* J(x,y,z)$ \hfill $ //résultat : (-(6.25,6.25,0))$
+\newline
+\STATE {//Calcul nouvelles valeurs des coordonnées}
+\STATE $ (x_{k+1},y_{k+1},z_{k+1}) = (x_k,y_k,z_k)+ d_k $ \hfill $ //résultat : (6.25,6.25,0)$
+\STATE {//Incrémentation de k}
+\newline
+\STATE $ k \leftarrow k+1$\hfill $ //k = 4$
+\STATE $ $
+
+\STATE {//Cinquième itération :}
+\STATE{//Calcule du gradient de $ J $ :}
+\STATE $ \nabla J(x,y,z) = (2x_k,2y_k,2z_k)$  \hfill $  // résultat : (12.5,12.5,0) $
+\newline
+\STATE {//Calcule du gradient de $ L $ : }
+\STATE $\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2)) = \nabla J(x_k,y_k,z_k) + \lambda_1 \nabla g_1(x_a,y_a,z_a) + \lambda_2 \nabla g_2(x_a,y_a,z_a)) $ \hfill $// résultat : (32.5, 32.5, 0)$
+\STATE $  \varepsilon _1 = \norme{\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2))}$
+\newline
+\STATE {//Calcule de la direction de la pente dk (méthode de Newton) : }
+\STATE $ d_k = -H[J](x,y,z)^{-1}* J(x,y,z)$ \hfill $ //résultat : (-(3.125,3.125,0))$
+\newline
+\STATE {//Calcul nouvelles valeurs des coordonnées}
+\STATE $ (x_{k+1},y_{k+1},z_{k+1}) = (x_k,y_k,z_k)+ d_k $ \hfill $ //résultat : (3.125,3.125,0)$
+\newline
+\STATE {//Incrémentation de k}
+\STATE $ k \leftarrow k+1$\hfill $ //k = 5$
+\newline
+
+\STATE {//Sixième itération :}
+\STATE{//Calcule du gradient de $ J $ :}
+\STATE $ \nabla J(x,y,z) = (2x_k,2y_k,2z_k)$  \hfill $  // résultat : (6.25,6.25,0) $
+\newline
+\STATE {//Calcule du gradient de $ L $ : }
+\STATE $\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2)) = \nabla J(x_k,y_k,z_k) + \lambda_1 \nabla g_1(x_a,y_a,z_a) + \lambda_2 \nabla g_2(x_a,y_a,z_a)) $ \hfill $// résultat : (26.25, 26.25, 0)$
+\STATE $  \varepsilon _1 = \norme{\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2))}$
+\newline
+\STATE {//Calcule de la direction de la pente dk (méthode de Newton) : }
+\STATE $ d_k = -H[J](x,y,z)^{-1}* J(x,y,z)$ \hfill $ //résultat : (-(1.5625,1.5625,0))$
+\STATE {//Calcul nouvelles valeurs des coordonnées}
+\newline
+\STATE $ (x_{k+1},y_{k+1},z_{k+1}) = (x_k,y_k,z_k)+ d_k $ \hfill $ //résultat : (1.5625,1.5625,0)$
+\STATE {//Incrémentation de k}
+\newline
+\STATE $ k \leftarrow k+1$\hfill $ //k = 6$
+\newline
+
+\STATE {//Septième itération :}
+\STATE{//Calcule du gradient de $ J $ :}
+\STATE $ \nabla J(x,y,z) = (2x_k,2y_k,2z_k)$  \hfill $  // résultat : (3.125, 3.125, 0) $
+\newline
+\STATE {//Calcule du gradient de $ L $ : }
+\STATE $\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2)) = \nabla J(x_k,y_k,z_k) + \lambda_1 \nabla g_1(x_a,y_a,z_a) + \lambda_2 \nabla g_2(x_a,y_a,z_a)) $ \hfill $// résultat : (23.125, 23.125, 0)$
+\STATE $  \varepsilon _1 = \norme{\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2))}$
+\newline
+\STATE {//Calcule de la direction de la pente dk (méthode de Newton) : }
+\STATE $ d_k = -H[J](x,y,z)^{-1}* J(x,y,z)$ \hfill $ //résultat : (-(0.78125,0.78125,0))$
+\STATE {//Calcul nouvelles valeurs des coordonnées}
+\newline
+\STATE $ (x_{k+1},y_{k+1},z_{k+1}) = (x_k,y_k,z_k)+ d_k $ \hfill $ //résultat : (0.78125,0.78125,0)$
+\STATE {//Incrémentation de k}
+\newline
+\STATE $ k \leftarrow k+1$\hfill $ //k = 7$
+\newline
+
+\STATE {//Huitième itération :}
+\STATE{//Calcule du gradient de $ J $ :}
+\STATE $ \nabla J(x,y,z) = (2x_k,2y_k,2z_k)$  \hfill $  // résultat : (1.5625, 1.5625, 0) $
+\newline
+\STATE {//Calcule du gradient de $ L $ : }
+\STATE $\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2)) = \nabla J(x_k,y_k,z_k) + \lambda_1 \nabla g_1(x_a,y_a,z_a) + \lambda_2 \nabla g_2(x_a,y_a,z_a)) $ \hfill $// résultat : (21.5625, 21.5625, 0)$
+\STATE $  \varepsilon _1 = \norme{\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2))}$
+\newline
+\STATE {//Calcule de la direction de la pente dk (méthode de Newton) : }
+\STATE $ d_k = -H[J](x,y,z)^{-1}* J(x,y,z)$ \hfill $ //résultat : (-(0.390625,0.390625,0))$
+\newline
+\STATE {//Calcul nouvelles valeurs des coordonnées}
+\STATE $ (x_{k+1},y_{k+1},z_{k+1}) = (x_k,y_k,z_k)+ d_k $ \hfill $ //résultat : (0.390625,0.390625,0)$
+\newline
+\STATE {//Incrémentation de k}
+\STATE $ k \leftarrow k+1$\hfill $ //k = 8$
+\newline
+
+\STATE {//neuvième itération :}
+\STATE{//Calcule du gradient de $ J $ :}
+\STATE $ \nabla J(x,y,z) = (2x_k,2y_k,2z_k)$  \hfill $  // résultat : (0.78125, 0.78125, 0) $
+\newline
+\STATE {//Calcule du gradient de $ L $ : }
+\STATE $\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2)) = \nabla J(x_k,y_k,z_k) + \lambda_1 \nabla g_1(x_a,y_a,z_a) + \lambda_2 \nabla g_2(x_a,y_a,z_a)) $ \hfill $// résultat : (20.78125, 20.78125, 0)$
+\STATE $  \varepsilon _1 = \norme{\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2))}$
+\newline
+\STATE {//Calcule de la direction de la pente dk (méthode de Newton) : }
+\STATE $ d_k = -H[J](x,y,z)^{-1}* J(x,y,z)$ \hfill $ //résultat : (-(0.1953125,0.1953125,0))$
+\newline
+\STATE {//Calcul nouvelles valeurs des coordonnées}
+\STATE $ (x_{k+1},y_{k+1},z_{k+1}) = (x_k,y_k,z_k)+ d_k $ \hfill $ //résultat : (0.1953125,0.1953125,0)$
+\newline
+\STATE {//Incrémentation de k}
+\STATE $ k \leftarrow k+1 \hfill  //k = 9$
+\newline
+
+\STATE {//Dixième itération :}
+\STATE{//Calcule du gradient de $ J $ :}
+\STATE $ \nabla J(x,y,z) = (2x_k,2y_k,2z_k)$  \hfill $  // résultat : (0.390625, 0.390625, 0) $
+\newline
+\STATE {//Calcule du gradient de $ L $ : }
+\STATE $\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2)) = \nabla J(x_k,y_k,z_k) + \lambda_1 \nabla g_1(x_a,y_a,z_a) + \lambda_2 \nabla g_2(x_a,y_a,z_a)) $ \hfill $// résultat : (20.390625, 20.390625, 0)$
+\STATE $  \varepsilon _1 = \norme{\nabla L((x_k,y_k,z_k),(\lambda_1,\lambda_2))}$
+\newline
+\STATE {//Calcule de la direction de la pente dk (méthode de Newton) : }
+\STATE $ d_k = -H[J](x,y,z)^{-1}* J(x,y,z)$ \hfill $ //résultat : (-(0.097665625,0.097665625,0))$
+\newline
+\STATE {//Calcul nouvelles valeurs des coordonnées}
+\STATE $ (x_{k+1},y_{k+1},z_{k+1}) = (x_k,y_k,z_k)+ d_k $ \hfill $ //résultat : (0.097665625,0.097665625,0)$
+\newline
+\STATE {//Incrémentation de k}
+\STATE $ k \leftarrow k+1$\hfill $ //k = 10$
+\newline
+\STATE {// Fin de la boucle "while" car nous avons atteint k =10, condition mettant fin à la //boucle}
+\newline
+
+ \ENDWHILE
+
+\end{algorithmic}
+\end{algorithmfloat}
+
+
+\hrulefill
+
 \bibliographystyle{plain}
 \bibliography{stdlib_sbphilo}