{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Nonlinear constrained optimization"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The nonlinear constrained optimization interface in\n",
    "`Optim` assumes that the user can write the optimization\n",
    "problem in the following way."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "$$\n",
    "\\min_{x\\in\\mathbb{R}^n} f(x) \\quad \\text{such that}\\\\\n",
    "l_x \\leq \\phantom{c(}x\\phantom{)} \\leq u_x \\\\\n",
    "l_c \\leq c(x) \\leq u_c.\n",
    "$$"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "For equality constraints on $x_j$ or $c(x)_j$ you set those\n",
    "particular entries of bounds to be equal, $l_j=u_j$.\n",
    "Likewise, setting $l_j=-\\infty$ or $u_j=\\infty$ means that the\n",
    "constraint is unbounded from below or above respectively."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Optim, NLSolversBase #hide\n",
    "import NLSolversBase: clear! #hide"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Constrained optimization with `IPNewton`"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We will go through examples on how to use the constraints interface\n",
    "with the interior-point Newton optimization algorithm [IPNewton](../../algo/ipnewton.md)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Throughout these examples we work with the standard Rosenbrock function.\n",
    "The objective and its derivatives are given by"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "fun(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2\n",
    "\n",
    "function fun_grad!(g, x)\n",
    "g[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]\n",
    "g[2] = 200.0 * (x[2] - x[1]^2)\n",
    "end\n",
    "\n",
    "function fun_hess!(h, x)\n",
    "h[1, 1] = 2.0 - 400.0 * x[2] + 1200.0 * x[1]^2\n",
    "h[1, 2] = -400.0 * x[1]\n",
    "h[2, 1] = -400.0 * x[1]\n",
    "h[2, 2] = 200.0\n",
    "end;"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Optimization interface\n",
    "To solve a constrained optimization problem we call the `optimize`\n",
    "method\n",
    "``` julia\n",
    "optimize(d::AbstractObjective, constraints::AbstractConstraints, initial_x::Tx, method::ConstrainedOptimizer, options::Options)\n",
    "```"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can create instances of `AbstractObjective` and\n",
    "`AbstractConstraints` using the types `TwiceDifferentiable` and\n",
    "`TwiceDifferentiableConstraints` from the package `NLSolversBase.jl`."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Box minimization\n",
    "We want to optimize the Rosenbrock function in the box\n",
    "$-0.5 \\leq x \\leq 0.5$, starting from the point $x_0=(0,0)$.\n",
    "Box constraints are defined using, for example,\n",
    "`TwiceDifferentiableConstraints(lx, ux)`."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": " * Status: success\n\n * Candidate solution\n    Final objective value:     2.500000e-01\n\n * Found with\n    Algorithm:     Interior Point Newton\n\n * Convergence measures\n    |x - x'|               = 4.39e-10 ≰ 0.0e+00\n    |x - x'|/|x'|          = 8.79e-10 ≰ 0.0e+00\n    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00\n    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00\n    |g(x)|                 = 1.00e+00 ≰ 1.0e-08\n\n * Work counters\n    Seconds run:   2  (vs limit Inf)\n    Iterations:    43\n    f(x) calls:    68\n    ∇f(x) calls:   68\n"
     },
     "metadata": {},
     "execution_count": 3
    }
   ],
   "cell_type": "code",
   "source": [
    "x0 = [0.0, 0.0]\n",
    "df = TwiceDifferentiable(fun, fun_grad!, fun_hess!, x0)\n",
    "\n",
    "lx = [-0.5, -0.5]; ux = [0.5, 0.5]\n",
    "dfc = TwiceDifferentiableConstraints(lx, ux)\n",
    "\n",
    "res = optimize(df, dfc, x0, IPNewton())"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "source": [
    "Like the rest of Optim, you can also use `autodiff=:forward` and just pass in\n",
    "`fun`."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "If we only want to set lower bounds, use `ux = fill(Inf, 2)`"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": " * Status: success (objective increased between iterations)\n\n * Candidate solution\n    Final objective value:     7.987239e-20\n\n * Found with\n    Algorithm:     Interior Point Newton\n\n * Convergence measures\n    |x - x'|               = 3.54e-10 ≰ 0.0e+00\n    |x - x'|/|x'|          = 3.54e-10 ≰ 0.0e+00\n    |f(x) - f(x')|         = 2.40e-19 ≰ 0.0e+00\n    |f(x) - f(x')|/|f(x')| = 3.00e+00 ≰ 0.0e+00\n    |g(x)|                 = 8.83e-09 ≤ 1.0e-08\n\n * Work counters\n    Seconds run:   0  (vs limit Inf)\n    Iterations:    35\n    f(x) calls:    63\n    ∇f(x) calls:   63\n"
     },
     "metadata": {},
     "execution_count": 4
    }
   ],
   "cell_type": "code",
   "source": [
    "ux = fill(Inf, 2)\n",
    "dfc = TwiceDifferentiableConstraints(lx, ux)\n",
    "\n",
    "clear!(df)\n",
    "res = optimize(df, dfc, x0, IPNewton())"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Defining \"unconstrained\" problems"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "An unconstrained problem can be defined either by passing\n",
    "`Inf` bounds or empty arrays.\n",
    "**Note that we must pass the correct type information to the empty `lx` and `ux`**"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": " * Status: success\n\n * Candidate solution\n    Final objective value:     5.998937e-19\n\n * Found with\n    Algorithm:     Interior Point Newton\n\n * Convergence measures\n    |x - x'|               = 1.50e-09 ≰ 0.0e+00\n    |x - x'|/|x'|          = 1.50e-09 ≰ 0.0e+00\n    |f(x) - f(x')|         = 1.80e-18 ≰ 0.0e+00\n    |f(x) - f(x')|/|f(x')| = 3.00e+00 ≰ 0.0e+00\n    |g(x)|                 = 7.92e-09 ≤ 1.0e-08\n\n * Work counters\n    Seconds run:   0  (vs limit Inf)\n    Iterations:    34\n    f(x) calls:    63\n    ∇f(x) calls:   63\n"
     },
     "metadata": {},
     "execution_count": 5
    }
   ],
   "cell_type": "code",
   "source": [
    "lx = fill(-Inf, 2); ux = fill(Inf, 2)\n",
    "dfc = TwiceDifferentiableConstraints(lx, ux)\n",
    "\n",
    "clear!(df)\n",
    "res = optimize(df, dfc, x0, IPNewton())\n",
    "\n",
    "lx = Float64[]; ux = Float64[]\n",
    "dfc = TwiceDifferentiableConstraints(lx, ux)\n",
    "\n",
    "clear!(df)\n",
    "res = optimize(df, dfc, x0, IPNewton())"
   ],
   "metadata": {},
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Generic nonlinear constraints"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We now consider the Rosenbrock problem with a constraint on\n",
    "$$\n",
    "   c(x)_1 = x_1^2 + x_2^2.\n",
    "$$"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We pass the information about the constraints to `optimize`\n",
    "by defining a vector function `c(x)` and its Jacobian `J(x)`."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The Hessian information is treated differently, by considering the\n",
    "Lagrangian of the corresponding slack-variable transformed\n",
    "optimization problem. This is similar to how the [CUTEst\n",
    "library](https://github.com/JuliaSmoothOptimizers/CUTEst.jl) works.\n",
    "Let $H_j(x)$ represent the Hessian of the $j$th component\n",
    "$c(x)_j$ of the generic constraints.\n",
    "and $\\lambda_j$ the corresponding dual variable in the\n",
    "Lagrangian. Then we want the `constraint` object to\n",
    "add the values of $H_j(x)$ to the Hessian of the objective,\n",
    "weighted by $\\lambda_j$."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The Julian form for the supplied function $c(x)$ and the derivative\n",
    "information is then added in the following way."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "con_c!(c, x) = (c[1] = x[1]^2 + x[2]^2; c)\n",
    "function con_jacobian!(J, x)\n",
    "    J[1,1] = 2*x[1]\n",
    "    J[1,2] = 2*x[2]\n",
    "    J\n",
    "end\n",
    "function con_h!(h, x, λ)\n",
    "    h[1,1] += λ[1]*2\n",
    "    h[2,2] += λ[1]*2\n",
    "end;"
   ],
   "metadata": {},
   "execution_count": 6
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Note that `con_h!` adds the `λ`-weighted Hessian value of each\n",
    "element of `c(x)` to the Hessian of `fun`.**"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can then optimize the Rosenbrock function inside the ball of radius\n",
    "$0.5$."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": " * Status: success\n\n * Candidate solution\n    Final objective value:     2.966216e-01\n\n * Found with\n    Algorithm:     Interior Point Newton\n\n * Convergence measures\n    |x - x'|               = 0.00e+00 ≤ 0.0e+00\n    |x - x'|/|x'|          = 0.00e+00 ≤ 0.0e+00\n    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00\n    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00\n    |g(x)|                 = 7.71e-01 ≰ 1.0e-08\n\n * Work counters\n    Seconds run:   0  (vs limit Inf)\n    Iterations:    28\n    f(x) calls:    109\n    ∇f(x) calls:   109\n"
     },
     "metadata": {},
     "execution_count": 7
    }
   ],
   "cell_type": "code",
   "source": [
    "lx = Float64[]; ux = Float64[]\n",
    "lc = [-Inf]; uc = [0.5^2]\n",
    "dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,\n",
    "                                     lx, ux, lc, uc)\n",
    "res = optimize(df, dfc, x0, IPNewton())"
   ],
   "metadata": {},
   "execution_count": 7
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can add a lower bound on the constraint, and thus\n",
    "optimize the objective on the annulus with\n",
    "inner and outer radii $0.1$ and $0.5$ respectively."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "┌ Warning: Initial guess is not an interior point\n",
      "└ @ Optim ~/work/Optim.jl/Optim.jl/src/multivariate/solvers/constrained/ipnewton/ipnewton.jl:116\n",
      "\n",
      "Stacktrace:\n",
      "  [1] initial_state(method::IPNewton{typeof(Optim.backtrack_constrained_grad), Symbol}, options::Optim.Options{Float64, Nothing}, d::TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, constraints::TwiceDifferentiableConstraints{typeof(Main.var\"##226\".con_c!), typeof(Main.var\"##226\".con_jacobian!), typeof(Main.var\"##226\".con_h!), Float64}, initial_x::Vector{Float64})\n",
      "    @ Optim ~/work/Optim.jl/Optim.jl/src/multivariate/solvers/constrained/ipnewton/ipnewton.jl:117\n",
      "  [2] optimize\n",
      "    @ ~/work/Optim.jl/Optim.jl/src/multivariate/solvers/constrained/ipnewton/interior.jl:229 [inlined]\n",
      "  [3] optimize(d::TwiceDifferentiable{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, constraints::TwiceDifferentiableConstraints{typeof(Main.var\"##226\".con_c!), typeof(Main.var\"##226\".con_jacobian!), typeof(Main.var\"##226\".con_h!), Float64}, initial_x::Vector{Float64}, method::IPNewton{typeof(Optim.backtrack_constrained_grad), Symbol})\n",
      "    @ Optim ~/work/Optim.jl/Optim.jl/src/multivariate/solvers/constrained/ipnewton/interior.jl:229\n",
      "  [4] top-level scope\n",
      "    @ ~/work/Optim.jl/Optim.jl/docs/src/examples/generated/ipnewton_basics.ipynb:4\n",
      "  [5] eval\n",
      "    @ ./boot.jl:385 [inlined]\n",
      "  [6] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)\n",
      "    @ Base ./loading.jl:2076\n",
      "  [7] #44\n",
      "    @ ~/.julia/packages/Literate/sLusI/src/Literate.jl:894 [inlined]\n",
      "  [8] (::IOCapture.var\"#4#7\"{Core.TypeofBottom, Literate.var\"#44#45\"{String, Bool, Module, String}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, Base.PipeEndpoint, Base.PipeEndpoint})()\n",
      "    @ IOCapture ~/.julia/packages/IOCapture/Rzdxd/src/IOCapture.jl:161\n",
      "  [9] with_logstate(f::Function, logstate::Any)\n",
      "    @ Base.CoreLogging ./logging.jl:515\n",
      " [10] with_logger\n",
      "    @ ./logging.jl:627 [inlined]\n",
      " [11] capture(f::Literate.var\"#44#45\"{String, Bool, Module, String}; rethrow::Type, color::Bool, passthrough::Bool, capture_buffer::IOBuffer)\n",
      "    @ IOCapture ~/.julia/packages/IOCapture/Rzdxd/src/IOCapture.jl:158\n",
      " [12] execute_block(sb::Module, block::String; inputfile::String, fake_source::String, softscope::Bool)\n",
      "    @ Literate ~/.julia/packages/Literate/sLusI/src/Literate.jl:892\n",
      " [13] execute_block\n",
      "    @ ~/.julia/packages/Literate/sLusI/src/Literate.jl:877 [inlined]\n",
      " [14] execute_notebook(nb::Dict{Any, Any}; inputfile::String, fake_source::String, softscope::Bool)\n",
      "    @ Literate ~/.julia/packages/Literate/sLusI/src/Literate.jl:795\n",
      " [15] (::Literate.var\"#38#40\"{Dict{String, Any}})()\n",
      "    @ Literate ~/.julia/packages/Literate/sLusI/src/Literate.jl:774\n",
      " [16] cd(f::Literate.var\"#38#40\"{Dict{String, Any}}, dir::String)\n",
      "    @ Base.Filesystem ./file.jl:112\n",
      " [17] jupyter_notebook(chunks::Vector{Literate.Chunk}, config::Dict{String, Any})\n",
      "    @ Literate ~/.julia/packages/Literate/sLusI/src/Literate.jl:773\n",
      " [18] notebook(inputfile::String, outputdir::String; config::Dict{Any, Any}, kwargs::@Kwargs{execute::Bool})\n",
      "    @ Literate ~/.julia/packages/Literate/sLusI/src/Literate.jl:710\n",
      " [19] top-level scope\n",
      "    @ ~/work/Optim.jl/Optim.jl/docs/generate.jl:18\n",
      " [20] include(fname::String)\n",
      "    @ Base.MainInclude ./client.jl:489\n",
      " [21] top-level scope\n",
      "    @ ~/work/Optim.jl/Optim.jl/docs/make.jl:10\n",
      " [22] include(mod::Module, _path::String)\n",
      "    @ Base ./Base.jl:495\n",
      " [23] exec_options(opts::Base.JLOptions)\n",
      "    @ Base ./client.jl:318\n",
      " [24] _start()\n",
      "    @ Base ./client.jl:552\n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": " * Status: success\n\n * Candidate solution\n    Final objective value:     2.966216e-01\n\n * Found with\n    Algorithm:     Interior Point Newton\n\n * Convergence measures\n    |x - x'|               = 0.00e+00 ≤ 0.0e+00\n    |x - x'|/|x'|          = 0.00e+00 ≤ 0.0e+00\n    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00\n    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00\n    |g(x)|                 = 7.71e-01 ≰ 1.0e-08\n\n * Work counters\n    Seconds run:   0  (vs limit Inf)\n    Iterations:    34\n    f(x) calls:    158\n    ∇f(x) calls:   158\n"
     },
     "metadata": {},
     "execution_count": 8
    }
   ],
   "cell_type": "code",
   "source": [
    "lc = [0.1^2]\n",
    "dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,\n",
    "                                     lx, ux, lc, uc)\n",
    "res = optimize(df, dfc, x0, IPNewton())"
   ],
   "metadata": {},
   "execution_count": 8
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Note that the algorithm warns that the Initial guess is not an\n",
    "interior point.** `IPNewton` can often handle this, however, if the\n",
    "initial guess is such that `c(x) = u_c`, then the algorithm currently\n",
    "fails. We may fix this in the future."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Multiple constraints\n",
    "The following example illustrates how to add an additional constraint.\n",
    "In particular, we add a constraint function\n",
    "$$\n",
    "   c(x)_2 = x_2\\sin(x_1)-x_1\n",
    "$$"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function con2_c!(c, x)\n",
    "    c[1] = x[1]^2 + x[2]^2     ## First constraint\n",
    "    c[2] = x[2]*sin(x[1])-x[1] ## Second constraint\n",
    "    c\n",
    "end\n",
    "function con2_jacobian!(J, x)\n",
    "    # First constraint\n",
    "    J[1,1] = 2*x[1]\n",
    "    J[1,2] = 2*x[2]\n",
    "    # Second constraint\n",
    "    J[2,1] = x[2]*cos(x[1])-1.0\n",
    "    J[2,2] = sin(x[1])\n",
    "    J\n",
    "end\n",
    "function con2_h!(h, x, λ)\n",
    "    # First constraint\n",
    "    h[1,1] += λ[1]*2\n",
    "    h[2,2] += λ[1]*2\n",
    "    # Second constraint\n",
    "    h[1,1] += λ[2]*x[2]*-sin(x[1])\n",
    "    h[1,2] += λ[2]*cos(x[1])\n",
    "    # Symmetrize h\n",
    "    h[2,1]  = h[1,2]\n",
    "    h\n",
    "end;"
   ],
   "metadata": {},
   "execution_count": 9
  },
  {
   "cell_type": "markdown",
   "source": [
    "We generate the constraint objects and call `IPNewton` with\n",
    "initial guess $x_0 = (0.25,0.25)$."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": " * Status: success\n\n * Candidate solution\n    Final objective value:     1.000000e+00\n\n * Found with\n    Algorithm:     Interior Point Newton\n\n * Convergence measures\n    |x - x'|               = 6.90e-10 ≰ 0.0e+00\n    |x - x'|/|x'|          = 3.55e+08 ≰ 0.0e+00\n    |f(x) - f(x')|         = 1.38e-09 ≰ 0.0e+00\n    |f(x) - f(x')|/|f(x')| = 1.38e-09 ≰ 0.0e+00\n    |g(x)|                 = 2.00e+00 ≰ 1.0e-08\n\n * Work counters\n    Seconds run:   0  (vs limit Inf)\n    Iterations:    29\n    f(x) calls:    215\n    ∇f(x) calls:   215\n"
     },
     "metadata": {},
     "execution_count": 10
    }
   ],
   "cell_type": "code",
   "source": [
    "x0 = [0.25, 0.25]\n",
    "lc = [-Inf, 0.0]; uc = [0.5^2, 0.0]\n",
    "dfc = TwiceDifferentiableConstraints(con2_c!, con2_jacobian!, con2_h!,\n",
    "                                     lx, ux, lc, uc)\n",
    "res = optimize(df, dfc, x0, IPNewton())"
   ],
   "metadata": {},
   "execution_count": 10
  },
  {
   "cell_type": "markdown",
   "source": [
    "---\n",
    "\n",
    "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
   ],
   "metadata": {}
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.2"
  },
  "kernelspec": {
   "name": "julia-1.10",
   "display_name": "Julia 1.10.2",
   "language": "julia"
  }
 },
 "nbformat": 4
}