Fix the SQP algorithm.
[Projet_Recherche_Operationnelle.git] / rapport / ProjetOptimRO.tex
CommitLineData
e37aec65
JB
1\documentclass[12pt,oneside,a4paper]{book}
2
3
4%%%%%Packages
5
de30386e 6
e37aec65 7\usepackage{latexsym}
66a4e4ad 8\usepackage{amsmath}
8000c039 9\usepackage{amsthm}
66a4e4ad 10\usepackage{mathtools}
e37aec65
JB
11\usepackage{amssymb}
12\usepackage[utf8]{inputenc}
13\usepackage[francais]{babel}
14\usepackage{color}
15\usepackage{geometry}
16\usepackage{graphicx}
17\usepackage{amsfonts}
18\usepackage[T1]{fontenc}
19\usepackage{multirow}
20\usepackage{fancyhdr}
21\usepackage{tocbibind}
22\usepackage{lmodern}
7fd51f8b 23\usepackage{enumitem}
e37aec65
JB
24
25
26%%%%%Marges & en-t\^etes
27
28\geometry{hmargin=2.3cm, vmargin=3cm}
29\fancyhf{} % supprime les en-t\^etes et pieds pr\'ed\'efinis
30\fancyhead[FC]{\bfseries\thepage} % N∞page centre bas
31\fancyhead[HC]{\footnotesize\leftmark} % chapitre centre haut
32\renewcommand{\headrulewidth}{0.2pt} % filet en haut
33\addtolength{\headheight}{0.5pt} % espace pour le filet
66a4e4ad 34\renewcommand{\footrulewidth}{0.2pt} % filet en bas
e37aec65
JB
35
36
37%%%%%Th\'eor\`eme et d\'efinitions
38
39\newtheorem{Def}{D\'efinition}
40\newtheorem{Not}[Def]{Notation}
41\newtheorem{Th}{Th\'eor\`eme}
42\newtheorem{Prop}[Th]{Proposition}
43\newtheorem{Cor}[Th]{Corollaire}
44\newtheorem{Rmq}{Remarque}
45
682e0379 46\newcommand{\norme}[1]{\left\Vert #1 \right\Vert}
e37aec65
JB
47
48%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
50
51\begin{document}
52
53%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
54%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
55
56%%%%%Page de garde
57
58\begin{center}
59
b0927d0a 60 %\includegraphics[scale=0.5]{logo_sciences_rvb.png}\\
cec1e8f8 61 \includegraphics[scale=0.5]{polytech.png}\\
b0927d0a
JB
62
63 \vspace*{0.5cm}
64
65 \footnotesize{
66 \large \bf D\'epartement d'Informatique, Réseaux et Multimédia\\
67 \large \bf 5ème année\\
68 }
69
70 \vspace*{0.5cm}
71
72 %\large{Master 2 Professionnel\\
73 %Math\'ematiques et Informatique des Nouvelles Technologies\\}
74
75 \large{Projet \\ en \\ Optimisation et Recherche Opérationnelle \\}
76
77 \vspace*{0.7cm}
78
79 \begin{tabular}{c}
80 \hline
9f054a9c 81 ~ \\
7fd51f8b 82 \LARGE\textbf {Programmation Quadratique Séquentielle ou PQS} \\
9f054a9c
JB
83 \LARGE\textbf {en} \\
84 \LARGE\textbf {Optimisation non linéraire sous contraintes} \\
85 ~ \\
b0927d0a
JB
86 \hline
87 \end{tabular}
88
89 \vspace*{0.7cm}
90
91 \includegraphics[scale=0.4]{CE.PNG}\\
92
93 \vspace*{0.5cm}
94
95 \large par\\
96
97 %\large \bsc{}\\
98 %\normalsize{M\'emoire encadr\'e par :} \large St\'ephane \bsc{Ballet}\\
99
100 \vspace*{0.2cm}
91df3de1 101 \large {\bf Jérôme \bsc{Benoit} et Sylvain \bsc{Papa}}\\
b0927d0a
JB
102
103 %\vspace*{0.1cm}
104
105 % \large sous la direction de \\
106
107 %\vspace*{0.1cm}
108
109 %Eric Audureau et Thierry Masson
110
111 %\vspace*{1cm}
112
113 \vspace*{1cm}
114
115 %\normalsize{Licence de Mathématiques 3ème année}
116 \normalsize{Année 2018-2019}
e37aec65
JB
117
118\end{center}
119
120\thispagestyle{empty}
121
122\newpage
123
124
125%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
126%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
127
128
129\pagestyle{plain}
130\frontmatter
131
132
133%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
134%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
135
136
e37aec65
JB
137%%%%%Table des mati\`eres
138
139\tableofcontents
140
141\begin{figure}[!b]
b0927d0a
JB
142 \begin{center}
143 %\includegraphics{logo_fac2}
144 \includegraphics[scale=0.04]{amu}
145 \end{center}
e37aec65
JB
146\end{figure}
147
148\newpage
149
150
151%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
152%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
153
154
155\mainmatter
156\pagestyle{fancy}
157
158
159%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
160\chapter{Introduction générale}
161
8000c039
JB
162L'objectif de ce chapitre est de faire un bref rappel des définitions, notions et résultats essentiels en recherche opérationnelle ainsi que en mathématiques nécessaires à l'étude de la méthode PQS.
163\newline
164Elle est loin d'être exhaustive mais devrait suffire dans le cadre de ce projet.
165
e37aec65
JB
166\vspace{.5em}
167
168\section{Qu'est-ce que la recherche opérationnelle?}
169
64f7c064
JB
170\subsection{Présentation rapide}
171
f899c72b 172La recherche opérationnelle est une discipline dite "hybride" au confluent de plusieurs disciplines dont principalement les mathématiques (l'analyse numérique, les probabilités, la statistique) et l'informatique (l'algorithmie).
e7e85554 173\newline
f899c72b 174On la considère usuellement comme une sous discipline des mathématiques de la décision. Elle a de nombreuses applications, particulièrement en intelligence artificielle.
91df3de1 175
64f7c064
JB
176\subsection{Définition de la problèmatique}
177
f899c72b 178Définissons le problème central $ \mathcal{P} $ que se propose de résoudre la recherche opérationnelle :
6ec0df37 179\begin{Def}
5e4341d1
JB
180 Soient $(n, p, q) \in \mathbb{N}^3$, $x \in \mathbb{R}^n$, une fonction $g: \mathbb{R}^n \longrightarrow \mathbb{R}^p$ représentant les contraintes d'inégalités, une fonction $h: \mathbb{R}^n \longrightarrow \mathbb{R}^q$ représentant les contraintes d'égalités et une fonction dite objectif $J: \mathbb{R}^n \longrightarrow \mathbb{R}$.
181 \newline
182 La problèmatique $ \mathcal{P} $ se définit par :
183 $$
aa023e1c 184 \mathcal{P} \left \{
7fd51f8b 185 \begin{array}{l}
aa023e1c
JB
186 \displaystyle\min_{x \in \mathbb{R}^n} J(x) \\
187 g(x) \leq 0 \\
188 h(x) = 0
189 \end{array}
190 \right .
5e4341d1 191 $$
6ec0df37 192\end{Def}
6ec0df37 193\begin{Def}
5e4341d1
JB
194 On définit $ \mathcal{C} $ l'ensemble des contraintes par :
195 $$ \mathcal{C} = \left \{ x \in \mathbb{R}^n \ | \ g(x) \leq 0 \land h(x) = 0 \right \} $$
6ec0df37 196\end{Def}
3b344e8c 197Elle se doit de résoudre les problèmes d'existence d'une solution ($ \mathcal{C} \neq \emptyset $ et $ \displaystyle\min_{x \in \mathbb{R}^n} J(x) $ défini dans $ \mathcal{C} $) ainsi que de construction d'une solution dans $ \mathcal{C} $.
64f7c064 198
e37aec65
JB
199\section{Qu'est-ce que l'optimisation?}
200
f899c72b
JB
201\subsection{Définition}
202
203La recherche d'une méthode permettant de trouver la solution au problème $ \mathcal{P} $ dans $ \mathcal{C} $ est l'activité principale de l'optimisation.
204\newline
8a00a107 205Si la modélisation de la problèmatique $ \mathcal{P} $ est considérée comme un art, la recherche d'une solution au problème $ \mathcal{P} $ dans $ \mathcal{C} $ est, elle, une science.
f899c72b
JB
206
207\subsection{Quelques définitions annexes}
208
209Définissons quelques notions supplémentaires de base nécessaires à la suite :
682e0379
JB
210\begin{Def}
211 On définit le Lagrangien associé à $ \mathcal{P} $ par :
212 $$ \begin{array}{r c l}
213 L : \mathbb{R}^n \times \mathbb{R}^q \times \mathbb{R}_+^p & \longrightarrow & \mathbb{R} \\
7fd51f8b 214 (x,\lambda,\mu) & \longmapsto & L(x,\lambda,\mu) = J(x) + \sum\limits_{i=1}^{q} \lambda_i h_i(x) + \sum\limits_{j=1}^{p} \mu_j g_j(x) \\
682e0379
JB
215 & & L(x,\lambda,\mu) = J(x) + \langle \lambda,h(x) \rangle_{\mathbb{R}^q} + \langle \mu,g(x) \rangle_{\mathbb{R}^p}
216 \end{array} $$
217 où l’on note $ \lambda $ et $ \mu $ les vecteurs de coordonnées respectives $ (\lambda_1,\ldots,\lambda_q) $ et $ (\mu_1,\ldots,\mu_p) $.
218\end{Def}
f899c72b 219\begin{Def}
9f054a9c
JB
220 Soient $ \mathbb{R}^n $ un espace topologique, $ A \subset \mathbb{R}^n $ et $ x^\ast \in \mathbb{R}^n $.
221 \newline
222 On dit que $ x^\ast $ est \textbf{intérieur} à $ A $ si $ A $ est un voisinage de $ x^\ast $. On appelle intérieur de $ A $ l'ensemble des points intérieurs à $ A $ et on le note $ \mathring{A} $.
f899c72b 223\end{Def}
9f054a9c
JB
224\begin{Rmq}
225 $ A \subset \mathbb{R}^n $ est un ouvert $ \iff A = \mathring{A} $.
226\end{Rmq}
f899c72b 227\begin{Def}
9f054a9c
JB
228 Soient $ \mathbb{R}^n $ un espace topologique, $ A \subset \mathbb{R}^n $ et $ x^\ast \in \mathbb{R}^n $.
229 \newline
230 On dit que $ x^\ast $ est \textbf{adhérent} à $ A $ si et seulement si $ \forall V \in \mathcal{V}(x^\ast) \ A \cap V \neq \emptyset $. On appelle adhérence de $ A $ l'ensemble des points adhérents à $ A $ et on le note $ \overline{A} $.
f899c72b 231\end{Def}
9f054a9c
JB
232\begin{Rmq}
233 $ A \subset \mathbb{R}^n $ est un fermé $ \iff A = \overline{A} $.
234\end{Rmq}
f899c72b 235\begin{Def}
7590f4bb
JB
236 Soient $ f : \mathbb{R}^n \longrightarrow \mathbb{R} $ et $ S \subset \mathbb{R}^n $. On définit $ \mathrm{argmin} $ de $ f $ sur $ S $ par :
237 $$ \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) \} $$
238\end{Def}
239\begin{Def}
240 Soient une fonction $ f : \mathbb{R}^n \longrightarrow \mathbb{R} $ et $ x^\ast \in \mathbb{R}^n $.
9f054a9c
JB
241 \newline
242 On dit que $ f $ est continue en $ x^\ast $ si
243 $$ \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 $$
f899c72b 244\end{Def}
d17ef079
JB
245\begin{Def}
246 Soient $ k \in \{ 1,\ldots,n \} $ et une fonction $ f: \mathbb{R}^n \longrightarrow \mathbb{R} $.
247 \newline
248 On dit que la $ k^{ième} $ dérivée partielle de $ f $ existe au point $ x^\ast \in \mathbb{R}^n $ si l’application
249 $$ t \longmapsto f(x^\ast_1,\ldots,x^\ast_{k-1},x^\ast_k + t,x^\ast_{k+1},\ldots,x^\ast_n) $$
250 définie sur un voisinage de $ 0 $ dans $ \mathbb{R} $ et à valeurs dans $ \mathbb{R} $ est dérivable en $ 0 $.
251 \newline
252 Dans ce cas on note
253 $$ \frac{\partial f}{\partial x_k}(x^\ast) $$ ou $$ \partial_k f(x^\ast) $$
254 cette dérivée.
255\end{Def}
66a4e4ad 256\begin{Def}
7590f4bb 257 Soient une fonction $ f : \mathbb{R}^n \longrightarrow \mathbb{R} $
66a4e4ad 258 et $ x^\ast, h \in \mathbb{R}^n $.
d17ef079 259 \newline
66a4e4ad
JB
260 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
261 \[
262 f(x^\ast + h) = f(x^\ast) + d_{x^\ast}f(h) + \underset{h \rightarrow 0}{\mathrm{o}}(\norme{h})
263 \]
264 Autrement dit il existe une application $ \varepsilon_{x^\ast} $ définie sur le voisinage de $ 0 $ dans $ \mathbb{R}^n $ et à valeurs dans $ \mathbb{R} $
265 telle que $ \lim\limits_{h \rightarrow 0} \varepsilon_{x^\ast}(h) = 0 $ et
266 \[
267 f(x^\ast + h) = f(x^\ast) + d_{x^\ast}f(h) + \norme{h}\varepsilon_{x^\ast}(h)
268 \]
d17ef079 269 On appelle $ d_{x^\ast}f $ la différentielle de $ f $ en $ x^\ast $.
66a4e4ad 270\end{Def}
d17ef079
JB
271\begin{Rmq}
272 On peut démontrer que : $$ d_{x^\ast}f = \sum_{i=1}^n\frac{\partial f}{\partial x_i}(x^\ast) $$.
273\end{Rmq}
5e4341d1
JB
274\begin{Def}
275 Soit une fonction $ f: \mathbb{R}^n \longrightarrow \mathbb{R} $ différentiable.
276 \newline
277 Le gradient de $ f $, noté $\nabla f$, en $ x^\ast \in \mathbb{R}^n$ se définit par :
278 \[
de30386e 279 \nabla f(x^\ast) = (\frac{\partial f}{\partial x_1}(x^\ast),\ldots,\frac{\partial f}{\partial x_n}(x^\ast))
5e4341d1
JB
280 \]
281\end{Def}
d17ef079 282\begin{Rmq}
a5f09ac4 283 $ \forall h \in \mathbb{R}^n \ d_{x^\ast}f(h) = \langle \nabla f(x^\ast),h \rangle = \nabla f(x^\ast)^\top h $
d17ef079 284\end{Rmq}
9f054a9c 285\begin{Def}
84af1fe2 286 Soit $ f: \mathbb{R}^n \longrightarrow \mathbb{R} $ une fonction de classe $ \mathcal{C}^2 $.
9f054a9c
JB
287 On définit la matrice hessienne de $ f $ en $ x^\ast $ par :
288 $$ H[f](x^\ast) =
289 \begin{pmatrix}
290 \frac{\partial^2 f}{\partial x_1^2}(x^\ast) & \cdots & \frac{\partial^2 f}{\partial x_1\partial x_n}(x^\ast) \\
291 \vdots & & \vdots \\
292 \frac{\partial^2 f}{\partial x_n\partial x_1}(x^\ast) & \cdots & \frac{\partial^2 f}{\partial x_n^2}(x^\ast)
293 \end{pmatrix} $$
294\end{Def}
295\begin{Prop}
296 \begin{enumerate}
297 \item $ H[f](x^\ast) $ est une matrice symétrique (Théorème de symétrie de Schwarz).
682e0379 298 \item On a le développement de Taylor-Young à l'ordre 2 en $ x^\ast $ suivant :
9f054a9c 299 $$ f(x^\ast + v) = f(x^\ast) + \langle \nabla f(x^\ast),v \rangle + \frac{1}{2} v^\top H[f](x^\ast) v + \varepsilon(v) $$
682e0379
JB
300 ou
301 $$ f(x^\ast + v) = f(x^\ast) + \langle \nabla f(x^\ast),v \rangle + \frac{1}{2} \langle H[f](x^\ast)v,v \rangle + \varepsilon(v) $$
9f054a9c
JB
302 avec $ \frac{|\varepsilon(v)|}{\norme{v}} \rightarrow 0 $ quand $ \norme{v} \rightarrow 0 $.
303 \end{enumerate}
304\end{Prop}
682e0379
JB
305\begin{proof}
306 Elle repose entièrement sur deux autres théorèmes dont les preuves sont connues et de la réécriture de formulation de résultat.
307\end{proof}
f899c72b
JB
308
309\subsection{Conditions d'existence d'un extremum}
310
7590f4bb
JB
311On peut démontrer que $ \mathcal{C}$ est un ensemble fermé de $ \mathbb{R}^n $ si $ g $ et $ h $ sont continues.
312On peut en déduire $ \mathcal{C} $ est un ensemble fermé et borné de $ \mathbb{R}^n $.
3b344e8c 313\begin{Th}[Théorème de Weierstrass]
9f054a9c
JB
314 Soient $ \mathcal{C} \neq \emptyset \subset \mathbb{R}^n $ un fermé borné et $ f : \mathcal{C} \longrightarrow \mathbb{R} $ une fonction continue.
315 \newline
316 Alors : $$ \exists x^\ast \in \mathcal{C} \ \forall x \in \mathcal{C} \ f(x) \geq f(x^\ast) $$
317 Autrement dit $ x^\ast $ est un minimum global de $ J $ sur $ \mathcal{C} $.
318 \newline
319 De la même façon, il existe un maximum global de $ J $ sur $ \mathcal{C} $.
3b344e8c 320\end{Th}
7590f4bb 321Si $ 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}.
b6cd6632 322
3b344e8c
JB
323\subsection{Conditions de caractérisation d'un extremum}
324
9f054a9c 325Dans le cas où $ J, g, h $ sont continûment différentiable et ses dérivées sont continues (c'est à dire de classe $ \mathcal{C}^1 $), la recherche du mimimum consiste à faire des descentes par gradient [section \ref{descente}] de $ J $ sur $ \mathcal{C} $ avec comme critère d'arrêt : $ x_i = \displaystyle\min_{x \in \mathcal{C}} J(x) \iff \forall \varepsilon \in \mathbb{R}_{+}^{*} \ \norme{\nabla J(x_i)} < \varepsilon $, $ i \in \mathbb{N} $ \cite{FEA}.
f899c72b 326\newline
9f054a9c 327On 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}.
f899c72b 328
7590f4bb 329\subsubsection{Conditions nécessaires de Karuch-Kuhn-Tucker ou \textit{KKT}}\label{KKT}
f899c72b
JB
330
331\begin{Th}
9f054a9c
JB
332 Soient $ x^\ast \in \mathbb{R}^n $, $ I = \{ 1,\ldots,p \} $ et $ J = \{ 1,\ldots,q \} $.
333 \newline
334 Les conditions nécessaires pour que $ x^\ast \in \mathcal{C}$ soit un minimum local de $ J $ sont :
335 \newline
336 \newline
337 \centerline{$ \{ \nabla g_1(x^\ast),\ldots,\nabla g_p(x^\ast),\nabla h_1(x^\ast),\ldots,\nabla h_q(x^\ast) \} $ sont linéairement indépendants.}
338 \newline
339 \newline
340 et
341 $$ \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 $$
342 On appelle $ (\mu_i)_{i \in I}$ les multiplicateurs de Kuhn-Tucker et $ (\lambda_j)_{j \in J}$ les multiplicateurs de Lagrange.
7590f4bb
JB
343 \newline
344 On nomme également les conditions \textit{KTT} conditions nécessaires d'optimalité du premier ordre.
f899c72b 345\end{Th}
8000c039 346\begin{proof}
9f054a9c 347 Elle repose sur le lemme de Farkas \cite{FEA,RON}.
8000c039 348\end{proof}
6079168f 349Il 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.
7590f4bb
JB
350\begin{Def}
351 On appelle un point admissible $ x^\ast \in \mathcal{C} $ un point critique de $ \mathcal{P} $ si il statisfait les conditions \textit{KKT}.
352\end{Def}
8000c039 353
7590f4bb
JB
354\subsubsection{Conditions suffisantes du deuxième ordre}
355
356\begin{Th}
357 Les conditions suffisantes en plus de celles \textit{KKT} pour que $ x^\ast \in \mathcal{C} $ soit un minimum local de $ J $ sont :
358 \begin{enumerate}[label=(\roman*)]
359 \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 $.
360 \item $ \forall v \in \mathbb{R}^n \land v \neq 0 \ \langle H_x[L](x^\ast,\lambda,\mu)v,v \rangle > 0 $.
361 \end{enumerate}
362\end{Th}
8000c039
JB
363
364%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
365
366\chapter{Méthode de programmation quadratique séquentielle}
367
368Dans ce projet, nous nous proposons d'étudier une des méthodes d'optimisation non linéaire avec contraintes nommée programmation quadratique séquentielle ou PQS.
369
370\vspace{.5em}
371
372\section{Methode de descente}\label{descente}
373
682e0379
JB
374Nous supposons que le domaine des contraintes de $ \mathcal{P} $ est un ouvert de $ \mathbb{R}^n $ (c'est à dire que nous n'avons pas de contraintes) et $ J $ est une fonction définie sur $ \mathbb{R}^n $ à valeurs réelles supposée différentiable, voire même deux fois différentiable. Les conditions nécessaires d’optimalité du premier et du second ordre expriment le fait qu’il n’est pas possible de “descendre” à partir d’un point de minimum (local ou global). Cette observation va servir de point de départ à l’élaboration des méthodes dites de descente.
375
8000c039 376Partant d’un point $ x_0 \in \mathbb{R}^n $ arbitrairement choisi, un algorithme de descente va chercher à générer une suite d’itérés $ (x_k)_{k \in \mathbb{N}} $ de $ \mathbb{R}^n $ définie par :
9f054a9c 377$$ x_{k+1} = x_k + s_kd_k $$ où $ s_k \in \mathbb{R}_{+}^{*},d_k \in \mathbb{R}^n $ et avec
8000c039
JB
378$$ \forall k \in \mathbb{N} \ J(x_{k+1}) \leq J(x_k) $$
379Un tel algorithme est ainsi déterminé par deux éléments à chaque étape $ k $ : le choix de la direction $ d_k $ appelée direction de descente, et le choix de la taille du pas $ s_k $ à faire dans la direction $ d_k $. Cette étape est appelée \textit{recherche linéaire}.
380
381\subsection{Définition d'une direction de descente}
382
9f054a9c 383Un vecteur $ d \in \mathbb{R}^n $ est une direction de descente pour $ J $ à partir d’un point $ x_0 \in \mathbb{R}^n $ si $ t \longmapsto f(x_0 + td) $ est strictement décroissante en $ t = 0 $, c’est-à-dire :
8000c039 384$$ \exists \eta \in \mathbb{R}_{+}^{*} \ \forall t \in ]0,\eta] \ J(x_0 + td) < J(x_0) $$
6079168f 385Il est donc important d’analyser le comportement de la fonction $ J $ dans certaines directions.
8000c039 386\begin{Prop}
9f054a9c
JB
387 Soient $ J : \mathbb{R}^n \longrightarrow \mathbb{R} $ différentiable et $ d \in \mathbb{R}^n $.
388 \newline
389 d est un vecteur de descente de $ J $ en $ x_0 \in \mathbb{R}^n $ ssi :
390 $$ \nabla J(x_0)^\top d < 0 $$
391 De plus
392 $$ \forall \beta < 1 \in \mathbb{R}_{+} \ \exists \eta \in \mathbb{R}_{+}^{*} \ \forall t \in ]0,\eta] \ J(x_0 + td) < J(x_0) + t\beta \nabla J(x_0)^\top d < J(x_0) $$
8000c039
JB
393\end{Prop}
394\begin{proof}
9f054a9c 395 Elle s'effectue en utilisant le développement de Taylor-Young de l’application $ t \longmapsto f(x_0 + td) $ à l’ordre 1.
8000c039 396\end{proof}
6079168f 397Cette dernière inégalité garantit une décroissance minimum de la fonction $ J $ dans la direction $ d $ et peut se traduire par : la décroissance de la fonction $ J $, en effectuant un pas de longueur $ t $ dans la direction $ d $, est au moins égale à la longueur du pas multipliée par une fraction de la pente. Le schéma général d’un algorithme de descente est alors le suivant :
8000c039
JB
398
399\hrulefill
400\newline
401ALGORITHME DE DESCENTE MODÈLE.
402\newline
403\textit{Entrées}: $ J : \mathbb{R}^n \longrightarrow \mathbb{R} $ différentiable, $ x_0 \in \mathbb{R}^n $ point initial arbitraire.
404\newline
e17152e0 405\textit{Sortie}: une approximation $ x_k $ de la solution $ x^\ast $ du problème : $ \displaystyle\min_{x \in \mathbb{R}^n} J(x) $.
8000c039 406\begin{enumerate}
e17152e0 407 \item $ k := 0 $.
8000c039 408 \item Tant que "test d’arrêt" non satisfait,
9f054a9c
JB
409 \begin{enumerate}
410 \item Trouver une direction de descente $ d_k $ telle que : $ \nabla J(x_k)^\top d_k < 0 $.
411 \item \textit{Recherche linéaire} : Choisir un pas $ s_k > 0 $ à faire dans cette direction et tel que : $$ J(x_k + s_kd_k) < J(x_k). $$
412 \item Mise à jour : $ x_{k+1} = x_k + s_kd_k; \ k := k + 1 $.
413 \end{enumerate}
414 \item Retourner $ x_k $.
8000c039
JB
415\end{enumerate}
416
417\hrulefill
418
419\subsection{Choix de la direction de descente}
420
e17152e0 421Une fois la théorie bien maîtrisée, calculer une direction de descente est relativement simple. Dans le cas différentiable, il existe deux grandes stratégies de choix de direction de descente :
8000c039
JB
422\begin{itemize}
423 \item la stratégie de Cauchy : $ d_k = -\nabla J(x_k) $, conduisant aux \textit{algorithmes de gradient}.
424 \item la stratégie de Newton : $ d_k = -H[J](x_k)^{-1} \nabla J(x_k) $, conduisant aux \textit{algorithmes Newtoniens}.
425\end{itemize}
9f054a9c 426Remarquons que si $ x_k $ est un point stationnaire ($ \iff \nabla J(x_k) = 0 $) non optimal alors toutes ces directions sont nulles et aucun de ces algorithmes ne pourra progresser. Ce problème peut être résolu en utilisant des approches de type région de confiance qui ne seront pas étudiées dans le cadre de ce projet.
8000c039
JB
427
428\subsection{Critère d’arrêt}
429
9f054a9c 430Soit $ x^\ast $ un minimum local de l'objectif $ J $ à optimiser. Supposons que l’on choisisse comme test d’arrêt dans l’algorithme de descente modèle, le critère idéal : "$ x_k = x^\ast $". Dans un monde idéal (i.e. en supposant tous les calculs exacts et la capacité de calcul illimitée), soit l’algorithme s’arrête après un nombre fini d’itérations, soit il construit (théoriquement) une suite infinie $ x_0,x_1,\ldots,x_k,\ldots $ de points de $ \mathbb{R}^n $ qui converge vers $ x^\ast $.
f899c72b 431\newline
6079168f 432En pratique, un test d’arrêt devra être choisi pour garantir que l’algorithme s’arrête toujours après un nombre fini d’itérations et que le dernier point calculé soit suffisamment proche de $ x^\ast $.
8000c039 433
6079168f 434Soit $ \varepsilon \in \mathbb{R}_{+}^{*} $ la précision demandée. Plusieurs critères sont à notre disposition : tout d’abord (et c’est le plus naturel), un critère d’optimalité basé sur les conditions nécessaires d’optimalité du premier ordre : en optimisation différentiable sans contrainte, on testera si
8000c039
JB
435$$ \norme{\nabla J(x_k)} < \varepsilon, $$
436auquel cas l’algorithme s’arrête et fournit l’itéré courant $ x_k $ comme solution.
437
e17152e0 438En pratique, le test d’optimalité n’est pas toujours satisfait et on devra faire appel à d’autres critères fondés sur l’expérience du numérique :
8000c039
JB
439\begin{itemize}
440 \item Stagnation de la solution : $ \norme{x_{k+1} - x_k} < \varepsilon(1 + \norme{x_k}) $.
441 \item Stagnation de la valeur courante : $ |J(x_{k+1}) - J(x_k)| < \varepsilon(1 + |J (x_k)|) $.
442 \item Nombre d’itérations dépassant un seuil fixé à l’avance : $ k < IterMax $.
443\end{itemize}
444et généralement une combinaison de ces critères :
de30386e 445\newline
8000c039
JB
446\newline
447Critère d’arrêt =
448\begin{tabular}{l}
9f054a9c 449 Test d’optimalité satisfait \\
8000c039
JB
450 OU (Stagnation de la valeur courante ET Stagnation de la solution) \\
451 OU Nombre d’itérations maximum autorisé dépassé
452\end{tabular}
6ec0df37 453
8000c039 454\subsection{La recherche linéaire}
e37aec65 455
8000c039
JB
456Supposons pour l’instant résolu le problème du choix de la direction de descente et intéressons nous uniquement au calcul du pas : c’est la phase de recherche linéaire.
457\newline
458Soit $ x_0 \in \mathbb{R}^n $ un point non critique et $ d $ une direction de descente de $ J $ en $ x_0 $. Nous cherchons à calculer un pas $ s > 0 $ de sorte que :
459$$ J(x_0 + sd) < J(x_0). $$
6079168f 460Le choix de ce pas répond généralement à deux objectifs souvent contradictoires : trouver le meilleur pas possible et effectuer le moins de calculs possibles. Ces deux objectifs ont donné naissance à deux grandes familles : les algorithmes à pas fixe et ceux à pas optimal.
8000c039
JB
461
462\hrulefill
463\newline
464RECHERCHE LINÉAIRE : PAS FIXE. $ s_k = s_{k-1} $
465
466\hrulefill
467
468\hrulefill
469\newline
470RECHERCHE LINÉAIRE : PAS OPTIMAL. $ s_k $ solution du problème $ \displaystyle\min_{s \in \mathbb{R}_{+}^{*}} J(x_k + sd_k) $
e37aec65 471
8000c039
JB
472\hrulefill
473\newline
a5f09ac4 474Illustrées par les méthodes de descente de gradient, aucune de ces deux stratégies ne s’est révélée réellement convaincante : si la première peut être “risquée” du point de vue de la convergence, la seconde est souvent loin d’être triviale à mettre en oeuvre (sauf dans le cas quadratique) et généralement inutilement coûteuse : en effet, à quoi bon calculer très précisément un pas optimal dans une direction qui n’est peut-être pas la bonne ? (comme c’est par exemple le cas pour la méthode de plus profonde descente). Les recherches linéaires modernes reposent sur l’idée qu’un pas de descente acceptable est un pas qui fait “suffisamment” décroître la fonction objectif. Reste alors à définir les pas qui sont acceptables et ceux qui ne le sont pas.
8000c039
JB
475\begin{Def}
476 On appelle $ \varphi : s \in \mathbb{R} \longmapsto J(x + sd)$ la fonction mérite associée à $ J $ en $ x $.
477\end{Def}
478\begin{Def}
9f054a9c 479 Dans le cas où $ J $ est différentiable sur $ \mathcal{C} $, on dit que un algorithme de descente converge ssi
7590f4bb 480 $$ \forall x_0 \in \mathbb{R}^n \lim\limits_{k \rightarrow +\infty} \norme{\nabla J(x_k)} = 0 $$
8000c039
JB
481\end{Def}
482
483\subsubsection{Principe de démonstration de convergence}
de30386e 484
e17152e0 485Une technique classique en optimisation pour obtenir des résultats de convergence globale consiste à montrer que l’algorithme de descente considéré vérifie une inégalité du type :
8000c039 486$$ J(x_k) - J(x_{k+1}) \geq c\norme{\nabla J(x_k)}^2, $$
e17152e0 487où $ c $ est une constante réelle.
8000c039
JB
488\newline
489En sommant ces inégalités pour $ k $ variant de $ 0 $ à $ N - 1 $, on obtient :
490$$ \forall N \in \mathbb{N} \ J(x_0) - J(x_N) \geq c \sum_{i=0}^{N-1}\norme{\nabla J(x_i)}^2 $$
9f054a9c 491Si $ 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 :
8000c039 492$$ \lim\limits_{k \rightarrow +\infty} \norme{\nabla J(x_k)} = 0 $$
7590f4bb
JB
493\begin{Def}
494 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 :
495 \begin{itemize}
496 \item linéaire si l'erreur décroît linéairement i.e. :
497 $$ \exists \tau \in ]0,1[ \ \lim_{k \rightarrow +\infty} \frac{\norme{x_{k+1} - x^\ast}}{\norme{x_k - x^\ast}} = \tau $$
498 \item superlinéaire si :
499 $$ \lim_{k \rightarrow +\infty} \frac{\norme{x_{k+1} - x^\ast}}{\norme{x_k - x^\ast}} = 0 $$
500 \item d'ordre $ p $ si :
501 $$ \exists \tau \geq 0 \ \lim_{k \rightarrow +\infty} \frac{\norme{x_{k+1} - x^\ast}}{\norme{x_k - x^\ast}^p} = \tau $$
502 En particulier, si $ p = 2 $, la convergence est dite quadratique.
503 \end{itemize}
504\end{Def}
505L'é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.
8000c039
JB
506
507\section{Méthode Newtonienne}
508
e17152e0 509Les hypothèses sur $ \mathcal{P} $ de la section précédente restent les mêmes dans cette section. L’algorithme de Newton en optimisation est une application directe de l’algorithme de Newton pour la résolution d’équations du type : $ F(x) = 0 $. En optimisation sans contrainte, l’algorithme de Newton cherche les solutions de l’équation :
8000c039
JB
510$$ \nabla J(x) = 0, $$
511autrement dit, les points critiques de la fonction $ J $ à minimiser.
512\newline
513En supposant $ J $ de classe $ \mathcal{C}^2 $ et la matrice hessienne $ H[J](x_k) $ inversible, une itération de l’algorithme de Newton s’écrit :
514$$ x_{k+1} = x_k - H[J](x_k)^{-1} \nabla J(x_k), $$
682e0379
JB
515où $ 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 :
516$$ \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 $$
e17152e0 517Autrement dit, $ d_k $ est le point de minimum global de l’approximation de second ordre de $ J $ au voisinage du point courant $ x_k $.
a5f09ac4 518À 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 $.
7fd51f8b
JB
519\newline
520Les propriétés remarquables de cet algorithme sont :
682e0379
JB
521
522\begin{tabular}{|p{20em}|p{20em}|}
523 \hline
524 Avantages & Inconvénients \\
525 \hline
526 sa convergence quadratique (le nombre de décimales exactes est multiplié par 2 à chaque itération). & \\
527 \hline
528 & les difficultés et le coût de calcul de la hessienne $ H[J](x_k) $ : l’expression analytique des dérivées secondes est rarement disponible dans les applications. \\
529 \hline
530 & le coût de résolution du système linéaire $ H[J](x_k )(x_{k+1} - x_k) = \nabla J(x_k) $. \\
531 \hline
532 & l’absence de convergence si le premier itéré est trop loin de la solution, ou si la hessienne est singulière. \\
533 \hline
534 & pas de distinction entre minima, maxima et points stationnaires. \\
535 \hline
536\end{tabular}
537\newline
6079168f 538La question que l’on se pose est donc : comment forcer la convergence globale de l’algorithme de Newton ? L’idée des méthodes de type Newton consiste à reprendre l’algorithme de Newton en remplaçant les itérations par :
682e0379
JB
539$$ x_{k+1} = x_k - s_k H_k^{-1} \nabla J(x_k), $$
540
541\begin{itemize}
542 \item la matrice $ H_k $ est une approximation de la hessienne $ H[J](x_k) $.
543 \item $ s_k > 0 $ est le pas calculé par une recherche linéaire bien choisie.
544\end{itemize}
545Plusieurs questions se posent alors :
546\begin{itemize}
547 \item Comment déterminer une matrice $ H_k $ qui soit une “bonne” approximation de la hessienne à l’itération $ k $ sans utiliser les informations de second ordre et garantir que $ H_k^{-1} \nabla J(x_k) $ soit bien une direction de descente de $ J $ en $ x_k $, sachant que la direction de Newton, si elle existe, n’en est pas nécessairement une ?
548 \item Comment conserver les bonnes propriétés de l’algorithme de Newton ?
549\end{itemize}
7fd51f8b 550Nous ne répondrons pas à ces questions qui sont hors du cadre de ce projet. Cette section permet d'introduire certains prérequis pour l'étude de la méthode PQS et de rendre compte de sa filiation.
de30386e 551
8000c039 552\section{Méthode PQS (ou SQP)}
e37aec65 553
7590f4bb 554Nous 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é}.
682e0379 555
329ebbc2
JB
556\subsection{Problème quadratique sous contraintes linéaires}
557
558Nous 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.
559\newline
560Ce type de problème quadratique se pose sous la forme :
561$$
562 \mathcal{PQ} \left \{
563 \begin{array}{l}
564 \displaystyle\min_{x \in \mathbb{R}^n} c^\top x + \frac{1}{2} x^\top \mathcal{Q} x \\
565 A^\top x + b \leq 0 \\
566 A^{\prime^\top} x + b^\prime = 0
567 \end{array}
568 \right .
569$$
570où $$ \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 $$
571Or
572$$ 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 $$
573Donc le problème se ramène à :
574
575\subsubsection{Algorithme 1}
576
577\subsubsection{Algorithme 2}
578
a5f09ac4 579\subsection{Algorithmes Newtoniens}
682e0379 580
6079168f 581Les 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} $.
682e0379
JB
582
583\subsection{Algorithme PQS}
584
585\subsubsection{Contraintes d’égalité}
586
587Considérons un problème d’optimisation différentiable $ \mathcal{P} $ avec contraintes d’égalité :
588$$
589 \mathcal{P} \left \{
7fd51f8b 590 \begin{array}{l}
682e0379
JB
591 \displaystyle\min_{x \in \mathbb{R}^n} J(x) \\
592 h(x) = 0
593 \end{array}
594 \right .
595$$
596où $ J: \mathbb{R}^n \longrightarrow \mathbb{R} $ et $h: \mathbb{R}^n \longrightarrow \mathbb{R}^q$ sont supposées au moins différentiables.
597\newline
598Les conditions d’optimalité de Lagrange (ou \textit{KKT}) s’écrivent :
7fd51f8b 599$$ \nabla J(x) + \sum\limits_{i=1}^{q} \lambda_i \nabla h_i(x) = 0 \iff \nabla L(x,\lambda) = 0 $$
682e0379
JB
600donc $ \mathcal{P} $ devient :
601$$ \begin{pmatrix}
7fd51f8b 602 \nabla J(x) + \sum\limits_{i=1}^{q} \lambda_i \nabla h_i(x) \\
682e0379
JB
603 h(x)
604 \end {pmatrix} = 0 $$
605Pour résoudre ce système d’équations, utilisons la méthode de Newton dont une itération s’écrit ici :
7fd51f8b
JB
606$$ H[L](x_k,\lambda_k)\begin{pmatrix}
607 x_{k+1} - x_k \\
608 \lambda_{k+1} - \lambda_k
609 \end{pmatrix} = -\nabla L(x_k,\lambda_k) $$
610soit :
611$$ \begin{pmatrix}
612 H_x[L](x_k,\lambda_k) & D_h(x_k)^\top \\
613 D_h(x_k) & 0
614 \end{pmatrix} \begin{pmatrix}
615 x_{k+1} - x_k \\
616 \lambda_{k+1} - \lambda_k
617 \end{pmatrix} = -\begin{pmatrix}
618 \nabla_x L(x_k,\lambda_k) \\
619 h(x_k)
620 \end{pmatrix} $$
621où $ D_h(x) $ désigne la matrice jacobienne de l’application $ h : \mathbb{R}^n \longrightarrow \mathbb{R}^q $ définie par :
7590f4bb 622$$ D_h(x)^\top = \begin{bmatrix} \nabla h_1(x)^\top\ldots\nabla h_q(x)^\top \end{bmatrix} $$
7fd51f8b
JB
623Posons : $ 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 :
624$$ \begin{pmatrix}
625 H_k & D_h(x_k)^\top \\
626 D_h(x_k) & 0
627 \end{pmatrix} \begin{pmatrix}
628 d \\
629 \mu - \lambda_k
630 \end{pmatrix} = -\begin{pmatrix}
631 \nabla_x L(x_k,\lambda_k) \\
632 h(x_k)
633 \end{pmatrix} $$
634et est bien définie à condition que la matrice $ H_x[L](x_k,\lambda_k) $ soit inversible. Ce sera le cas si :
635\begin{enumerate}[label=(\roman*)]
7590f4bb 636 \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.
7fd51f8b
JB
637 \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é.
638\end{enumerate}
639Revenons à l’itération. Elle s’écrit encore :
640$$
641 \left \{
642 \begin{array}{r c l}
643 H_kd + \sum\limits_{i=1}^q(\mu_i - \lambda_{k_i})\nabla h_i(x_k) & = & -\nabla_x L(x_k,\lambda_k) \\
644 \nabla h_i(x_k)^\top d + h_i(x_k) & = & 0, \ \forall i \in \{1,\ldots,q\}
645 \end{array}
646 \right .
647$$
e17152e0 648Or $ \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ù :
7fd51f8b
JB
649$$
650 \left \{
651 \begin{array}{r c l}
652 H_kd + \sum\limits_{i=1}^q\mu_i\nabla h_i(x_k) & = & -\nabla J(x_k) \\
653 \nabla h_i(x_k)^\top d + h_i(x_k) & = & 0, \ \forall i \in \{1,\ldots,q\}
654 \end{array}
655 \right .
656$$
6079168f 657On reconnait dans le système ci-dessus les conditions d’optimalité de Lagrange du problème quadratique suivant :
7fd51f8b
JB
658$$
659 \mathcal{PQ}_k \left \{
660 \begin{array}{l}
661 \displaystyle\min_{d \in \mathbb{R}^n} \nabla J(x_k)^\top d + \frac{1}{2}d^\top H_k d \\
662 h_i(x_k) + \nabla h_i(x_k)^\top d = 0, \ \forall i \in \{1,\ldots,q\}
663 \end{array}
664 \right .
665$$
666Le 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.
667\newline
668Comme 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 :
669
e17152e0
JB
670\hrulefill
671\newline
6079168f 672ALGORITHME PQS AVEC CONSTRAINTES D'ÉGALITÉ.
e17152e0
JB
673\newline
674\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.
675\newline
676\textit{Sortie}: une approximation $ x_k $ de la solution $ x^\ast $ du problème $ \mathcal{P} $.
677\begin{enumerate}
678 \item $ k := 0 $.
679 \item Tant que $ \norme{\nabla L(x_k,\lambda_k)} > \varepsilon $,
680 \begin{enumerate}
681 \item Résoudre le sous-problème quadratique :
682 $$
683 \mathcal{PQ}_k \left \{
684 \begin{array}{l}
685 \displaystyle\min_{d \in \mathbb{R}^n} \nabla J(x_k)^\top d + \frac{1}{2}d^\top H_k d \\
686 h_i(x_k) + \nabla h_i(x_k)^\top d = 0, \ \forall i \in \{1,\ldots,q\}
687 \end{array}
688 \right .
689 $$
690 et obtenir la solution primale $ d_k $ et le multiplicateur $ \lambda^{\prime} $ associé à la contrainte d’égalité.
691 \item $ x_{k+1} = x_k + d_k; \ \lambda_{k+1} = \lambda^{\prime}; \ k := k + 1 $.
692 \end{enumerate}
693 \item Retourner $ x_k $.
694\end{enumerate}
695
696\hrulefill
697
7fd51f8b 698\subsubsection{Contraintes d’inégalité}
682e0379 699
e17152e0
JB
700Intéressons nous maintenant aux problèmes avec contraintes d’égalité et d’inégalité :
701$$
702 \mathcal{P} \left \{
703 \begin{array}{l}
704 \displaystyle\min_{x \in \mathbb{R}^n} J(x) \\
705 g(x) \leq 0 \\
706 h(x) = 0
707 \end{array}
708 \right .
709$$
710où $ 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.
711\newline
7590f4bb 712Selon 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 :
e17152e0 713$$ 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 $$
7590f4bb
JB
714Soit à l'ordre 2 pour le Lagrangien :
715$$ 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) $$
716et à l'ordre 1 pour les contraintes :
717$$ g(x) \approx g(x_k) + \nabla g(x_k)^\top(x - x_k) $$
718$$ h(x) \approx h(x_k) + \nabla h(x_k)^\top(x - x_k) $$
719En 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 $ :
e17152e0
JB
720
721\hrulefill
722\newline
6079168f 723ALGORITHME PQS AVEC CONSTRAINTES D'ÉGALITÉ ET D'INEGALITÉ.
e17152e0
JB
724\newline
725\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.
726\newline
727\textit{Sortie}: une approximation $ x_k $ de la solution $ x^\ast $ du problème $ \mathcal{P} $.
728\begin{enumerate}
729 \item $ k := 0 $.
730 \item Tant que $ \norme{\nabla L(x_k,\lambda_k,\mu_k)} > \varepsilon $,
731 \begin{enumerate}
732 \item Résoudre le sous-problème quadratique :
733 $$
734 \mathcal{PQ}_k \left \{
735 \begin{array}{l}
736 \displaystyle\min_{d \in \mathbb{R}^n} \nabla J(x_k)^\top d + \frac{1}{2}d^\top H_k d \\
19959015 737 g_j(x_k) + \nabla g_j(x_k)^\top d \leq 0, \ \forall j \in \{1,\ldots,p\} \\
e17152e0
JB
738 h_i(x_k) + \nabla h_i(x_k)^\top d = 0, \ \forall i \in \{1,\ldots,q\}
739 \end{array}
740 \right .
741 $$
742 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.
743 \item $ x_{k+1} = x_k + d_k; \ \lambda_{k+1} = \lambda^{\prime}; \ \mu_{k+1} = \mu^{\prime}; \ k := k + 1 $.
744 \end{enumerate}
745 \item Retourner $ x_k $.
746\end{enumerate}
747
748\hrulefill
749\newline
750Afin 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.
751\newline
752Etant 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.
753
7590f4bb
JB
754\subsection{Stratégie d'approximation de la hessienne}
755
756\subsubsection{Équation de sécante et approximation}
757
758L'approximation $ H_k $ de la hessienne du Lagrangien peut être obtenu par la relation :
759$$ \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) $$
760On construit une approximation $ H_{k+1} $ de $ H[L](x_{k+1},\lambda_{k+1},\mu_{k+1}) $ comme solution de l’équation :
761$$ 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}) $$
762appelée équation de sécante ou équation de quasi-Newton.
763\newline
764De 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 :
765$$ 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 $$
766Dans 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.
767\newline
768Une 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 :
769$$ H_{k+1} = H_k + U_k $$
770
771\subsubsection{Mises à jour DFP et BFGS}
772
773\subsection{Exemple d'utilisation de PQS}
774
775Considérons le problème $ \mathcal{P} $ suivant :
776$$
777 \mathcal{P} \left \{
778 \begin{array}{l}
779 \displaystyle\min_{(x,y,z) \in \mathbb{R}^3} J(x,y,z) = x^2 + y^2 + z^2 -r^2 \\
780 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 \\
781 \end{array}
782 \right .
783$$
19959015 784où $$ (r,r_1,r_2) \in \mathbb{R}_+^{*^3} \land r < r_1 \land r < r_2. $$
84af1fe2 785\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.
7590f4bb 786\newline
84af1fe2 787Le 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). $$
7590f4bb 788\newline
84af1fe2 789Le 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). $$
7590f4bb 790\newline
84af1fe2
JB
791Le gradient de $ g $ : $$ \nabla g(x,y,z) = (\nabla g_1(x,y,z),\nabla g_2(x,z,z)) $$
792$$ = ((\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)) $$
793$$ = ((2x,2y,0),(2x,0,2z)). $$
7590f4bb 794\newline
84af1fe2
JB
795Le gradient du Lagrangien $ L $ :
796$$ \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)) $$
797\newline
798La matrice hessienne de $ J $ : $$ H[J](x,y,z) =
7590f4bb 799 \begin{pmatrix}
46973afb
JB
800 \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) \\
801 \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) \\
802 \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) \\
7590f4bb
JB
803 \end{pmatrix} =
804 \begin{pmatrix}
84af1fe2
JB
805 2 & 0 & 0 \\
806 0 & 2 & 0 \\
807 0 & 0 & 2 \\
808 \end{pmatrix} = 2Id_{\mathbb{R}^3} $$
329ebbc2 809On 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} $.
7590f4bb 810
e37aec65
JB
811\bibliographystyle{plain}
812\bibliography{stdlib_sbphilo}
813
814%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
815
816\end{document}