NLopt+ForwardDiff single call to simulator (using caching and DiffResults)

Optim NLopt combined objective constraints - final example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "379ac43f",
   "metadata": {},
   "source": [
    "This is a summarized extract of the [Optim NLopt combined objective constraints](Optim%20NLopt%20combined%20objective%20constraints.ipynb) notebook, with only the final proposition as a selfcontained example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0da1c39d-66a5-41d0-b977-83102bd7e8a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "using NLopt, ForwardDiff"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7641270-a814-48bd-92da-93d824c6db0a",
   "metadata": {},
   "source": [
    "## Example (toy) problem: Rosen-Suzuki Problem (HS43)\n",
    "\n",
    "Taken from *SAS IML documentation* [Chapter 11. Nonlinear Optimization Examples](http://www.math.wpi.edu/saspdf/iml/chap11.pdf)  (or [online version](http://documentation.sas.com/doc/fr/imlug/15.2/imlug_nonlinearoptexpls_sect005.htm)) which is itself citing Problem 43 of (Hock & Schittkowski 1981)\n",
    "\n",
    "Reference: Hock, W., and Schittkowski, K. (1981). *Test Examples for Nonlinear Programming Codes*. Vol. 187 of Lecture Notes in Economics and Mathematical Systems. Berlin: Springer-Verlag. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ec1d2a1-00ab-4930-8837-09c9a48af983",
   "metadata": {},
   "source": [
    "Objective and constraints:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e129abda-fb36-4ccb-94d5-185082cfec26",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "g3 (generic function with 1 method)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f(x) = x[1]^2 + x[2]^2 + 2x[3]^2 + x[4]^2 - 5x[1] - 5x[2] - 21x[3] + 7x[4]\n",
    "g1(x) = -8 + x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2 + x[1] - x[2] + x[3] - x[4]\n",
    "g2(x) = -10 + x[1]^2 + 2x[2]^2 + x[3]^2 + 2x[4]^2 - x[1] - x[4]\n",
    "g3(x) = -5 + 2x[1]^2 + x[2]^2 + x[3]^2 + 2x[1] - x[2] - x[4]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "424dcfa8-e99c-4848-b323-965eb0d3d68f",
   "metadata": {},
   "source": [
    "Now let's pretend these are all computed in one master function `simulate(x)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "56108d4a-1658-4706-be20-fdedea7534a5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(f = -44.0, g1 = 0.0, g2 = -1.0, g3 = 0.0)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Global counter of simulator calls:\n",
    "simcount = 0\n",
    "\n",
    "function simulate(x)\n",
    "    global simcount += 1\n",
    "    return (f=f(x), g1=g1(x), g2=g2(x), g3=g3(x))\n",
    "end\n",
    "xopt = [0.,1,2,-1]\n",
    "simulate(xopt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "577f8c76-e667-49af-9351-727d1fa2e055",
   "metadata": {},
   "source": [
    "Check the simulation counter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e934789a-61ba-4509-8b64-c0747604b154",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simcount"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61ae4af2-8451-4490-b4f3-28a20b9caa64",
   "metadata": {},
   "source": [
    "## 4) Combined objective and constraints, using memoization (caching) + single ForwardDiff call to get value and jacobian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "81ab1280-cdc3-4847-b660-f7c75e9c0ae2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Vector{Float64}:\n",
       " -44.0\n",
       "   0.0\n",
       "  -1.0\n",
       "   0.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"simulate returning a vector instead of a named tuple\"\"\"\n",
    "function simulate_vec(x)\n",
    "    [v for v in simulate(x)]\n",
    "end\n",
    "simulate_vec(xopt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8fc5b743-3e79-49cd-af72-809cddff4ef7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "cached_objcons_single"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\" create myf, mygvec functions which share a common cache of the last simulator calls\n",
    "\n",
    "only calls the simulator once to get value and jacobian, using DiffResults\n",
    "\"\"\"\n",
    "function cached_objcons_single()\n",
    "    x_prev = zeros(4)*NaN\n",
    "    val_prev = zeros(4)*NaN\n",
    "    Jac_prev = zeros(4,4)*NaN\n",
    "    result = DiffResults.JacobianResult(zeros(4), zeros(4)) # (output, input)\n",
    "    \n",
    "    function runsim(x)\n",
    "        x_prev[:] = x\n",
    "        result = ForwardDiff.jacobian!(result, simulate_vec, x);\n",
    "        val_prev[:] = DiffResults.value(result)\n",
    "        Jac_prev[:] = DiffResults.jacobian(result)\n",
    "    end\n",
    "    \n",
    "    function myf(x::Vector, grad::Vector)\n",
    "        if x != x_prev\n",
    "            # update cache\n",
    "            runsim(x)\n",
    "        end\n",
    "        if length(grad) > 0\n",
    "            grad[:] = Jac_prev[1,:]\n",
    "        end\n",
    "        return val_prev[1]\n",
    "    end\n",
    "\n",
    "    function mygvec(result::Vector, x::Vector, grad::Matrix)\n",
    "        if x != x_prev\n",
    "            # update cache\n",
    "            runsim(x)\n",
    "        end\n",
    "        if length(grad) > 0\n",
    "            # Watch out the transpose due to convention mismatch between NLopt and ForwardDiff.\n",
    "            # Without it, we get unfortunately no error,\n",
    "            # but the convergence is, of course, poor\n",
    "            grad[:] = transpose(Jac_prev[2:4,:])\n",
    "        end\n",
    "        result[:] = val_prev[2:4]\n",
    "    end\n",
    "    \n",
    "    return myf, mygvec\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc6909f5-81f6-4f65-962f-f2f339c14c6a",
   "metadata": {},
   "source": [
    "Create the cached objective and constraints functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "84f75caa-27f7-4e76-b98b-41c52ad4ce51",
   "metadata": {},
   "outputs": [],
   "source": [
    "myf_scache, mygvec_scache = cached_objcons_single();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf277e80-3e69-4ed7-a3d6-1d92377fb467",
   "metadata": {},
   "source": [
    "Build the NLopt Opt structure and solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "bc377aa3-e5bb-4e06-ba46-e2738dd9b7d7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "f* = -44.00000000000344 after 20 iterations (returned XTOL_REACHED)\n",
      "xopt = [-7.800235413268972e-9, 1.000000002578578, 2.000000004835217, -0.9999999936822117]\n",
      "simulator calls: 16 (0.8 per iteration)\n"
     ]
    }
   ],
   "source": [
    "# Solver choice: SLSQP, MMA\n",
    "opt = Opt(:LD_SLSQP, 4)\n",
    "#opt = Opt(:LD_MMA, 4)\n",
    "\n",
    "opt.xtol_rel = 1e-6\n",
    "\n",
    "# Set objective and constraints\n",
    "opt.min_objective = myf_scache\n",
    "inequality_constraint!(opt, mygvec_scache, [0., 0., 0.])\n",
    "\n",
    "# Starting point\n",
    "x0 = [2, 1., 1., 1.] # ⚠ SLSQP can get stuck when starting at [1,1,1,1]\n",
    "\n",
    "# Solve and display\n",
    "simcount = 0 # reset simulato call counter\n",
    "(minf,minx,ret) = optimize(opt, x0)\n",
    "numevals = opt.numevals # the number of function evaluations\n",
    "println(\"f* = $minf after $numevals iterations (returned $ret)\")\n",
    "println(\"xopt = $minx\")\n",
    "println(\"simulator calls: $simcount ($(simcount/numevals) per iteration)\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.0",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.0"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}

Optim NLopt combined objective constraints.ipynb