{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Maximisation de l'utilité sous contrainte\n",
    "\n",
    "_F-X. Dehon - dehon[@]unice.fr - 14 avr 2020_\n",
    "\n",
    "Utilité $U(x,y)=(x+2)(x+3y)$ sous les contraintes $x\\geq 0$, $y\\geq 0$ et $5x+3y\\leq r$.\n",
    "\n",
    "On dessine :\n",
    "\n",
    "- en orange le graphe de la fonction utilité : l'ensemble des points $(x\\geq 0,y\\geq 0,z)$ tels que $z=U(x,y)$\n",
    "- en vert le plan donné par la contrainte (saturée) de revenu disponible et du coût des biens : $5x+3y=r(=20)$\n",
    "- en bleu clair le plan $z=u (=42)$.\n",
    "\n",
    "On visualise ainsi\n",
    "\n",
    "- les $(x,y)$ tels que $U(x,y)$ soit maximal sous les contraintes $x\\geq 0$, $y\\geq 0$ et $5x+3y\\leq r$,\n",
    "- le graphe de la restriction de $U$ à l'ensemble des $(x,y)$ tels que $x\\geq 0$, $y\\geq 0$ et $5x+3y=r$, de même que ceux des restrictions de $U$ aux $(x,0),\\ x\\geq 0$ et aux $(0,y),\\ y\\geq 0$,\n",
    "- la ligne de niveau $U(x,y)=u (=42)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "##### 1. Graphe de la fonction d'utilité"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "r=20;u=42\n",
    "x,y,z = var('x,y,z')\n",
    "xmax,ymax,zmax=(8,8,80)\n",
    "\n",
    "b(x,y)=r-(5*x+3*y) #contrainte revenu - prix\n",
    "U(x,y)=(x+2)*(x+3*y) #utilité"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "A=line3d([(0,0,0),(xmax,0,0),(0,0,0)])+line3d([(0,0,0),(0,ymax,0),(0,0,0)])+line3d([(0,0,0),(0,0,zmax),(0,0,0)])#bug avec line3d et threejs\n",
    "T=text3d(' x ', (xmax/2, -.5, -.5))+text3d(' y ', (-.5, ymax/2, -.5))+text3d(' z ', (-.5, -.5,zmax/2))+text3d('0',(-.5, -.5,-.5))+text3d(str(xmax),(xmax, -.5,-.5))+text3d(str(ymax),(-.5,ymax,-.5))+text3d(str(zmax),(-.5,-.5,zmax))\n",
    "B=implicit_plot3d(b,(x,0,r/5),(y,0,r/3),(z,0,zmax),color='green', opacity=0.8)\n",
    "I=implicit_plot3d(z==u,(x,0,xmax),(y,0,ymax),(z,0,zmax),color='lightblue', opacity=0.5)+text3d('z='+str(u),(-.5,-.5,u))\n",
    "P = implicit_plot3d(z-U(x,y),(x,0,xmax),(y,0,ymax),(z,0,zmax),color='orange', opacity=0.7)#,region=b\n",
    "\n",
    "show(P+B+A+T+I,frame=false,aspect_ratio=[1,1,.1],viewer='threejs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "##### 2. Lignes de niveau (appelées courbes d'indifférence en microéconomie) de la fonction d'utilité $U$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "JJ=contour_plot(U(x,y),(x,0,xmax),(y,0,ymax),contours=15,fill=False,labels=True,label_inline=False)\n",
    "show(JJ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Courbe d'indifférence $U(x,y)=42$ (en orange), région $\\big( x\\geq 0\\textrm{ et }y\\geq 0\\textrm{ et }5x+3y\\leq 20\\big) $ en vert, $1/10$ème du gradiant de $U$ en $(2,y)$ où $U(2,y)=42$ :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "var('y')\n",
    "x0=2;y0=solve(U(x0,y)==42,y)[0].rhs()\n",
    "print y0, U(x0,y0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "var('x,y')\n",
    "R=region_plot([b(x,y)>=0,x>=0,y>=0],(x,-.1,xmax),(y,-.1,ymax),incol='lightgreen',bordercol='green')+implicit_plot(U(x,y)==u,(x,0,xmax),(y,0,ymax),color='orange')\n",
    "G=arrow((x0,y0),vector((x0,y0))+1/10*U.gradient()(x0,y0),color='orange')\n",
    "show(R+G,aspect_ratio=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "##### 3. La courbe d'indifférence portée sur le graphe de la fonction d'utilité $U$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "N=implicit_plot(U(x,y)==u,(x,0,xmax),(y,0,ymax)).matplotlib().get_children()[1].collections[0].get_paths()[0].to_polygons(closed_only=False)[0]\n",
    "N3d=line3d([[p[0],p[1],u] for p in N],color='red')+text3d('z='+str(u),(-.5,-.5,u))\n",
    "N3d0=line3d([[p[0],p[1],0] for p in N],color='black')\n",
    "show(P+B+A+T+I+N3d+N3d0+line3d([[r/5,0,0],[0,r/3,0],[r/5,0,0]],color='black'),frame=false,aspect_ratio=[1,1,.1],viewer='threejs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "##### 4. Documentation Sagemath\n",
    "\n",
    "https://ask.sagemath.org/question/9262/get-list-of-coordinates-from-plot-object/\n",
    "\n",
    "https://matplotlib.org/3.1.0/api/path_api.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "N3d=[[p[0],p[1],u] for p in N]#;print N3d\n",
    "show(line3d(N3d,color='black'),aspect_ratio=[1,1,.1],viewer='threejs')\n",
    "\n",
    "#I=implicit_plot3d((z-1)^2+(x-1)^2,(x,0,xmax),(y,0,ymax),(z,0,zmax),color='black')#n'est pas dessiné\n",
    "#I=implicit_plot(U(x,y)==u,(x,0,xmax),(y,0,ymax),color='black').plot3d()#produit une erreur\n",
    "#NN=line2d(N,color='black').plot3d(u)+text3d('z='+str(u),(-.5,-.5,u))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "##### 5. Max de l'utilité sous contrainte : fonction Sagemath native"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "minimize_constrained?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a = (1.500057102546868, 4.16657149575522) , b(a) = 0.0 , U(a) = 48.9999999869572\n"
     ]
    }
   ],
   "source": [
    "var('x,y')\n",
    "bb=lambda x,y:r-(5*x+3*y)\n",
    "#def bb(x,y):return(r-(5*x+3*y)\n",
    "f=lambda (x,y):-U(x,y)\n",
    "c0=lambda (x,y):bb(x,y)#contrainte non prise en compte si b(x,y) au lieu de bb(x,y) !\n",
    "c1=lambda (x,y):x\n",
    "c2=lambda (x,y):y\n",
    "a=minimize_constrained(f,[c1,c2,c0],(1,2))\n",
    "print 'a =',a,', b(a) =',c0(a),', U(a) =',U(*a)#comparer avec ce qu'on voit sur le dessin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath (default)",
   "language": "sagemath",
   "metadata": {
    "cocalc": {
     "description": "Open-source mathematical software system",
     "priority": 10,
     "url": "https://www.sagemath.org/"
    }
   },
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}