{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Penalty Methods\n",
    "For this training course we'll be using the optimization routines in SciPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.optimize as opt\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need to define a callback function so that we can plot the optimization steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def callback(x):\n",
    "    global xprev\n",
    "    plt.plot([xprev[0],x[0]],[xprev[1],x[1]],'b.-')\n",
    "    xprev = x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.1 Quadratic Penalty Method\n",
    "\n",
    "Recall that the quadratic penalty function, for a single equality constraint $c$, is given by\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "\\Phi_{\\mu}(x) = f(x) + \\frac{1}{2\\mu} \\left(c(x)\\right)^2\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "where $\\mu >0$ is a penalty parameter. The gradient for the quadratic penalty is\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "\\nabla \\Phi_{\\mu}(x) = \\nabla f(x) + \\frac{1}{\\mu}c(x)\\nabla c(x)\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "and the Hessian is given by\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "\\nabla^2 \\Phi_{\\mu}(x) = \\nabla^2 f(x) + \\frac{1}{\\mu}c(x)\\nabla^2 c(x)  + \\frac{1}{\\mu}\\nabla c(x)\\nabla c(x)^T\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "#### Coding Task:\n",
    "\n",
    "For different fixed values of $\\mu$, apply your favourite SciPy optimization routine to minimize $\\Phi_{\\mu}(x)$ for the following objectives:\n",
    "\n",
    "1. $f(x) = \\frac{1}{2}(ax_1^2 + x_2^2) \\qquad\\qquad\\qquad$ with $\\; c(x) = x_1 + x_2^2 - 1$ \n",
    "\n",
    "2. $f(x) = (a - x_1)^2 + b(x_2 - x_1^2)^2 \\qquad$ with $\\; c(x) = x_1 + x_2^2 - 1$ \n",
    "\n",
    "starting at $(x_1,x_2) = (1,1)$ for suitable choices of parameters $a,b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "a = 1\n",
    "mu = \n",
    "\n",
    "# Initial guess\n",
    "x0 = np.array(...)\n",
    "\n",
    "# Objective function, gradient and Hessian\n",
    "obj = lambda x: \n",
    "gobj = lambda x: \n",
    "hobj = lambda x: \n",
    "\n",
    "# Constraint function, gradient and Hessian\n",
    "con = lambda x: \n",
    "gcon = lambda x: \n",
    "hcon = lambda x: \n",
    "\n",
    "# Penalty function, gradient and Hessian\n",
    "fun = lambda x: \n",
    "jac = lambda x: \n",
    "hess = lambda x: \n",
    "\n",
    "# Plot function contours\n",
    "plt.figure()\n",
    "X = np.linspace(...)\n",
    "Y = np.linspace(...)\n",
    "Z = np.meshgrid(X,Y)\n",
    "plt.contour(...)\n",
    "plt.colorbar()\n",
    "\n",
    "# Plot constraint\n",
    "plt.contour(X,Y,con(Z),[0],colors='r')\n",
    "\n",
    "# Call SciPy's Newton-CG\n",
    "xprev = x0 # for plotting\n",
    "res = opt.minimize(..., method=..., callback=callback)\n",
    "\n",
    "# Print results and show plot\n",
    "print(res)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "a = 1\n",
    "b = 100\n",
    "mu = \n",
    "\n",
    "# Initial guess\n",
    "x0 = np.array(...)\n",
    "\n",
    "# Objective function, gradient and Hessian\n",
    "obj = lambda x: \n",
    "gobj = lambda x: \n",
    "hobj = lambda x: \n",
    "\n",
    "# Constraint function, gradient and Hessian\n",
    "con = lambda x: \n",
    "gcon = lambda x: \n",
    "hcon = lambda x: \n",
    "\n",
    "# Penalty function, gradient and Hessian\n",
    "fun = lambda x: \n",
    "jac = lambda x: \n",
    "hess = lambda x: \n",
    "\n",
    "# Plot function contours\n",
    "plt.figure()\n",
    "X = np.linspace(...)\n",
    "Y = np.linspace(...)\n",
    "Z = np.meshgrid(X,Y)\n",
    "plt.contour(...)\n",
    "plt.colorbar()\n",
    "\n",
    "# Plot constraint\n",
    "plt.contour(X,Y,con(Z),[0],colors='r')\n",
    "\n",
    "# Call SciPy's Newton-CG\n",
    "xprev = x0 # for plotting\n",
    "res = opt.minimize(..., method=..., callback=callback)\n",
    "\n",
    "# Print results and show plot\n",
    "print(res)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.2 Augmented Lagrangian Method\n",
    "\n",
    "Recall that the augmented Lagrangian function, for a single equality constraint $c$, is given by\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "\\Psi_{u,\\mu}(x) = f(x) - uc(x) + \\frac{1}{2\\mu} \\left(c(x)\\right)^2\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "where $u,\\mu >0$ are auxiliary parameters. The gradient for the augmented Lagrangian is\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "\\nabla \\Psi_{u,\\mu}(x) = \\nabla f(x) - u\\nabla c(x) + \\frac{1}{\\mu}c(x)\\nabla c(x)\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "and the Hessian is given by\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "\\nabla^2 \\Psi_{u,\\mu}(x) = \\nabla^2 f(x) -u\\nabla^2 c(x) + \\frac{1}{\\mu}c(x)\\nabla^2 c(x)  + \\frac{1}{\\mu}\\nabla c(x)\\nabla c(x)^T\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "#### Coding Task:\n",
    "\n",
    "For different fixed values of $u,\\mu$, apply your favourite SciPy optimization routine to minimize $\\Psi_{u,\\mu}(x)$ for the following objectives:\n",
    "\n",
    "1. $f(x) = \\frac{1}{2}(ax_1^2 + x_2^2) \\qquad\\qquad\\qquad$ with $\\; c(x) = x_1 + x_2^2 - 1$ \n",
    "\n",
    "2. $f(x) = (a - x_1)^2 + b(x_2 - x_1^2)^2 \\qquad$ with $\\; c(x) = x_1 + x_2^2 - 1$ \n",
    "\n",
    "starting at $(x_1,x_2) = (1,1)$ for suitable choices of parameters $a,b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "a = 1\n",
    "u = \n",
    "mu = \n",
    "\n",
    "# Initial guess\n",
    "x0 = np.array(...)\n",
    "\n",
    "# Objective function, gradient and Hessian\n",
    "obj = lambda x: \n",
    "gobj = lambda x: \n",
    "hobj = lambda x: \n",
    "\n",
    "# Constraint function, gradient and Hessian\n",
    "con = lambda x: \n",
    "gcon = lambda x: \n",
    "hcon = lambda x: \n",
    "\n",
    "# Penalty function, gradient and Hessian\n",
    "fun = lambda x: \n",
    "jac = lambda x: \n",
    "hess = lambda x: \n",
    "\n",
    "# Plot function contours\n",
    "plt.figure()\n",
    "X = np.linspace(...)\n",
    "Y = np.linspace(...)\n",
    "Z = np.meshgrid(X,Y)\n",
    "plt.contour(...)\n",
    "plt.colorbar()\n",
    "\n",
    "# Plot constraint\n",
    "plt.contour(X,Y,con(Z),[0],colors='r')\n",
    "\n",
    "# Call SciPy's Newton-CG\n",
    "xprev = x0 # for plotting\n",
    "res = opt.minimize(..., method=..., callback=callback)\n",
    "\n",
    "# Print results and show plot\n",
    "print(res)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "a = 1\n",
    "b = 100\n",
    "u = \n",
    "mu = \n",
    "\n",
    "# Initial guess\n",
    "x0 = np.array(...)\n",
    "\n",
    "# Objective function, gradient and Hessian\n",
    "obj = lambda x: \n",
    "gobj = lambda x: \n",
    "hobj = lambda x: \n",
    "\n",
    "# Constraint function, gradient and Hessian\n",
    "con = lambda x: \n",
    "gcon = lambda x: \n",
    "hcon = lambda x: \n",
    "\n",
    "# Penalty function, gradient and Hessian\n",
    "fun = lambda x: \n",
    "jac = lambda x: \n",
    "hess = lambda x: \n",
    "\n",
    "# Plot function contours\n",
    "plt.figure()\n",
    "X = np.linspace(...)\n",
    "Y = np.linspace(...)\n",
    "Z = np.meshgrid(X,Y)\n",
    "plt.contour(...)\n",
    "plt.colorbar()\n",
    "\n",
    "# Plot constraint\n",
    "plt.contour(X,Y,con(Z),[0],colors='r')\n",
    "\n",
    "# Call SciPy's Newton-CG\n",
    "xprev = x0 # for plotting\n",
    "res = opt.minimize(..., method=..., callback=callback)\n",
    "\n",
    "# Print results and show plot\n",
    "print(res)\n",
    "plt.show()"
   ]
  },
  {
   "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.8.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}