{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["# 1A.algo - Optimisation sous contrainte\n", "\n", "L'optimisation sous contrainte est un probl\u00e8me r\u00e9solu. Ce notebook utilise une librairie externe et la compare avec l'algorithme *Arrow-Hurwicz* qu'il faudra impl\u00e9menter. Plus de pr\u00e9cision dans cet article [Damped Arrow-Hurwicz algorithm for sphere packing](https://arxiv.org/abs/1605.05473)."]}, {"cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [{"data": {"text/html": ["<div id=\"my_id_menu_nb\">run previous cell, wait for 2 seconds</div>\n", "<script>\n", "function repeat_indent_string(n){\n", "    var a = \"\" ;\n", "    for ( ; n > 0 ; --n)\n", "        a += \"    \";\n", "    return a;\n", "}\n", "var update_menu_string = function(begin, lfirst, llast, sformat, send, keep_item, begin_format, end_format) {\n", "    var anchors = document.getElementsByClassName(\"section\");\n", "    if (anchors.length == 0) {\n", "        anchors = document.getElementsByClassName(\"text_cell_render rendered_html\");\n", "    }\n", "    var i,t;\n", "    var text_menu = begin;\n", "    var text_memo = \"<pre>\\nlength:\" + anchors.length + \"\\n\";\n", "    var ind = \"\";\n", "    var memo_level = 1;\n", "    var href;\n", "    var tags = [];\n", "    var main_item = 0;\n", "    var format_open = 0;\n", "    for (i = 0; i <= llast; i++)\n", "        tags.push(\"h\" + i);\n", "\n", "    for (i = 0; i < anchors.length; i++) {\n", "        text_memo += \"**\" + anchors[i].id + \"--\\n\";\n", "\n", "        var child = null;\n", "        for(t = 0; t < tags.length; t++) {\n", "            var r = anchors[i].getElementsByTagName(tags[t]);\n", "            if (r.length > 0) {\n", "child = r[0];\n", "break;\n", "            }\n", "        }\n", "        if (child == null) {\n", "            text_memo += \"null\\n\";\n", "            continue;\n", "        }\n", "        if (anchors[i].hasAttribute(\"id\")) {\n", "            // when converted in RST\n", "            href = anchors[i].id;\n", "            text_memo += \"#1-\" + href;\n", "            // passer \u00e0 child suivant (le chercher)\n", "        }\n", "        else if (child.hasAttribute(\"id\")) {\n", "            // in a notebook\n", "            href = child.id;\n", "            text_memo += \"#2-\" + href;\n", "        }\n", "        else {\n", "            text_memo += \"#3-\" + \"*\" + \"\\n\";\n", "            continue;\n", "        }\n", "        var title = child.textContent;\n", "        var level = parseInt(child.tagName.substring(1,2));\n", "\n", "        text_memo += \"--\" + level + \"?\" + lfirst + \"--\" + title + \"\\n\";\n", "\n", "        if ((level < lfirst) || (level > llast)) {\n", "            continue ;\n", "        }\n", "        if (title.endsWith('\u00b6')) {\n", "            title = title.substring(0,title.length-1).replace(\"<\", \"&lt;\")\n", "         .replace(\">\", \"&gt;\").replace(\"&\", \"&amp;\");\n", "        }\n", "        if (title.length == 0) {\n", "            continue;\n", "        }\n", "\n", "        while (level < memo_level) {\n", "            text_menu += end_format + \"</ul>\\n\";\n", "            format_open -= 1;\n", "            memo_level -= 1;\n", "        }\n", "        if (level == lfirst) {\n", "            main_item += 1;\n", "        }\n", "        if (keep_item != -1 && main_item != keep_item + 1) {\n", "            // alert(main_item + \" - \" + level + \" - \" + keep_item);\n", "            continue;\n", "        }\n", "        while (level > memo_level) {\n", "            text_menu += \"<ul>\\n\";\n", "            memo_level += 1;\n", "        }\n", "        text_menu += repeat_indent_string(level-2);\n", "        text_menu += begin_format + sformat.replace(\"__HREF__\", href).replace(\"__TITLE__\", title);\n", "        format_open += 1;\n", "    }\n", "    while (1 < memo_level) {\n", "        text_menu += end_format + \"</ul>\\n\";\n", "        memo_level -= 1;\n", "        format_open -= 1;\n", "    }\n", "    text_menu += send;\n", "    //text_menu += \"\\n\" + text_memo;\n", "\n", "    while (format_open > 0) {\n", "        text_menu += end_format;\n", "        format_open -= 1;\n", "    }\n", "    return text_menu;\n", "};\n", "var update_menu = function() {\n", "    var sbegin = \"\";\n", "    var sformat = '<a href=\"#__HREF__\">__TITLE__</a>';\n", "    var send = \"\";\n", "    var begin_format = '<li>';\n", "    var end_format = '</li>';\n", "    var keep_item = -1;\n", "    var text_menu = update_menu_string(sbegin, 2, 4, sformat, send, keep_item,\n", "       begin_format, end_format);\n", "    var menu = document.getElementById(\"my_id_menu_nb\");\n", "    menu.innerHTML=text_menu;\n", "};\n", "window.setTimeout(update_menu,2000);\n", "            </script>"], "text/plain": ["<IPython.core.display.HTML object>"]}, "execution_count": 2, "metadata": {}, "output_type": "execute_result"}], "source": ["from jyquickhelper import add_notebook_menu\n", "add_notebook_menu()"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Le langage Python propose des modules qui permettent de r\u00e9soudre des probl\u00e8mes d'optimisation sous contraintes et il n'est pas forc\u00e9ment n\u00e9cessaire de conna\u00eetre la th\u00e9orie derri\u00e8re les algorithmes de r\u00e9solution pour s'en servir. Au cours de cette s\u00e9ance, on verra comment faire. M\u00eame si comprendre comment utiliser une fonction d'un module tel que  [cvxopt](http://cvxopt.org/) requiert parfois un peu de temps et de lecture. \n", "\n", "On verra \u00e9galement un algorithme simple d'optimisation. C'est une bonne fa\u00e7on de comprendre que cela prend du temps si on veut impl\u00e9menter soi-m\u00eame ce type de solution tout en \u00e9tant aussi rapide et efficace."]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Exercice 1 : optimisation avec cvxopt\n", "\n", "On souhaite r\u00e9soudre le probl\u00e8me d'optimisation suivant :\n", "\n", "$\\left \\{ \\begin{array}{l} \\min_{x,y} \\left \\{ x^2 + y^2 - xy + y \\right \\}  \\\\ sous \\; contrainte \\; x + 2y = 1 \\end{array}\\right .$\n", "\n", "Le module [cvxopt](http://cvxopt.org/userguide/solvers.html\\#problems-with-nonlinear-objectives) est un des modules les plus indiqu\u00e9s pour r\u00e9soudre ce probl\u00e8me. Voici quelques instructions qui l'utilisent :"]}, {"cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [{"data": {"text/plain": ["(1, 2)"]}, "execution_count": 3, "metadata": {}, "output_type": "execute_result"}], "source": ["from cvxopt import solvers, matrix\n", "m = matrix( [ [2.0, 1.1] ] )  # mettre des r\u00e9els (float) et non des entiers\n", "                              # cvxopt ne fait pas de conversion implicite\n", "t = m.T                       # la transpos\u00e9e\n", "t.size                        # affiche les dimensions de la matrice"]}, {"cell_type": "markdown", "metadata": {}, "source": ["La documentation ``cvxopt`` est parfois peu explicite. Il ne faut pas h\u00e9siter \u00e0 regarder les exemples d'abord et \u00e0 la lire avec attention les lignes qui d\u00e9crivent les valeurs que doivent prendre chaque param\u00e8tre de la fonction. Le plus int\u00e9ressant pour le cas qui nous int\u00e9resse est celui-ci (tir\u00e9 de la page [problems with nonlinear objectives](http://cvxopt.org/userguide/solvers.html#problems-with-nonlinear-objectives)) :"]}, {"cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["     pcost       dcost       gap    pres   dres\n", " 0:  0.0000e+00  0.0000e+00  1e+00  1e+00  1e+00\n", " 1:  9.9000e-01  4.6251e+00  1e-02  2e+00  7e+01\n", " 2:  3.6389e+00  3.9677e+00  1e-04  1e-01  3e+01\n", " 3:  3.0555e+00  3.3406e+00  1e-06  1e-01  2e+01\n", " 4:  2.5112e+00  2.7758e+00  1e-08  1e-01  8e+00\n", " 5:  2.1118e+00  2.3358e+00  1e-10  1e-01  4e+00\n", " 6:  1.9684e+00  2.1118e+00  1e-12  6e-02  1e+00\n", " 7:  2.0493e+00  2.0796e+00  1e-14  1e-02  1e-01\n", " 8:  2.0790e+00  2.0794e+00  1e-16  2e-04  2e-03\n", " 9:  2.0794e+00  2.0794e+00  1e-18  2e-06  2e-05\n", "10:  2.0794e+00  2.0794e+00  1e-20  2e-08  2e-07\n", "11:  2.0794e+00  2.0794e+00  1e-22  2e-10  2e-09\n", "Optimal solution found.\n", "[ 5.00e-01]\n", "[ 2.50e-01]\n", "\n"]}], "source": ["from cvxopt import solvers, matrix, spdiag, log\n", "\n", "def acent(A, b):\n", "    m, n = A.size\n", "    def F(x=None, z=None):\n", "        if x is None: \n", "            # l'algorithme fonctionne de mani\u00e8re it\u00e9rative\n", "            # il faut choisir un x initial, c'est ce qu'on fait ici\n", "            return 0, matrix(1.0, (n,1))\n", "        if min(x) <= 0.0: \n", "            return None   # cas impossible\n", "        \n", "        # ici commence le code qui d\u00e9finit ce qu'est une it\u00e9ration\n", "        f = -sum(log(x))\n", "        Df = -(x**-1).T\n", "        if z is None: return f, Df\n", "        H = spdiag(z[0] * x**(-2))\n", "        return f, Df, H\n", "    \n", "    return solvers.cp(F, A=A, b=b)['x']\n", "\n", "A = matrix ( [[1.0,2.0]] ).T\n", "b = matrix ( [[ 1.0 ]] )\n", "print(acent(A,b))"]}, {"cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["[ 5.00e-01]\n", "[ 2.50e-01]\n", "\n"]}], "source": ["# il existe un moyen d'\u00e9viter l'affichage des logs (pratique si on doit faire\n", "# un grand nombre d'optimisation)\n", "from cvxopt import solvers\n", "solvers.options['show_progress'] = False\n", "print(acent(A,b))\n", "solvers.options['show_progress'] = True"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Cet exemple r\u00e9soud le probl\u00e8me de minimisation suivant :\n", "\n", "$\\left\\{ \\begin{array}{l} \\min_{X} \\left\\{ - \\sum_{i=1}^n \\ln x_i \\right \\} \\\\ sous \\; contrainte \\;  AX = b \\end{array} \\right.$\n", "\n", "Les deux modules [numpy](http://www.numpy.org/) et [cvxopt](http://cvxopt.org/) n'utilisent pas les m\u00eames matrices (les m\u00eames objets ``matrix``) bien qu'elles portent le m\u00eame 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\u00e9crite par une autre classe."]}, {"cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["[ 0.00e+00  4.50e+00]\n", "[ 1.50e+00 -6.00e+00]\n", "\n"]}], "source": ["import cvxopt\n", "m = cvxopt.matrix( [[ 0, 1.5], [ 4.5, -6] ] )\n", "print(m)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Ces informations devraient vous permettre de r\u00e9soudre le premier probl\u00e8me."]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Exercice 2 : l'algorithme de Arrow-Hurwicz\n", "\n", "On cherche \u00e0 r\u00e9soudre le m\u00eame probl\u00e8me avec l'algorithme de [Arrow-Hurwicz](http://dido.wss.yale.edu/P/ccdp/ec/e-0296.pdf) (tir\u00e9 de [Introduction \u00e0 l'optimisation](http://www.editions-ellipses.fr/product_info.php?products_id=8830) de Jean-Christophe Culioli, voir \u00e9galement l'[algorithme d'Uzawa](http://en.wikipedia.org/wiki/Uzawa_iteration)) et les notations suivantes :\n", "\n", "$\\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.$\n", "\n", "L'algorithme est le suivant :\n", "\n", "\n", "- On choisit $\\epsilon > 0$, $\\rho > 0$, des vecteurs al\u00e9atoires $U_0 \\in \\mathbb{R}^2$ et $P_0 \\in \\mathbb{R}^d$ ($d$ est le nombre de contraintes).\n", "- A l'it\u00e9ration $k$, on met \u00e0 jour :\n", "\n", "  $\\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}$\n", "\n", "- On retourne \u00e0 l'\u00e9tape pr\u00e9c\u00e9dente jusqu'\u00e0 ce que la suite $(U_k)$ n'\u00e9volue plus.\n", "\n", "\n", "\n", "Le coefficient $P_t$ correspond au coefficient de Lagrange pour un Lagrangien d\u00e9fini comme suit :\n", "\n", "$L(U,P) = J(U) + P' \\theta(U)$\n", "\n", "La suite $(U_t)_t$ converge vers la solution en se d\u00e9pla\u00e7ant le long du gradient de la fonction $J$ lorsque la contrainte est v\u00e9rifi\u00e9e. Lorsqu'elle ne l'est pas, cette direction est un m\u00e9lange du gradient de la fonction \u00e0 optimiser et de celui de la contrainte \u00e0 respecter.\n", "\n", "Impl\u00e9menter cet algorithme et v\u00e9rifier qu'il converge vers la m\u00eame solution que celle obtenue pour le premier probl\u00e8me."]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Exercice 3 : Le lagrangien augment\u00e9\n", "\n", "On reprend l'algorithme pr\u00e9c\u00e9dent mais appliqu\u00e9 \u00e0 la fonction :\n", "\n", "$A(U) = J(U) + \\frac{c}{2} \\theta^2(U)$\n", "\n", "L'objectif est ici de remplacer la fonction $J$ par la fonction $A$ (ou Lagrangien augment\u00e9) dans l'algorithme pr\u00e9c\u00e9dent et de v\u00e9rifier qu'il converge en un nombre moindre d'it\u00e9rations lorsqu'on choisit les m\u00eames conditions initiales."]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Prolongement 1 : in\u00e9galit\u00e9\n", "\n", "On souhaite maintenant ajouter une contrainte d'in\u00e9galit\u00e9 : $Du -d \\leqslant 0$. L'id\u00e9e est de cr\u00e9er une matrice $\\Theta$ d\u00e9finie par bloc :\n", "\n", "$\\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)$\n", "\n", "On consid\u00e8re que $\\Theta$ est une matrice. On \u00e9tend ce d\u00e9coupage en bloc au formule de mise \u00e0 jour de l'algorithme de Arrow-Hurwicz :\n", "\n", "$\\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}$\n", "\n", "Les coefficients $r_i$ valent 0 ou 1 selon que la contrainte d'in\u00e9galit\u00e9 $i$ est respect\u00e9e ou non. On applique la modification pour le probl\u00e8me :\n", "\n", "$\\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.$\n", "\n", "\n"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Prolongement 2 : optimisation d'une fonction lin\u00e9aire\n", "\n", "On s'int\u00e9resse \u00e0 un probl\u00e8me diff\u00e9rent :\n", "\n", "$\\left\\{ \\begin{array}{l} \\min_X C'X  \\\\ sous \\; contrainte \\;  AX \\leqslant B \\end{array}\\right.$\n", "\n", "Une fa\u00e7on de r\u00e9soudre ce probl\u00e8me est l'[algorithme de Karmarkar](http://fr.wikipedia.org/wiki/Algorithme\\_de\\_Karmarkar). On optimise maintenant une fonction lin\u00e9aire et non une fonction quadratique. Par cons\u00e9quent, la solution est forc\u00e9ment sur un bord : des contraintes sont satur\u00e9es.\n", "\n", "Vous pouvez d\u00e9couvrir un autre exemple d'optimisation dynamique au travers de ce blog : [Optimisation sous contraintes appliqu\u00e9e au calcul du report des voix](http://www.xavierdupre.fr/blog/2013-12-07_nojs.html)."]}, {"cell_type": "code", "execution_count": 6, "metadata": {"collapsed": true}, "outputs": [], "source": []}], "metadata": {"kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"}, "language_info": {"codemirror_mode": {"name": "ipython", "version": 3}, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.1"}}, "nbformat": 4, "nbformat_minor": 2}