+Considérons un problème d’optimisation différentiable $ \mathcal{P} $ avec contraintes d’égalité :
+$$
+ \mathcal{P} \left \{
+ \begin{array}{l}
+ \displaystyle\min_{x \in \mathbb{R}^n} J(x) \\
+ h(x) = 0
+ \end{array}
+ \right .
+$$
+où $ J: \mathbb{R}^n \longrightarrow \mathbb{R} $ et $h: \mathbb{R}^n \longrightarrow \mathbb{R}^q$ sont supposées au moins différentiables.
+\newline
+Les conditions d’optimalité de Lagrange (ou \textit{KKT}) s’écrivent :
+$$ \nabla J(x) + \sum\limits_{i=1}^{q} \lambda_i \nabla h_i(x) = 0 \iff \nabla L(x,\lambda) = 0 $$
+donc $ \mathcal{P} $ devient :
+$$ \begin{pmatrix}
+ \nabla J(x) + \sum\limits_{i=1}^{q} \lambda_i \nabla h_i(x) \\
+ h(x)
+ \end{pmatrix} = 0 $$
+Pour résoudre ce système d’équations, utilisons la méthode de Newton dont une itération s’écrit ici :
+$$ H[L](x_k,\lambda_k)\begin{pmatrix}
+ x_{k+1} - x_k \\
+ \lambda_{k+1} - \lambda_k
+ \end{pmatrix} = -\nabla L(x_k,\lambda_k) $$
+soit :
+$$ \begin{pmatrix}
+ H_x[L](x_k,\lambda_k) & D_h(x_k)^\top \\
+ D_h(x_k) & 0
+ \end{pmatrix} \begin{pmatrix}
+ x_{k+1} - x_k \\
+ \lambda_{k+1} - \lambda_k
+ \end{pmatrix} = -\begin{pmatrix}
+ \nabla_x L(x_k,\lambda_k) \\
+ 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)^\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 \\
+ D_h(x_k) & 0
+ \end{pmatrix} \begin{pmatrix}
+ d \\
+ \mu - \lambda_k
+ \end{pmatrix} = -\begin{pmatrix}
+ \nabla_x L(x_k,\lambda_k) \\
+ h(x_k)
+ \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)^\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 :
+$$
+ \left \{
+ \begin{array}{r c l}
+ H_kd + \sum\limits_{i=1}^q(\mu_i - \lambda_{k_i})\nabla h_i(x_k) & = & -\nabla_x L(x_k,\lambda_k) \\
+ \nabla h_i(x_k)^\top d + h_i(x_k) & = & 0, \ \forall i \in \{1,\ldots,q\}
+ \end{array}
+ \right .
+$$
+Or $ \nabla_x L(x_k,\lambda_k) = \nabla J(x_k) + \sum\limits_{i=1}^{q} \lambda_{k_i} \nabla h_i(x_k) $, d'où :
+$$
+ \left \{
+ \begin{array}{r c l}
+ H_kd + \sum\limits_{i=1}^q\mu_i\nabla h_i(x_k) & = & -\nabla J(x_k) \\
+ \nabla h_i(x_k)^\top d + h_i(x_k) & = & 0, \ \forall i \in \{1,\ldots,q\}
+ \end{array}
+ \right .
+$$
+On reconnait dans le système ci-dessus les conditions d’optimalité de Lagrange du problème quadratique suivant :
+$$
+ \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 \\
+ h_i(x_k) + \nabla h_i(x_k)^\top d = 0, \ \forall i \in \{1,\ldots,q\}
+ \end{array}
+ \right .
+$$
+Le problème $ \mathcal{PQ}_k $ peut être vu comme la minimisation d’une approximation quadratique du Lagrangien de $ \mathcal{P} $ avec une approximation linéaire des contraintes.
+\newline
+Comme son nom l’indique, la méthode PQS consiste à remplacer le problème initial par une suite de problèmes quadratiques sous contraintes linéaires plus faciles à résoudre. L’algorithme est le suivant :
+
+\hrulefill
+\newline
+ALGORITHME PQS AVEC CONSTRAINTES D'ÉGALITÉ.
+\newline
+\textit{Entrées}: $ J : \mathbb{R}^n \longrightarrow \mathbb{R} $, $ h : \mathbb{R}^n \longrightarrow \mathbb{R}^q $ différentiables, $ x_0 \in \mathbb{R}^n $ point initial arbitraire, $ \lambda_0 \in \mathbb{R}^q $ multiplicateur initial, $ \varepsilon > 0 $ précision demandée.
+\newline
+\textit{Sortie}: une approximation $ x_k $ de la solution $ x^\ast $ du problème $ \mathcal{P} $.
+\begin{enumerate}
+ \item $ k := 0 $.
+ \item Tant que $ \norme{\nabla L(x_k,\lambda_k)} > \varepsilon $,
+ \begin{enumerate}
+ \item Résoudre le sous-problème quadratique :
+ $$
+ \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 \\
+ h_i(x_k) + \nabla h_i(x_k)^\top d = 0, \ \forall i \in \{1,\ldots,q\}
+ \end{array}
+ \right .
+ $$
+ et obtenir la solution primale $ d_k $ et le multiplicateur $ \lambda^{\prime} $ associé à la contrainte d’égalité.
+ \item $ x_{k+1} = x_k + d_k; \ \lambda_{k+1} = \lambda^{\prime}; \ k := k + 1 $.
+ \end{enumerate}
+ \item Retourner $ x_k $.
+\end{enumerate}
+
+\hrulefill
+
+\subsubsection{Contraintes d’inégalité}
+
+Intéressons nous maintenant aux problèmes avec contraintes d’égalité et d’inégalité :
+$$
+ \mathcal{P} \left \{
+ \begin{array}{l}
+ \displaystyle\min_{x \in \mathbb{R}^n} J(x) \\
+ g(x) \leq 0 \\
+ h(x) = 0
+ \end{array}
+ \right .
+$$
+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 à 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 $ :
+$$
+ \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(x_k) + \nabla g(x_k)^\top d \leq 0 \\
+ h(x_k) + \nabla h(x_k)^\top d = 0
+ \end{array}
+ \right .
+$$
+Comme pour la cas précédent, on résoud alors le sous problème quadratique $ \mathcal{PQ}_k $ qui est plus simple avec une méthode vue plus haut en posant $ \mathcal{Q} = H_k = H[L](x_k,\lambda_k,\mu_k) $.
+\newpage
+\hrulefill
+\newline
+ALGORITHME PQS AVEC CONSTRAINTES D'ÉGALITÉ ET D'INÉGALITÉ.
+\newline
+\textit{Entrées}: $ J : \mathbb{R}^n \longrightarrow \mathbb{R} $, $g: \mathbb{R}^n \longrightarrow \mathbb{R}^p$, $ h : \mathbb{R}^n \longrightarrow \mathbb{R}^q $ différentiables, $ x_0 \in \mathbb{R}^n $ point initial arbitraire, $ \lambda_0 \in \mathbb{R}_+^p $ et $ \mu_0 \in \mathbb{R}_+^q $ multiplicateurs initiaux, $ \varepsilon > 0 $ précision demandée.
+\newline
+\textit{Sortie}: une approximation $ x_k $ de la solution $ x^\ast $ du problème $ \mathcal{P} $.
+\begin{enumerate}
+ \item $ k := 0 $.
+ \item Tant que $ \norme{\nabla L(x_k,\lambda_k,\mu_k)} > \varepsilon $,
+ \begin{enumerate}
+ \item Résoudre le sous-problème quadratique :
+ $$
+ \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 \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 .
+ $$
+ et obtenir la solution primale $ d_k $ et les multiplicateurs $ \lambda^{\prime} $ et $ \mu^{\prime} $ associé aux contraintes d’inégalité et d’égalité respectivement.
+ \item $ x_{k+1} = x_k + d_k; \ \lambda_{k+1} = \lambda^{\prime}; \ \mu_{k+1} = \mu^{\prime}; \ k := k + 1 $.
+ \end{enumerate}
+ \item Retourner $ x_k $.
+\end{enumerate}
+
+\hrulefill
+\newline
+Afin que le sous-programme quadratique $ \mathcal{PQ}_k $ admette une unique solution, la plupart des implémentations actuelles de PQS utilisent une approximation du hessien $ H_k $ du Lagrangien qui soit définie positive, en particulier celle fournie par les techniques quasi-newtonienne (BFGS) par exemple.
+\newline
+Étant 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 $$
+Les méthodes de mises à jour PSB et BFGS suivent par exemple cette stratégie.
+
+\subsubsection{Mises à jour PSB}
+
+$$ H_{k+1} = H_k + \frac{1}{d^\top d} ((y - H_k d)d^\top + d(y - H_k d)^\top) - \frac{(y - H_k d)^\top d}{(d^\top d)^2}d^\top d $$
+où
+$$ d = x_{k+1} - x_k $$
+$$ y = \nabla L(x_{k+1},\lambda_{k+1},\mu_{k+1}) - \nabla L(x_{k},\lambda_{k+1},\mu_{k+1}) $$
+
+\subsubsection{Mises à jour BFGS}
+
+$$ H_{k+1} = H_k + \frac{y y^\top}{y^\top d} - \frac{H_k d d^\top H_k}{d^\top H_k d}$$
+où
+$$ d = x_{k+1} - x_k $$
+$$ y = \nabla L(x_{k+1},\lambda_{k+1},\mu_{k+1}) - \nabla L(x_{k},\lambda_{k+1},\mu_{k+1}) $$
+
+\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} \land r < r_1 \land r < r_2. $$
+\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)) $$
+$$ = (2x(1 + \lambda_1 + \lambda_2),2y(1 + \lambda_1),2z(1 + \lambda_2)) $$
+\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} $.