{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["# 1A.algo - Optimisation sous contrainte (correction)\n", "\n", "Un peu plus de d\u00e9tails 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": ["Plan\n", "
run previous cell, wait for 2 seconds
\n", ""], "text/plain": [""]}, "execution_count": 2, "metadata": {}, "output_type": "execute_result"}], "source": ["from jyquickhelper import add_notebook_menu\n", "add_notebook_menu()"]}, {"cell_type": "markdown", "metadata": {}, "source": ["On rappelle le probl\u00e8me d'optimisation \u00e0 r\u00e9soudre :\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.5 \\end{array}\\right .$\n", "\n", "Les impl\u00e9mentations de l'algorithme Arrow-Hurwicz propos\u00e9es ici ne sont pas g\u00e9n\u00e9riques. Il n'est pas sugg\u00e9r\u00e9 de les r\u00e9utiliser \u00e0 moins d'utiliser pleinement le calcul matriciel de [numpy](http://www.numpy.org/)."]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Exercice 1 : optimisation avec cvxopt\n", "\n", "Le module [cvxopt](http://cvxopt.org/) utilise une fonction qui retourne la valeur de la fonction \u00e0 optimiser, sa d\u00e9riv\u00e9e, sa d\u00e9riv\u00e9e seconde.\n", "\n", "$\\begin{array}{rcl} f(x,y) &=& x^2 + y^2 - xy + y \\\\ \\frac{\\partial f(x,y)}{\\partial x} &=& 2x - y \\\\ \\frac{\\partial f(x,y)}{\\partial y} &=& 2y - x + 1 \\\\ \\frac{\\partial^2 f(x,y)}{\\partial x^2} &=& 2 \\\\ \\frac{\\partial^2 f(x,y)}{\\partial y^2} &=& 2 \\\\ \\frac{\\partial^2 f(x,y)}{\\partial x\\partial y} &=& -1 \\end{array}$"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Le param\u00e8tre le plus complexe est la fonction ``F`` pour lequel il faut lire la documentation de la fonction [solvers.cp](http://cvxopt.org/userguide/solvers.html#problems-with-nonlinear-objectives) qui d\u00e9taille les trois cas d'utilisation de la fonction ``F`` :\n", "\n", "* ``F()`` ou ``F(None,None)``, ce premier cas est sans doute le plus d\u00e9routant puisqu'il faut retourner le nombre de contraintes non lin\u00e9aires et le premier $x_0$\n", "* ``F(x)`` ou ``F(x,None)``\n", "* ``F(x,z)``\n", "\n", "L'algorithme de r\u00e9solution est it\u00e9ratif : on part d'une point $x_0$ qu'on d\u00e9place dans les directions oppos\u00e9s aux gradients de la fonction \u00e0 minimiser et des contraintes jusqu'\u00e0 ce que le point $x_t$ n'\u00e9volue plus. C'est pourquoi le premier d'utilisation de la focntion $F$ est en fait une initialisation. L'algorithme d'optimisation a besoin d'un premier point $x_0$ dans le domaine de d\u00e9fintion de la fonction $f$."]}, {"cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": [" pcost dcost gap pres dres\n", " 0: 0.0000e+00 4.3222e-01 1e+00 1e+00 1e+00\n", " 1: 2.6022e-01 4.2857e-01 1e-02 1e-01 1e-02\n", " 2: 4.2687e-01 4.2857e-01 1e-04 1e-03 1e-04\n", " 3: 4.2855e-01 4.2857e-01 1e-06 1e-05 1e-06\n", " 4: 4.2857e-01 4.2857e-01 1e-08 1e-07 1e-08\n", " 5: 4.2857e-01 4.2857e-01 1e-10 1e-09 1e-10\n", "Optimal solution found.\n", "{'dual objective': 0.42857142857142855, 'primal objective': 0.4285714268720223, 'primal slack': 1.000000000000004e-10, 'snl': <0x1 matrix, tc='d'>, 'relative gap': 2.333333333333341e-10, 'sl': <0x1 matrix, tc='d'>, 'dual slack': 0.9999999999999991, 'status': 'optimal', 'y': <1x1 matrix, tc='d'>, 'x': <2x1 matrix, tc='d'>, 'zl': <0x1 matrix, tc='d'>, 'znl': <0x1 matrix, tc='d'>, 'dual infeasibility': 9.995026102717158e-11, 'gap': 1.0000000000000031e-10, 'primal infeasibility': 1.2214984990086318e-09}\n", "solution: [ 4.29e-01 2.86e-01]\n", "\n"]}], "source": ["from cvxopt import solvers, matrix\n", "import random\n", "\n", "def fonction(x=None,z=None) : \n", " if x is None :\n", " x0 = matrix ( [[ random.random(), random.random() ]])\n", " return 0,x0\n", " f = x[0]**2 + x[1]**2 - x[0]*x[1] + x[1]\n", " d = matrix ( [ x[0]*2 - x[1], x[1]*2 - x[0] + 1 ] ).T\n", " if z is None: \n", " return f, d\n", " else : \n", " h = z[0] * matrix ( [ [ 2.0, -1.0], [-1.0, 2.0] ])\n", " return f, d, h\n", " \n", "A = matrix([ [ 1.0, 2.0 ] ]).trans()\n", "b = matrix ( [[ 1.0] ] )\n", " \n", "sol = solvers.cp ( fonction, A = A, b = b)\n", "print (sol)\n", "print (\"solution:\",sol['x'].T)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Exercice 2 : l'algorithme de Arrow-Hurwicz"]}, {"cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["i 100 diff 0.001 : 0.4247897864911925 [0.42689355878508073, 0.28435244828471445] -0.5749941851846841 -0.004401544645490363\n", "i 200 diff 4e-06 : 0.42858505488192533 [0.4285774472613499, 0.285720126646874] -0.5714194393546865 1.7700555097865944e-05\n", "i 300 diff 1e-08 : 0.42857138230973435 [0.4285714081852857, 0.28571426356925744] -0.5714285933247762 -6.467619939609648e-08\n", "8.65032490083e-11 353 -0.5714285707449258\n", "solution: [0.42857142760818157, 0.28571428421894984]\n"]}], "source": ["def fonction(X) : \n", " x,y = X\n", " f = x**2 + y**2 - x*y + y\n", " d = [ x*2 - y, y*2 - x + 1 ] \n", " return f, d\n", " \n", "def contrainte(X) : \n", " x,y = X\n", " f = x+2*y-1\n", " d = [ 1,2]\n", " return f, d\n", " \n", "X0 = [ random.random(),random.random() ]\n", "p0 = random.random()\n", "epsilon = 0.1\n", "rho = 0.1\n", "\n", "diff = 1\n", "iter = 0\n", "while diff > 1e-10 :\n", " f,d = fonction( X0 )\n", " th,dt = contrainte( X0 )\n", " Xt = [ X0[i] - epsilon*(d[i] + dt[i] * p0) for i in range(len(X0)) ]\n", "\n", " th,dt = contrainte( Xt )\n", " pt = p0 + rho * th\n", " \n", " iter += 1\n", " diff = sum ( [ abs(Xt[i] - X0[i]) for i in range(len(X0)) ] )\n", " X0 = Xt\n", " p0 = pt\n", " if iter % 100 == 0 :\n", " print (\"i {0} diff {1:0.000}\".format(iter,diff),\":\", f,X0,p0,th)\n", "\n", "print (diff,iter,p0)\n", "print(\"solution:\",X0)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["La code propos\u00e9 ici a \u00e9t\u00e9 repris et modifi\u00e9 de fa\u00e7on \u00e0 l'inclure dans une fonction qui s'adapte \u00e0 n'importe quel type de fonction et contrainte d\u00e9rivables : [Arrow_Hurwicz](http://www.xavierdupre.fr/app/ensae_teaching_cs/helpsphinx/td_1a/optimisation_contrainte.html?highlight=arrow#ensae_teaching_cs.td_1a.optimisation_contrainte.Arrow_Hurwicz). Il faut distinguer l'algorithme en lui-m\u00eame et la preuve de sa convergence. Cet algorithme fonctionne sur une grande classe de fonctions mais sa convergence n'est assur\u00e9e que lorsque les fonctions sont quadratiques."]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Exercice 3 : le lagrangien augment\u00e9"]}, {"cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["i 100 diff 9e-06 : 0.4284835397042614 [0.4285257942646543, 0.28566702773052666] -0.5712849940942254 -0.00014015027429237215\n", "i 200 diff 8e-10 : 0.42857142064500353 [0.4285714244564873, 0.2857142814529341] -0.5714285584819018 -1.263764448644622e-08\n", "9.59901602648e-11 223 -0.5714285699086377\n", "solution: [0.42857142808833615, 0.28571428521400477]\n"]}], "source": ["def fonction(X,c) : \n", " x,y = X\n", " f = x**2 + y**2 - x*y + y\n", " d = [ x*2 - y, y*2 - x + 1 ] \n", " \n", " v = x+2*y-1\n", " v = c/2 * v**2\n", " \n", " # la fonction retourne maintenant dv (ce qu'elle ne faisait pas avant)\n", " dv = [ 2*(x+2*y-1), 4*(x+2*y-1) ]\n", " dv = [ c/2 * dv[0], c/2 * dv[1] ]\n", " return f + v, d, dv\n", " \n", "def contrainte(X) : \n", " x,y = X\n", " f = x+2*y-1\n", " d = [ 1,2]\n", " return f, d\n", " \n", "X0 = [ random.random(),random.random() ]\n", "p0 = random.random()\n", "epsilon = 0.1\n", "rho = 0.1\n", "c = 1\n", "\n", "diff = 1\n", "iter = 0\n", "while diff > 1e-10 :\n", " f,d,dv = fonction( X0,c )\n", " th,dt = contrainte( X0 )\n", " # le dv[i] est nouveau\n", " Xt = [ X0[i] - epsilon*(d[i] + dt[i] * p0 + dv[i]) for i in range(len(X0)) ]\n", "\n", " th,dt = contrainte( Xt )\n", " pt = p0 + rho * th\n", " \n", " iter += 1\n", " diff = sum ( [ abs(Xt[i] - X0[i]) for i in range(len(X0)) ] )\n", " X0 = Xt\n", " p0 = pt\n", " if iter % 100 == 0 :\n", " print (\"i {0} diff {1:0.000}\".format(iter,diff),\":\", f,X0,p0,th)\n", " \n", "print (diff,iter,p0)\n", "print(\"solution:\",X0)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Prolongement 1 : in\u00e9galit\u00e9\n", "\n", "Le probl\u00e8me \u00e0 r\u00e9soudre est le suivant : \n", "\n", "$\\left\\{ \\begin{array}{l} \\min_U J(U) = u_1^2 + u_1^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.$"]}, {"cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": [" pcost dcost gap pres dres\n", " 0: 0.0000e+00 5.1249e-01 2e+00 1e+00 9e-01\n", " 1: 3.9825e-01 4.1007e-01 7e-02 4e-02 9e-03\n", " 2: 4.3458e-01 4.2532e-01 1e-02 3e-03 2e-16\n", " 3: 4.3118e-01 4.2927e-01 4e-03 9e-04 4e-16\n", " 4: 4.3033e-01 4.2993e-01 8e-04 2e-04 5e-16\n", " 5: 4.3005e-01 4.3000e-01 1e-04 2e-05 6e-16\n", " 6: 4.3000e-01 4.3000e-01 3e-06 9e-07 1e-15\n", " 7: 4.3000e-01 4.3000e-01 4e-08 1e-08 1e-15\n", "Optimal solution found.\n", "{'status': 'optimal', 'primal infeasibility': 9.636674935018775e-09, 'y': <1x1 matrix, tc='d'>, 'snl': <0x1 matrix, tc='d'>, 'sl': <1x1 matrix, tc='d'>, 'primal slack': 8.714188716285685e-11, 'dual slack': 0.20000259350365998, 'relative gap': 8.636668517465485e-08, 'znl': <0x1 matrix, tc='d'>, 'primal objective': 0.43000001865573895, 'x': <2x1 matrix, tc='d'>, 'gap': 3.713767462508084e-08, 'dual objective': 0.4299999999997598, 'dual infeasibility': 1.1909050386352223e-15, 'zl': <1x1 matrix, tc='d'>}\n", "solution: [ 4.00e-01 3.00e-01]\n", "\n"]}], "source": ["from cvxopt import solvers, matrix\n", "import random\n", "\n", "def fonction(x=None,z=None) : \n", " if x is None :\n", " x0 = matrix ( [[ random.random(), random.random() ]])\n", " return 0,x0\n", " f = x[0]**2 + x[1]**2 - x[0]*x[1] + x[1]\n", " d = matrix ( [ x[0]*2 - x[1], x[1]*2 - x[0] + 1 ] ).T\n", " h = matrix ( [ [ 2.0, -1.0], [-1.0, 2.0] ])\n", " if z is None: return f, d\n", " else : return f, d, h\n", " \n", "A = matrix([ [ 1.0, 2.0 ] ]).trans()\n", "b = matrix ( [[ 1.0] ] )\n", "\n", "G = matrix ( [[0.0, -1.0] ]).trans()\n", "h = matrix ( [[ -0.3] ] )\n", " \n", "sol = solvers.cp ( fonction, A = A, b = b, G=G, h=h)\n", "print (sol)\n", "print (\"solution:\",sol['x'].T) "]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Version avec l'algorithme de Arrow-Hurwicz"]}, {"cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["iteration 100 0.438928769823\n", "iteration 200 0.431511185295\n", "iteration 300 0.430474413611\n", "iteration 400 0.430161042987\n", "iteration 500 0.430060328795\n", "iteration 600 0.430025004165\n", "iteration 700 0.430011288151\n", "iteration 800 0.430005417575\n", "iteration 900 0.430002702711\n", "988\n", "solution: [[ 0.39982692 0.30007332]]\n"]}], "source": ["import numpy,random\n", "\n", "X0 = numpy.matrix ( [[ random.random(), random.random() ]]).transpose()\n", "P0 = numpy.matrix ( [[ random.random(), random.random() ]]).transpose()\n", " \n", "A = numpy.matrix([ [ 1.0, 2.0 ], [ 0.0, -1.0] ])\n", "tA = A.transpose()\n", "b = numpy.matrix ( [[ 1.0], [-0.30] ] )\n", "\n", "epsilon = 0.1\n", "rho = 0.1\n", "c = 1 \n", "\n", "first = True\n", "iter = 0\n", "while first or abs(J - oldJ) > 1e-8 :\n", " if first :\n", " J = X0[0,0]**2 + X0[1,0]**2 - X0[0,0]*X0[1,0] + X0[1,0]\n", " oldJ = J+1\n", " first = False\n", " else :\n", " oldJ = J\n", " J = X0[0,0]**2 + X0[1,0]**2 - X0[0,0]*X0[1,0] + X0[1,0]\n", " \n", " dj = numpy.matrix ( [ X0[0,0]*2 - X0[1,0], X0[1,0]*2 - X0[0,0] + 1 ] ).transpose()\n", "\n", " Xt = X0 - ( dj + tA * P0 ) * epsilon\n", " Pt = P0 + ( A * Xt - b) * rho\n", " \n", " if Pt [1,0] < 0 : Pt[1,0] = 0\n", " \n", " X0,P0 = Xt,Pt\n", " iter += 1\n", " if iter % 100 == 0 :\n", " print (\"iteration\",iter, J)\n", " \n", "print (iter)\n", "print (\"solution:\",Xt.T)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Prolongement 2 : optimisation d'une fonction lin\u00e9aire\n", "\n", "Correction \u00e0 venir."]}, {"cell_type": "code", "execution_count": 7, "metadata": {}, "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}