Skip to content

Instantly share code, notes, and snippets.

@pierre-haessig

pierre-haessig/Optim NLopt combined objective constraints - final example.ipynb

Created July 19, 2021 12:11
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save pierre-haessig/2963d5b3ec69d2451dff286d55ff8990 to your computer and use it in GitHub Desktop.
Save pierre-haessig/2963d5b3ec69d2451dff286d55ff8990 to your computer and use it in GitHub Desktop.
Download ZIP
NLopt+ForwardDiff single call to simulator (using caching and DiffResults)
Display the source blob
Display the rendered blob
Raw
Optim NLopt combined objective constraints - final example.ipynb
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
Optim NLopt combined objective constraints.ipynb
{
"cells": [
{
"cell_type": "markdown",
"id": "7e777e4c",
"metadata": {},
"source": [
"# How to solve a constrained optimization problem with NLopt+ForwardDiff when objective and constraints are returned by one single function\n",
"\n",
"**Abstract**. This notebooks is an iterative exploration on how to **reduce the number of calls of the \"simulator\"** (the function which jointly computes the objective and the constraints) when using optimization algorithms from NLopt and gradients computed by ForwardDiff.\n",
"\n",
"There are two challenges:\n",
"\n",
"1. NLopt API for describing an optimization problem expects two separate functions: one for the objective and one (at least) for constraints\n",
"2. ForwardDiff doesn't returns by defaut the value of the function being differentiated (although the value is computed internally)\n",
"\n",
"Proposed solution for each issue:\n",
"1. Use a caching mechanism to call the simulator only once and then share results in separate objective and constraint functions\n",
"2. Use [DiffResults](https://github.com/JuliaDiff/DiffResults.jl/), a companion library which is specifically made to retreive all available results from one ForwardDiff call (e.g. value and gradient).\n",
"\n",
"PH, July 2021"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "563db6f6",
"metadata": {},
"outputs": [],
"source": [
"using NLopt, ForwardDiff"
]
},
{
"cell_type": "markdown",
"id": "b8dda1fb",
"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. \n",
"\n",
"$$\\min f(x) = x_1^2 + x_2^2 + 2x_3^2 + x_4^2 - 5x_1 - 5x_2 - 21x_3 + 7x_4 $$\n",
"\n",
"such that $g_1(x), g_2(x), g_3(x) \\leq 0$, with:\n",
"\n",
"$$ g_1(x) = -8 + x_1^2 + x_2^2 + x_3^2 + x_4^2 + x_1 - x_2 + x_3 - x_4 $$\n",
"$$ g_2(x) = -10 + x_1^2 + 2x_2^2 + x_3^2 + 2x_4^2 - x_1 - x_4 $$\n",
"$$ g_3(x) = -5 + 2x_1^2 + x_2^2 + x_3^2 + 2x_1 - x_2 - x_4 $$\n",
"\n",
"Known solution:\n",
"\n",
"starting from point (1,1,1,1), $x^*=(0,1,2,-1)$ with $f(x^*) = -44$"
]
},
{
"cell_type": "markdown",
"id": "0bb7b9a3",
"metadata": {},
"source": [
"## 1) Standard NLopt implementation (separate objective and constraints)\n",
"\n",
"- separate constraints \n",
"- gradient with ForwarDiff"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "4e04951d",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-44.0"
]
},
"execution_count": 9,
"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",
"# Check value at known optimum\n",
"xopt = [0.,1,2,-1]\n",
"f(xopt)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "6aa7f808",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.0, -1.0, 0.0)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"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]\n",
"\n",
"g1(xopt), g2(xopt), g3(xopt) # active, *not active, active at optimum"
]
},
{
"cell_type": "markdown",
"id": "99a6daeb",
"metadata": {},
"source": [
"Now let's pretend these are all computed in one master function `simulate(x)`:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "2c5de28c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(f = -44.0, g1 = 0.0, g2 = -1.0, g3 = 0.0)"
]
},
"execution_count": 11,
"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",
"simulate(xopt)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "b22ade84",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(f = -19.0, g1 = -4.0, g2 = -6.0, g3 = -1.0)"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"simulate([1., 1., 1., 1.])"
]
},
{
"cell_type": "markdown",
"id": "3a2d6405",
"metadata": {},
"source": [
"Check the simulation counter:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "2091b7a3",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"simcount"
]
},
{
"cell_type": "markdown",
"id": "945de1ef",
"metadata": {},
"source": [
"NLopt objective and constraint functions wrappers:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "ca7241ac",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"myg3 (generic function with 1 method)"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function myf(x::Vector, grad::Vector)\n",
" if length(grad) > 0\n",
" ForwardDiff.gradient!(grad, x -> simulate(x).f, x)\n",
" end\n",
" return simulate(x).f\n",
"end\n",
"\n",
"function myg1(x::Vector, grad::Vector)\n",
" if length(grad) > 0\n",
" ForwardDiff.gradient!(grad, x -> simulate(x).g1, x)\n",
" end\n",
" return simulate(x).g1\n",
"end\n",
"\n",
"function myg2(x::Vector, grad::Vector)\n",
" if length(grad) > 0\n",
" ForwardDiff.gradient!(grad, x -> simulate(x).g2, x)\n",
" end\n",
" return simulate(x).g2\n",
"end\n",
"\n",
"function myg3(x::Vector, grad::Vector)\n",
" if length(grad) > 0\n",
" ForwardDiff.gradient!(grad, x -> simulate(x).g3, x)\n",
" end\n",
" return simulate(x).g3\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "35ca9389",
"metadata": {},
"source": [
"Build the NLopt Opt structure and solve"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "37c1af3a",
"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: 144 (7.2 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\n",
"inequality_constraint!(opt, myg1, 0.)\n",
"inequality_constraint!(opt, myg2, 0.)\n",
"inequality_constraint!(opt, myg3, 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)\")"
]
},
{
"cell_type": "markdown",
"id": "9da08afc",
"metadata": {},
"source": [
"Observations:\n",
"\n",
"- 7.2 simulator calls per iteration for SLSQ\n",
"- 8.0 simulator calls per iteration for MMA\n",
"\n",
"ideally, we would like this number to be:\n",
"- 1 if we can get ForwardDiff to also return the value in addition to the gradient\n",
"- 2 if we can't (i.e. we would need one call for the values of f,g1,g2,g3 and one extra for the jacobian)"
]
},
{
"cell_type": "markdown",
"id": "946230d0",
"metadata": {},
"source": [
"## 2) Standard NLopt implementation (separate objective and one vector constraint)\n",
"\n",
"- [vector-valued constraints](https://github.com/JuliaOpt/NLopt.jl#vector-valued-constraints) instead of separate constraints\n",
"- gradient (and Jacobian) with ForwarDiff"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "c27800b1",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"([0.0, -1.0, 0.0], [-4.0, -6.0, -1.0])"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"vectorized constraints\"\n",
"function gvec(x)\n",
" simres = simulate(x)\n",
" return [simres.g1, simres.g2, simres.g3]\n",
"end\n",
"gvec(xopt), gvec([1., 1., 1., 1.])"
]
},
{
"cell_type": "markdown",
"id": "0336e807",
"metadata": {},
"source": [
"Wrapper for vectorized constraint"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "4d02c3fc",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"mygvec (generic function with 1 method)"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function mygvec(result::Vector, x::Vector, grad::Matrix)\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",
" ForwardDiff.jacobian!(transpose(grad), gvec, x)\n",
" end\n",
" result[:] = gvec(x)\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "fe310ee1",
"metadata": {},
"source": [
"Build the NLopt Opt structure:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "92f1b9a4",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"f* = -5.862096054300197 after 195 iterations (returned XTOL_REACHED)\n",
"xopt = [1.270156211871642, 0.6180339887498949, 0.4999999999999999, 1.381966011250105]\n",
"simulator calls: 780 (4.0 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\n",
"inequality_constraint!(opt, mygvec, [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)\")"
]
},
{
"cell_type": "markdown",
"id": "7830851d",
"metadata": {},
"source": [
"Observations:\n",
"\n",
"- 3.1 simulator calls per iteration for SLSQ\n",
"- 4.0 simulator calls per iteration for MMA"
]
},
{
"cell_type": "markdown",
"id": "0fe30931",
"metadata": {},
"source": [
"Conclusion: the **number of simulator calls is divided by two**.\n",
"\n",
"Reasonable explanation:\n",
"- previously: calling simulator separately for f, g1, g2 , g3 (and their gradients)\n",
"- now: calling only for f and gvec (half as much)\n",
"\n",
"Also, it seems that ForwardDiff is able to get the gradient/Jacobian in just one call, that is using a chunk size of 4."
]
},
{
"cell_type": "markdown",
"id": "0385fdf3",
"metadata": {},
"source": [
"Attempt to check the fact that ForwardDiff chunk size is 4: not very clear!"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "8ad20876",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"ForwardDiff.GradientConfig{ForwardDiff.Tag{var\"#17#18\", Float64}, Float64, 4, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var\"#17#18\", Float64}, Float64, 4}}}((Partials(1.0, 0.0, 0.0, 0.0), Partials(0.0, 1.0, 0.0, 0.0), Partials(0.0, 0.0, 1.0, 0.0), Partials(0.0, 0.0, 0.0, 1.0)), ForwardDiff.Dual{ForwardDiff.Tag{var\"#17#18\", Float64}, Float64, 4}[Dual{ForwardDiff.Tag{var\"#17#18\", Float64}}(8.487983165e-314,2.33419537056e-313,2.75859452865e-313,1.69759663346e-313,3.8195924246e-313), Dual{ForwardDiff.Tag{var\"#17#18\", Float64}}(4.66839074116e-313,5.94158821604e-313,1.6975966342e-313,1.69759663317e-313,1.485397054e-313), Dual{ForwardDiff.Tag{var\"#17#18\", Float64}}(1.4853970552e-313,1.27319747493e-313,7.63918484777e-313,8.9123823242e-313,1.18831764321e-312), Dual{ForwardDiff.Tag{var\"#17#18\", Float64}}(1.230757558997e-312,1.294417432514e-312,6.3659873744e-314,1.506617011926e-312,5.0e-324)])"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cfg = ForwardDiff.GradientConfig(x -> simulate(x).f, xopt);\n",
"cfg"
]
},
{
"cell_type": "markdown",
"id": "b6b16956",
"metadata": {},
"source": [
"## 3) Combined objective and constraints, using memoization (caching)\n",
"\n",
"using a closure to share a memory of the last call of the simulator"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "e1e2ae08",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4-element Vector{Float64}:\n",
" -44.0\n",
" 0.0\n",
" -1.0\n",
" 0.0"
]
},
"execution_count": 20,
"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": 21,
"id": "324b5296",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Matrix{Float64}:\n",
" -5.0 -3.0 -13.0 5.0\n",
" 1.0 1.0 5.0 -3.0\n",
" -1.0 4.0 4.0 -5.0\n",
" 2.0 1.0 4.0 -1.0"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ForwardDiff.jacobian(simulate_vec, xopt) # n outputs × n inputs"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "b1598225",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3×4 Matrix{Float64}:\n",
" 1.0 1.0 5.0 -3.0\n",
" -1.0 4.0 4.0 -5.0\n",
" 2.0 1.0 4.0 -1.0"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ForwardDiff.jacobian(gvec, xopt) # n outputs × n inputs"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "9d63aa9a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"cached_objcons"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\" create myf, mygvec functions which share a common cache of the last simulator calls\n",
"\"\"\"\n",
"function cached_objcons()\n",
" x_prev = zeros(4)*NaN\n",
" val_prev = zeros(4)*NaN\n",
" Jac_prev = zeros(4,4)*NaN\n",
"\n",
" \n",
" function runsim(x)\n",
" x_prev[:] = x\n",
" val_prev[:] = simulate_vec(x)\n",
" Jac_prev[:] = ForwardDiff.jacobian(simulate_vec, x)\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": "56884ef7",
"metadata": {},
"source": [
"Create the cached objective and constraints functions:"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "1594cbe5",
"metadata": {},
"outputs": [],
"source": [
"myf_cache, mygvec_cache = cached_objcons();"
]
},
{
"cell_type": "markdown",
"id": "39dab698",
"metadata": {},
"source": [
"Test of the cached wrappers"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "bdb15c67",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-44.0"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"myf_cache(xopt, [])"
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "d049d21f",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4-element Vector{Float64}:\n",
" -5.0\n",
" -3.0\n",
" -13.0\n",
" 5.0"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"J = zeros(4)\n",
"myf_cache(xopt, J)\n",
"J"
]
},
{
"cell_type": "markdown",
"id": "c6397f96",
"metadata": {},
"source": [
"Empty matrix (for not requesting the Jacobian):"
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "01798caf",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0×0 Matrix{Any}"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"reshape([],0,0)"
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "b6efb618",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 0.0\n",
" -1.0\n",
" 0.0"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g = zeros(3)\n",
"mygvec_cache(g, xopt, reshape([],0,0))\n",
"g"
]
},
{
"cell_type": "code",
"execution_count": 31,
"id": "7fd920f5",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×3 Matrix{Float64}:\n",
" 1.0 -1.0 2.0\n",
" 1.0 4.0 1.0\n",
" 5.0 4.0 4.0\n",
" -3.0 -5.0 -1.0"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"J = zeros(4,3)\n",
"mygvec_cache(g, xopt, J)\n",
"J"
]
},
{
"cell_type": "markdown",
"id": "d9c755d5",
"metadata": {},
"source": [
"Build the NLopt Opt structure and solve"
]
},
{
"cell_type": "code",
"execution_count": 56,
"id": "1601d0ff",
"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: 32 (1.6 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_cache\n",
"inequality_constraint!(opt, mygvec_cache, [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)\")"
]
},
{
"cell_type": "markdown",
"id": "519a5314",
"metadata": {},
"source": [
"Observations:\n",
"\n",
"- 1.6 simulator calls per iteration for SLSQ (3.1 without caching)\n",
"- 1.7 simulator calls per iteration for MMA (4.0 without caching)\n",
"\n",
"This is equal or better than a factor two improvement. Perhaps the solver is sometimes repeatidly calling the functions at the same spot"
]
},
{
"cell_type": "markdown",
"id": "2b0ce2da",
"metadata": {},
"source": [
"## 4) Combined objective and constraints, using memoization (caching) + single ForwardDiff call to get value and jacobian\n",
"\n",
"Same as the above, but now the `runsim(x)` closure should only calls the simulator **once** to get both the value and the jacobian. With ForwardDiff, this can be done using [DiffResults.jl](https://github.com/JuliaDiff/DiffResults.jl) ([doc](https://juliadiff.org/DiffResults.jl/stable/))."
]
},
{
"cell_type": "code",
"execution_count": 45,
"id": "bb37988a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"cached_objcons_single"
]
},
"execution_count": 45,
"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": "04febac0",
"metadata": {},
"source": [
"Create the cached objective and constraints functions:"
]
},
{
"cell_type": "code",
"execution_count": 47,
"id": "0ffb8f57",
"metadata": {},
"outputs": [],
"source": [
"myf_scache, mygvec_scache = cached_objcons_single();"
]
},
{
"cell_type": "markdown",
"id": "bfdcc638",
"metadata": {},
"source": [
"Test of the cached wrappers"
]
},
{
"cell_type": "code",
"execution_count": 43,
"id": "ada37309",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-44.0"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"myf_scache(xopt, [])"
]
},
{
"cell_type": "code",
"execution_count": 44,
"id": "13c2f533",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4-element Vector{Float64}:\n",
" 16.0\n",
" 18.0\n",
" 8.0\n",
" 26.0"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"J = zeros(4)\n",
"myf_scache(xopt, J)\n",
"J"
]
},
{
"cell_type": "code",
"execution_count": 50,
"id": "b284e566",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 0.0\n",
" -1.0\n",
" 0.0"
]
},
"execution_count": 50,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g = zeros(3)\n",
"mygvec_scache(g, xopt, reshape([],0,0))\n",
"g"
]
},
{
"cell_type": "code",
"execution_count": 51,
"id": "bf07cbdb",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×3 Matrix{Float64}:\n",
" 1.0 -1.0 2.0\n",
" 1.0 4.0 1.0\n",
" 5.0 4.0 4.0\n",
" -3.0 -5.0 -1.0"
]
},
"execution_count": 51,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"J = zeros(4,3)\n",
"mygvec_scache(g, xopt, J)\n",
"J"
]
},
{
"cell_type": "markdown",
"id": "135386ad",
"metadata": {},
"source": [
"Build the NLopt Opt structure and solve"
]
},
{
"cell_type": "code",
"execution_count": 58,
"id": "1587ee8f",
"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)\")"
]
},
{
"cell_type": "markdown",
"id": "55fe7302",
"metadata": {},
"source": [
"Observations:\n",
"\n",
"- 0.8 simulator calls per iteration for SLSQ (3.1 without caching)\n",
"- 0.87 simulator calls per iteration for MMA (4.0 without caching)\n",
"\n",
"Conclusion: **we have now hit the bottom line of 1 simulator call per iteration**. For some reason (related to the way algorithm iteration works?), we are in fact slightly below (0.8–0.9)."
]
}
],
"metadata": {
"@webio": {
"lastCommId": null,
"lastKernelId": null
},
"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": false,
"sideBar": true,
"skip_h1_title": true,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": true
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment