.. _td1acenoncesession9rst: ====================================== 1A.algo - Optimisation sous contrainte ====================================== .. only:: html **Links:** :download:`notebook `, :downloadlink:`html `, :download:`python `, :downloadlink:`slides `, :githublink:`GitHub|_doc/notebooks/td1a_algo/td1a_cenonce_session9.ipynb|*` L’optimisation sous contrainte est un problème résolu. Ce notebook utilise une librairie externe et la compare avec l’algorithme *Arrow-Hurwicz* qu’il faudra implémenter. Plus de précision dans cet article `Damped Arrow-Hurwicz algorithm for sphere packing `__. .. code:: ipython3 from jyquickhelper import add_notebook_menu add_notebook_menu() .. contents:: :local: Le langage Python propose des modules qui permettent de résoudre des problèmes d’optimisation sous contraintes et il n’est pas forcément nécessaire de connaître la théorie derrière les algorithmes de résolution pour s’en servir. Au cours de cette séance, on verra comment faire. Même si comprendre comment utiliser une fonction d’un module tel que `cvxopt `__ requiert parfois un peu de temps et de lecture. On verra également un algorithme simple d’optimisation. C’est une bonne façon de comprendre que cela prend du temps si on veut implémenter soi-même ce type de solution tout en étant aussi rapide et efficace. Exercice 1 : optimisation avec cvxopt ------------------------------------- On souhaite résoudre le problème d’optimisation suivant : :math:`\left \{ \begin{array}{l} \min_{x,y} \left \{ x^2 + y^2 - xy + y \right \} \\ sous \; contrainte \; x + 2y = 1 \end{array}\right .` Le module `cvxopt `__ est un des modules les plus indiqués pour résoudre ce problème. Voici quelques instructions qui l’utilisent : .. code:: ipython3 from cvxopt import solvers, matrix m = matrix( [ [2.0, 1.1] ] ) # mettre des réels (float) et non des entiers # cvxopt ne fait pas de conversion implicite t = m.T # la transposée t.size # affiche les dimensions de la matrice .. parsed-literal:: (1, 2) La documentation ``cvxopt`` est parfois peu explicite. Il ne faut pas hésiter à regarder les exemples d’abord et à la lire avec attention les lignes qui décrivent les valeurs que doivent prendre chaque paramètre de la fonction. Le plus intéressant pour le cas qui nous intéresse est celui-ci (tiré de la page `problems with nonlinear objectives `__) : .. code:: ipython3 from cvxopt import solvers, matrix, spdiag, log def acent(A, b): m, n = A.size def F(x=None, z=None): if x is None: # l'algorithme fonctionne de manière itérative # il faut choisir un x initial, c'est ce qu'on fait ici return 0, matrix(1.0, (n,1)) if min(x) <= 0.0: return None # cas impossible # ici commence le code qui définit ce qu'est une itération f = -sum(log(x)) Df = -(x**-1).T if z is None: return f, Df H = spdiag(z[0] * x**(-2)) return f, Df, H return solvers.cp(F, A=A, b=b)['x'] A = matrix ( [[1.0,2.0]] ).T b = matrix ( [[ 1.0 ]] ) print(acent(A,b)) .. parsed-literal:: pcost dcost gap pres dres 0: 0.0000e+00 0.0000e+00 1e+00 1e+00 1e+00 1: 9.9000e-01 4.6251e+00 1e-02 2e+00 7e+01 2: 3.6389e+00 3.9677e+00 1e-04 1e-01 3e+01 3: 3.0555e+00 3.3406e+00 1e-06 1e-01 2e+01 4: 2.5112e+00 2.7758e+00 1e-08 1e-01 8e+00 5: 2.1118e+00 2.3358e+00 1e-10 1e-01 4e+00 6: 1.9684e+00 2.1118e+00 1e-12 6e-02 1e+00 7: 2.0493e+00 2.0796e+00 1e-14 1e-02 1e-01 8: 2.0790e+00 2.0794e+00 1e-16 2e-04 2e-03 9: 2.0794e+00 2.0794e+00 1e-18 2e-06 2e-05 10: 2.0794e+00 2.0794e+00 1e-20 2e-08 2e-07 11: 2.0794e+00 2.0794e+00 1e-22 2e-10 2e-09 Optimal solution found. [ 5.00e-01] [ 2.50e-01] .. code:: ipython3 # il existe un moyen d'éviter l'affichage des logs (pratique si on doit faire # un grand nombre d'optimisation) from cvxopt import solvers solvers.options['show_progress'] = False print(acent(A,b)) solvers.options['show_progress'] = True .. parsed-literal:: [ 5.00e-01] [ 2.50e-01] Cet exemple résoud le problème de minimisation suivant : :math:`\left\{ \begin{array}{l} \min_{X} \left\{ - \sum_{i=1}^n \ln x_i \right \} \\ sous \; contrainte \; AX = b \end{array} \right.` Les deux modules `numpy `__ et `cvxopt `__ n’utilisent pas les mêmes matrices (les mêmes objets ``matrix``) bien qu’elles portent le même nom dans les deux modules. Les fonctions de ``cvxopt`` ne fonctionnent qu’avec les matrices de ce module. Il ne faut pas oublier de convertir la matrice quand elle est décrite par une autre classe. .. code:: ipython3 import cvxopt m = cvxopt.matrix( [[ 0, 1.5], [ 4.5, -6] ] ) print(m) .. parsed-literal:: [ 0.00e+00 4.50e+00] [ 1.50e+00 -6.00e+00] Ces informations devraient vous permettre de résoudre le premier problème. Exercice 2 : l’algorithme de Arrow-Hurwicz ------------------------------------------ On cherche à résoudre le même problème avec l’algorithme de `Arrow-Hurwicz `__ (tiré de `Introduction à l’optimisation `__ de Jean-Christophe Culioli, voir également l’\ `algorithme d’Uzawa `__) et les notations suivantes : :math:`\left\{ \begin{array}{l} \min_U J(U) = u_1^2 + u_2^2 - u_1 u_2 + u_2 \\ sous \; contrainte \; \theta(U) = u_1 + 2u_2 - 1 = 0 \end{array} \right.` L’algorithme est le suivant : - On choisit :math:`\epsilon > 0`, :math:`\rho > 0`, des vecteurs aléatoires :math:`U_0 \in \mathbb{R}^2` et :math:`P_0 \in \mathbb{R}^d` (:math:`d` est le nombre de contraintes). - A l’itération :math:`k`, on met à jour : :math:`\begin{array}{l} U_{t+1} = U_t - \epsilon \left( \nabla J (U_t) + \left[ \nabla \theta(U_t) \right] ' P_t \right) \\ P_{t+1} = P_t + \rho \theta( U_{t+1} ) \end{array}` - On retourne à l’étape précédente jusqu’à ce que la suite :math:`(U_k)` n’évolue plus. Le coefficient :math:`P_t` correspond au coefficient de Lagrange pour un Lagrangien défini comme suit : :math:`L(U,P) = J(U) + P' \theta(U)` La suite :math:`(U_t)_t` converge vers la solution en se déplaçant le long du gradient de la fonction :math:`J` lorsque la contrainte est vérifiée. Lorsqu’elle ne l’est pas, cette direction est un mélange du gradient de la fonction à optimiser et de celui de la contrainte à respecter. Implémenter cet algorithme et vérifier qu’il converge vers la même solution que celle obtenue pour le premier problème. Exercice 3 : Le lagrangien augmenté ----------------------------------- On reprend l’algorithme précédent mais appliqué à la fonction : :math:`A(U) = J(U) + \frac{c}{2} \theta^2(U)` L’objectif est ici de remplacer la fonction :math:`J` par la fonction :math:`A` (ou Lagrangien augmenté) dans l’algorithme précédent et de vérifier qu’il converge en un nombre moindre d’itérations lorsqu’on choisit les mêmes conditions initiales. Prolongement 1 : inégalité -------------------------- On souhaite maintenant ajouter une contrainte d’inégalité : :math:`Du -d \leqslant 0`. L’idée est de créer une matrice :math:`\Theta` définie par bloc : :math:`\Theta = \left( \begin{array}{c} \theta \\ D \end{array} \right) \; , \; \Pi = \left( \begin{array}{c} P \\ D \end{array} \right) \; , \; F = diag \left( \begin{array}{c} 1 \\ r_i \end{array} \right)` On considère que :math:`\Theta` est une matrice. On étend ce découpage en bloc au formule de mise à jour de l’algorithme de Arrow-Hurwicz : :math:`\begin{array}{l} U_{t+1} = U_t - \epsilon \left( \nabla J (U_t) + \left[ \nabla \Theta(U_t) \right] ' Fi_t \right) \\ \Pi_{t+1} = F \left( \Pi_t + \rho \theta ( U_{t+1} ) \right) \end{array}` Les coefficients :math:`r_i` valent 0 ou 1 selon que la contrainte d’inégalité :math:`i` est respectée ou non. On applique la modification pour le problème : :math:`\left\{ \begin{array}{l} \min_U J(U) = u_1^2 + u_2^2 - u_1 u_2 + u_2 \\ sous \; contrainte \; \theta(U) = u_1 + 2u_2 - 1 = 0 \; et \; u_1 \geqslant 0.3 \end{array}\right.` Prolongement 2 : optimisation d’une fonction linéaire ----------------------------------------------------- On s’intéresse à un problème différent : :math:`\left\{ \begin{array}{l} \min_X C'X \\ sous \; contrainte \; AX \leqslant B \end{array}\right.` Une façon de résoudre ce problème est l’\ `algorithme de Karmarkar `__. On optimise maintenant une fonction linéaire et non une fonction quadratique. Par conséquent, la solution est forcément sur un bord : des contraintes sont saturées. Vous pouvez découvrir un autre exemple d’optimisation dynamique au travers de ce blog : `Optimisation sous contraintes appliquée au calcul du report des voix `__.