06_Optimisation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimisation, `solve()` and Numerical Evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import numpy as np\n",
    "import sympy as sp\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Usually, when we are performing Lagrange optimisation, the problem is formulated as such:\n",
    "\n",
    "$$\n",
    "\\mathcal{L}(x,y,\\lambda) = f(x,y) + \\lambda \\ g(x,y)\n",
    "$$\n",
    "\n",
    "Here, $f(x,y)$ is the function we which to optimise _on_. Where as $g(x,y)$ is the function we want to apply as a _constraint_.\n",
    "\n",
    "A useful example is the function \n",
    "\n",
    "$f(x,y) = x^2 y^2 +5$\n",
    "\n",
    "We will optimise this for a maxima over a cirlce of radius $\\frac{1}{2}$ centred at $(1,1)$.\n",
    "\n",
    "This gives:\n",
    "\n",
    "$g(x,y) = (x-1)^2 + (y-1)^2 - 0.25$ \n",
    "\n",
    "## Graphical Optimisation\n",
    "\n",
    "## Question 1: \n",
    "\n",
    "Examine the code below. \n",
    "\n",
    "Using the visuals and adjusting values, attempt to guess values for the maximum of the function.\n",
    "\n",
    "_This approach can be useful when a problem is time-consuming (or impossible) to solve analytically and precision isn't paramount._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3D plot setup\n",
    "fig = plt.figure()\n",
    "ax = fig.gca(projection='3d')\n",
    "\n",
    "# Create boundaries of f(x,y)\n",
    "x = np.linspace(0,2)\n",
    "y = np.linspace(0,2)\n",
    "\n",
    "# Create 2D domain of f(x,y)\n",
    "xf, yf = np.meshgrid(x,y)\n",
    "\n",
    "# Discrete version of f(x,y) over 2D domain\n",
    "fxy = (xf**2)*(yf**2) + 5\n",
    "\n",
    "# Plot our function to be optimised f(x,y)\n",
    "ax.plot_surface(xf, yf, fxy, alpha=0.5)\n",
    "ax.view_init(20,-45)\n",
    "\n",
    "# Generate points for constraint function g(x,y)\n",
    "theta = np.linspace(0,2*np.pi)\n",
    "r = 0.5  # radius\n",
    "xg = r*np.cos(theta) + 1\n",
    "yg = r*np.sin(theta) + 1\n",
    "\n",
    "# Create g(x,y)\n",
    "ax.plot(xg, yg, color='g', alpha=0.5, label='$g(x,y)$')\n",
    "gxy = (xg**2)*(yg**2) + 5\n",
    "ax.plot(xg, yg, gxy, label='$f(x,y)$ evaluated over $g(x,y)$', color='r')\n",
    "\n",
    "\n",
    "#####################################\n",
    "# Your guess here\n",
    "xp = [1.250]\n",
    "yp = [1.45]\n",
    "zp = [(xp[0]**2)*(yp[0]**2) + 5]\n",
    "ax.scatter(xp, yp, zp, color='k', label='Your guess at constrained maximum')\n",
    "#####################################\n",
    "\n",
    "# Make it pretty\n",
    "ax.grid(False)\n",
    "ax.legend()\n",
    "plt.locator_params(nbins=5) # Amount of numbers per axis\n",
    "# Label our axes\n",
    "ax.set_xlabel(\"x\")\n",
    "ax.set_ylabel(\"y\")\n",
    "ax.set_zlabel(\"z\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Symbolic Optimisation\n",
    "\n",
    "This method makes use of methods covered in previous notebooks and the `sympy.solve()` function.\n",
    "\n",
    "The documentation for this new function can be found online for your particular version of `sympy`, but a pragmatic run-through will be given here.\n",
    "\n",
    "### `sympy.solve()`\n",
    "\n",
    "Equations (or even relational expressions) are passed to `solve()` in a list. If the expression is an equation, it should be expressed such that it is equal to zero.\n",
    "\n",
    "### One Equation \n",
    "\n",
    "e.g. Let's solve\n",
    "\n",
    "$$\n",
    "x^2 = 64\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import solve, Symbol\n",
    "x = Symbol('x')\n",
    "f = x**2 - 6  # Notice that the expression is set up to be equal to zero.\n",
    "\n",
    "sols = solve(f,x)\n",
    "display(sols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### More than One Equation\n",
    "\n",
    "Take the system of equations:\n",
    "\n",
    "$$\n",
    "x + 2y -z = 5\n",
    "$$\n",
    "\n",
    "$$\n",
    "2x - y -z = 0\n",
    "$$\n",
    "\n",
    "$$\n",
    "-x -y +3z = 8\n",
    "$$\n",
    "\n",
    "We have to express these as equal to zero, so:\n",
    "\n",
    "$$\n",
    "x + 2y -z - 5 = 0\n",
    "$$\n",
    "\n",
    "$$\n",
    "2x - y -z = 0\n",
    "$$\n",
    "\n",
    "$$\n",
    "-x -y +3z -8 = 0\n",
    "$$\n",
    "\n",
    "and then we can use them in `sympy.solve()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Adding our new symbols\n",
    "y = Symbol('y')\n",
    "z = Symbol('z')\n",
    "\n",
    "# All equations equal to zero\n",
    "f1 = x + 2*y - z - 5\n",
    "f2 = 2*x - y - z \n",
    "f3 = -x -y +3*z - 8\n",
    "\n",
    "sols = solve([f1,f2,f3], (x,y,z))  # Returns a dictionary (dict)\n",
    "\n",
    "display(sols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the same equation used in the 'Matrices' notebook -- first looked at when inverting matrices to solve sets of linear equations. These two methods can be looked at as an example of two different ways to obtain the same solution. You could also theoretically create a matrix version of Lagrange Optimisation (however the logic for generating the matrix from partial derivatives requires some time)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Actually Doing Symbolic Optimisation in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at:\n",
    "\n",
    "$$\n",
    "f(x,y) = x^2 + y^2\n",
    "$$\n",
    "\n",
    "Constrained by:\n",
    "\n",
    "$$\n",
    "xy = 1\n",
    "$$\n",
    "\n",
    "Meaning that:\n",
    "\n",
    "$$ \n",
    "g(x,y) = xy - 1\n",
    "$$\n",
    "\n",
    "Hence:\n",
    "\n",
    "$$\n",
    "\\mathcal{L} = f(x,y) + \\lambda \\ g(x,y)\n",
    "$$\n",
    "\n",
    "where $\\lambda$ is our unknown Lagrange multiplier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating our lagrangian\n",
    "lmbda = sp.Symbol('λ')\n",
    "sp.init_printing()\n",
    "x = sp.Symbol('x')\n",
    "y = sp.Symbol('y')\n",
    "f = (x**2)+(y**2)\n",
    "g = x*y - 1\n",
    "\n",
    "display('Our Lagrangian')\n",
    "L  = f + lmbda*g\n",
    "display(L)\n",
    "\n",
    "display('Partial Derivatives')\n",
    "Lx = sp.diff(L,x)\n",
    "Ly = sp.diff(L,y)\n",
    "Ll = sp.diff(L, lmbda)\n",
    "display(Lx)\n",
    "display(Ly)\n",
    "display(Ll)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "extrema = sp.solve([Lx, Ly, Ll], (x,y,lmbda))\n",
    "\n",
    "for i in range(len(extrema)):\n",
    "    display(extrema[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we are examining the function in the real domain ($\\mathbb{R}$) we can ignore the solutions that contain\n",
    "$i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Question 2: \n",
    "\n",
    "Create two plots.\n",
    "\n",
    "## (a)\n",
    "\n",
    "Using `matplotlib` (`plt`) to display the function $f(x,y)$, $g(x,y)$ as well as the **real** extrema obtained by symbolic calculation above.\n",
    "\n",
    "## (b)\n",
    "\n",
    "A contour plot (using `plt.contour()`) of $f(x,y)$ and the constrained maxima points obtained by symbolic calculation. \n",
    "\n",
    "Note: If you do not recall how to create a contour plot, look at the `plt` documentation for `contour()` or the other notebook titled 'Partial Derivatives'.\n",
    "\n",
    "# Answer 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3D plot setup\n",
    "fig = plt.figure()\n",
    "ax = fig.gca(projection='3d')\n",
    "\n",
    "# Create boundaries of f(x,y)\n",
    "x = np.linspace(-2,2)\n",
    "y = np.linspace(-2,2)\n",
    "\n",
    "# Create 2D domain of f(x,y)\n",
    "xf, yf = np.meshgrid(x,y)\n",
    "\n",
    "# Discrete version of f(x,y) over 2D domain\n",
    "fxy = (xf**2)+(yf**2) \n",
    "\n",
    "# Plot our function to be optimised f(x,y)\n",
    "ax.plot_surface(xf, yf, fxy, alpha=0.1)\n",
    "ax.view_init(15, 20)\n",
    "\n",
    "# Generate points for constraint function g(x,y)\n",
    "xg = x\n",
    "xg1 = np.linspace(-2, -0.1)\n",
    "xg2 = np.linspace(0.1, 2)\n",
    "yg1 = 1/xg1\n",
    "yg2 = 1/xg2\n",
    "ax.plot(xg1, yg1, color='g', label='$g(x,y)$')\n",
    "ax.plot(xg2, yg2, color='g')\n",
    "\n",
    "# Plot extrema\n",
    "\n",
    "x_extrema = np.array([-1,1])\n",
    "y_extrema = np.array([-1,1])\n",
    "z_extrema = (x_extrema**2)+(y_extrema**2) \n",
    "\n",
    "ax.scatter(x_extrema, y_extrema, z_extrema,  color='r', label='maxima')\n",
    "\n",
    "# Make it pretty\n",
    "ax.grid(False)\n",
    "ax.legend()\n",
    "plt.locator_params(nbins=5) # Amount of numbers per axis\n",
    "# Label our axes\n",
    "ax.set_xlabel(\"x\")\n",
    "ax.set_ylabel(\"y\")\n",
    "ax.set_zlabel(\"z\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Question 3:\n",
    "\n",
    "_(Using what you now have seen about the symbolic method of solving Lagrange Optimisation problems)_\n",
    "\n",
    "For the function \n",
    "\n",
    "$$\n",
    "f(x,y) = x^2 y \n",
    "$$\n",
    "\n",
    "subject to the constraint \n",
    "\n",
    "$$ \n",
    "x^2 + 2y^2 = 6\n",
    "$$\n",
    "\n",
    "find the lagrange multipliers and the extrema.\n",
    "\n",
    "## Answer 3:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluating Partial Derivatives at Points and Along Curves\n",
    "\n",
    "Say we have some function $f(x(t),y(t))$. Then the _total derivative_ with respect to $t$ would be:\n",
    "\n",
    "$$\n",
    "\\frac{df}{dt} = \n",
    "\\frac{\\partial f}{\\partial x} \\frac{d x}{dt} +\n",
    "\\frac{\\partial f}{\\partial y} \\frac{d y}{dt} \n",
    "$$\n",
    "\n",
    "Let's calculate the total derivative of:\n",
    "\n",
    "$$\n",
    "f(x,y) = 2x^2 + y^2\n",
    "$$\n",
    "\n",
    "along the parametric path:\n",
    "\n",
    "$$\n",
    "\\langle x(t), y(t) \\rangle =\n",
    "\\langle \\cos(t), \\sin(t) \\rangle\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy.plotting import plot3d, plot_parametric\n",
    "t = Symbol('t')\n",
    "\n",
    "# Function and its partial derivatives.\n",
    "f = 2*x**2 + y**2 \n",
    "fx = sp.diff(f,x)\n",
    "fy = sp.diff(f,y)\n",
    "\n",
    "# Contour to evaluate df/dt over\n",
    "xc = sp.cos(t)\n",
    "yc = sp.sin(t)\n",
    "\n",
    "# Substitution\n",
    "fxc = fx.subs(x, xc)\n",
    "fyc = fy.subs(y,yc)\n",
    "\n",
    "# dx/dt and dy/dt of path in xy plane.\n",
    "dxc_dt = sp.diff(xc, t)\n",
    "dyc_dt = sp.diff(yc, t)\n",
    "\n",
    "\n",
    "display('Total Derivative')\n",
    "dfdt = fxc*dxc_dt+fyc*dyc_dt\n",
    "display(dfdt)\n",
    "\n",
    "\n",
    "display('Plot of f(x,y)')\n",
    "plot3d(f)\n",
    "\n",
    "\n",
    "display('The path we are taking the derivative along in xy plane.')\n",
    "plot_parametric(xc,yc )\n",
    "\n",
    "\n",
    "display('Total derivative df/dt along curve <x(t), y(t)>')\n",
    "sp.plot(dfdt)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerically Evaluating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Should you need to evaluate certain points, (i.e. a numerical value of the gradient is required) you can do so using a combination of methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.N(dfdt.subs(t,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the `N()` function means to numerically calculate a function. Before we did that though, $t$ was substituted for 2 with the `subs()` method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 4:\n",
    "\n",
    "_Using what we now know about total derivatives and how to calculate them in Python_\n",
    "\n",
    "Given:\n",
    "\n",
    "$$\n",
    "f(x,y,z) = x^2 + y^2 + z^2\n",
    "$$\n",
    "\n",
    "compute the symbolic expression for total derivative $\\frac{df}{dt}$ along the path:\n",
    "\n",
    "$$\n",
    "\\langle 2 e^t, 2 e^{-t}, 1 + t^2 \\rangle\n",
    "$$\n",
    "\n",
    "then numerically evaluate it at the point $t=3$.\n",
    "\n",
    "## Answer 4:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}