Merge pull request #5 from RXGottlieb/master

Updated examples to EAGO v0.7.1

Author: RXGottlieb

Manifest.toml

@@ -0,0 +1,19 @@
+[deps]
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+EAGO = "bb8be931-2a91-5aca-9f87-79e1cb69959a"
+Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b"
+IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
+
+[compat]
+CSV = "0.10.4"
+DataFrames = "1.3.4"
+EAGO = "0.7.1"
+Gurobi = "0.11.3"
+IntervalArithmetic = "0.20.6"
+JuMP = "1.1.1"
+MathOptInterface = "1.5.0"

Project.toml

@@ -7,18 +7,21 @@
     "# Customizing EAGO to Solve a Quasiconvex Problem \n",
     "\n",
     "[Matthew Wilhelm](https://psor.uconn.edu/person/matthew-wilhelm/)  \n",
+    "Department of Chemical and Biomolecular Engineering, University of Connecticut\n",
+    "\n",
+    "[Robert Gottlieb](https://psor.uconn.edu/person/robert-gottlieb/)  \n",
     "Department of Chemical and Biomolecular Engineering, University of Connecticut"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "In the following example, we illustrate how EAGO's basic branch and bound routine can be adapted for use in other algorithms.\n",
+    "In the following example, we illustrate how EAGO's basic branch-and-bound routine can be adapted for use in other algorithms.\n",
     "\n",
     "## An Algorithm for Solving Quasiconvex Problems\n",
     "\n",
-    "We'll adapt EAGO to implement the bisection based algorithm used to solve quasiconvex optimization problem presented in [1] to solve the below problem:\n",
+    "We'll adapt EAGO to implement the bisection-based algorithm used to solve the quasiconvex optimization problem presented in [1]:\n",
     "\n",
     "$\n",
     "\\begin{align}\n",
@@ -50,35 +53,35 @@
     "\\end{align}\n",
     "$\n",
     "\n",
-    "Let $\\phi_\\tau(\\mathbf y) = f(\\mathbf y) - \\tau$ such that $\\tau = (t^L + t^U)/2$. We solve for $\\mathbf y$ subject to constraints (24)-(27) were $\\phi_\\tau(\\mathbf y) \\leq 0$. If this is feasible, $t^* \\in [t^L, \\tau]$, else $t^* \\in [\\tau, t^U]$. The other interval containg $t^*$ is kept and the other is fathomed. This manner of bisection is repeated until an interval containing a feasible solution with a width of at most $\\epsilon$ is located [2]."
+    "Let $\\phi_\\tau(\\mathbf y) = f(\\mathbf y) - \\tau$ such that $\\tau = (t^L + t^U)/2$. We solve for $\\mathbf y$ subject to constraints (24)-(27) where $\\phi_\\tau(\\mathbf y) \\leq 0$. If this is feasible, $t^* \\in [t^L, \\tau]$, else $t^* \\in [\\tau, t^U]$. The interval containg $t^*$ is kept and the other is fathomed. This manner of bisection is repeated until an interval containing a feasible solution with a width of at most $\\epsilon$ is located [2]."
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "## Customizing EAGO's script\n",
-    "First, the preprocessing step, upper problem, and postprocessing routines are short circuited as only a single optimization problem need be solved at each iteration."
+    "First, the preprocessing step, upper problem, and postprocessing routines are short-circuited as only a single optimization problem needs to be solved at each iteration."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [],
    "source": [
     "using MathOptInterface, EAGO, JuMP\n",
-    "import EAGO: Optimizer\n",
+    "import EAGO: Optimizer, GlobalOptimizer\n",
     "\n",
     "struct QuasiConvex <: EAGO.ExtensionType end\n",
     "import EAGO: preprocess!, upper_problem!, postprocess!\n",
-    "function EAGO.preprocess!(t::QuasiConvex, x::Optimizer)\n",
+    "function EAGO.preprocess!(t::QuasiConvex, x::GlobalOptimizer)\n",
     "    x._preprocess_feasibility = true\n",
     "    end\n",
-    "function EAGO.upper_problem!(t::QuasiConvex, x::Optimizer)\n",
+    "function EAGO.upper_problem!(t::QuasiConvex, x::GlobalOptimizer)\n",
     "    x._upper_feasibility = true\n",
     "    end\n",
-    "function EAGO.postprocess!(t::QuasiConvex, x::Optimizer)\n",
+    "function EAGO.postprocess!(t::QuasiConvex, x::GlobalOptimizer)\n",
     "    x._postprocess_feasibility = true\n",
     "end"
    ]
@@ -87,23 +90,24 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Next we specify that only an absolute tolerance should be checked for convergence and termination."
+    "Next, we specify that only an absolute tolerance should be checked for convergence and termination."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
-    "import EAGO: convergence_check, termination_check, repeat_check\n",
-    "function EAGO.convergence_check(t::QuasiConvex, x::Optimizer)\n",
+    "import EAGO: convergence_check, termination_check\n",
+    "function EAGO.convergence_check(t::QuasiConvex, x::GlobalOptimizer)\n",
     "    gap = (x._upper_objective_value - x._lower_objective_value)\n",
     "    return (gap <= x._parameters.absolute_tolerance)\n",
     "end\n",
-    "function EAGO.termination_check(t::QuasiConvex, x::Optimizer)\n",
+    "function EAGO.termination_check(t::QuasiConvex, x::GlobalOptimizer)\n",
     "    flag = EAGO.convergence_check(t, x)\n",
     "    if flag\n",
+    "        x._end_state = EAGO.GS_OPTIMAL\n",
     "        x._termination_status_code = MathOptInterface.OPTIMAL\n",
     "        x._result_status_code = MathOptInterface.FEASIBLE_POINT\n",
     "    end\n",
@@ -115,47 +119,48 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Only these sixth variable representing $t$ need be branched on."
+    "We then indicate that only the sixth variable, representing $t$, should be branched on. Since we will apply our knowledge about which $t^*$ should be kept in the lower problem definition, we also short-circuit EAGO's repeat_check function here to tell EAGO  not to branch this node, but instead to repeatedly evaluate it."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [],
    "source": [
+    "import EAGO: repeat_check\n",
     "branch_variable = [i == 6 for i=1:6]\n",
-    "EAGO.repeat_check(t::QuasiConvex, x::Optimizer) = true"
+    "EAGO.repeat_check(t::QuasiConvex, x::GlobalOptimizer) = true"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "In the lower problem, we then specify that the problem is solved locally for a fixed $t$ value. The objective value is \n",
-    "updated and the problem is contracted in order to discard the region which is known not to contain the optimal value."
+    "In the lower problem, we then specify that the problem is to be solved locally for a fixed $t$ value. The objective value is then updated and the problem is contracted in order to discard the region which is known to not contain the optimal value."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [],
    "source": [
     "import EAGO: lower_problem!\n",
-    "function EAGO.lower_problem!(t::QuasiConvex, x::Optimizer)\n",
+    "function EAGO.lower_problem!(t::QuasiConvex, x::GlobalOptimizer)\n",
     "    y = x._current_node\n",
-    "    lower = y.lower_variable_bounds[6]\n",
-    "    upper = y.upper_variable_bounds[6]\n",
+    "    indx = x._sol_to_branch_map[6]\n",
+    "    lower = y.lower_variable_bounds[indx]\n",
+    "    upper = y.upper_variable_bounds[indx]\n",
     "    midy = (lower + upper)/2.0\n",
-    "    y.lower_variable_bounds[6] = midy\n",
-    "    y.upper_variable_bounds[6] = midy\n",
+    "    y.lower_variable_bounds[indx] = midy\n",
+    "    y.upper_variable_bounds[indx] = midy\n",
     "    EAGO.solve_local_nlp!(x)\n",
     "    feas = x._upper_feasibility\n",
-    "    y.lower_variable_bounds[6] = feas ?  lower : midy\n",
-    "    y.upper_variable_bounds[6] = feas ?  midy : upper\n",
-    "    x._lower_objective_value = y.lower_variable_bounds[6]\n",
-    "    x._upper_objective_value = y.upper_variable_bounds[6]\n",
+    "    y.lower_variable_bounds[indx] = feas ?  lower : midy\n",
+    "    y.upper_variable_bounds[indx] = feas ?  midy : upper\n",
+    "    x._lower_objective_value = y.lower_variable_bounds[indx]\n",
+    "    x._upper_objective_value = y.upper_variable_bounds[indx]\n",
     "    x._lower_feasibility = true\n",
     "    return\n",
     "end"
@@ -165,45 +170,77 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We now build the JuMP model representing this problem, solve it and then retrieve the solution."
+    "We now define the optimizer factory to extend the core EAGO optimizer for this special problem. The SubSolvers() constructor is used to set the extension type (t), as well as the relaxed optimizer (r) and upper-bounding optimizer (u), if necessary. In this case, we will use the default solvers and only set the extension type."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "-------------------------------------------------------------------------------------------------------\n",
-      "|  Iteration #  |   Nodes  | Lower Bound | Upper Bound  |    Gap    |   Ratio   |  Time   | Time Left |\n",
-      "-------------------------------------------------------------------------------------------------------\n",
-      "First Solution Found at Node 1\n",
-      "UBD = -1.71690346673131\n",
-      "Solution is :\n",
-      "    X[1] = 0.6521533635051513\n",
-      "    X[2] = 0.6687269925064182\n",
-      "    X[3] = 0.1827815731700123\n",
-      "    X[4] = 0.24121532455559402\n",
-      "    X[5] = 0.2994373267499198\n",
-      "    X[6] = -1.71690346673131\n"
-     ]
-    },
     {
      "data": {
       "text/plain": [
-       "-1.71690346673131"
+       "#3 (generic function with 1 method)"
       ]
      },
-     "execution_count": 10,
+     "execution_count": 5,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "opt = with_optimizer(Optimizer, absolute_tolerance = 1E-8, branch_variable = branch_variable, ext_type = QuasiConvex())\n",
+    "factory = () -> Optimizer(SubSolvers(; t = QuasiConvex()))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We now build the JuMP model representing this problem, solve it, and retrieve the solution."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\n",
+      "******************************************************************************\n",
+      "This program contains Ipopt, a library for large-scale nonlinear optimization.\n",
+      " Ipopt is released as open source code under the Eclipse Public License (EPL).\n",
+      "         For more information visit https://github.com/coin-or/Ipopt\n",
+      "******************************************************************************\n",
+      "\n",
+      "-----------------------------------------------------------------------------------------------------------------------------\n",
+      "|  Iteration #  |     Nodes    | Lower Bound  |  Upper Bound  |      Gap     |     Ratio    |     Time     |    Time Left   |\n",
+      "-----------------------------------------------------------------------------------------------------------------------------\n",
+      " \n",
+      "Empty Stack: Exhaustive Search Finished\n",
+      "First Solution Found at Node 1\n",
+      "LBD = -1.7169029079377651\n",
+      "UBD = -1.7169028893113136\n",
+      "Solution is :\n",
+      "    X[1] = 0.65234819549671\n",
+      "    X[2] = 0.6685518965185859\n",
+      "    X[3] = 0.18303164806114738\n",
+      "    X[4] = 0.24118652470572963\n",
+      "    X[5] = 0.29934139369195156\n",
+      "    X[6] = -1.7169028893113136\n",
+      " \n",
+      "Global solution at y*=[0.65234819549671, 0.6685518965185859, 0.18303164806114738, 0.24118652470572963, 0.29934139369195156] with a value of f*=-1.7169028893113136"
+     ]
+    }
+   ],
+   "source": [
+    "opt = optimizer_with_attributes(factory, \n",
+    "                                \"absolute_tolerance\" => 1E-8, \n",
+    "                                \"branch_variable\" => branch_variable,\n",
+    "                                \"iteration_limit\" => 1000)\n",
     "m = Model(opt)\n",
     "@variable(m, ((i<6) ? 0 : -5) <= y[i=1:6] <= ((i<6) ? 5 : 0))\n",
     "@constraint(m, sum(i*y[i] for i=1:5) - 5 == 0)\n",
@@ -215,26 +252,20 @@
     "@NLconstraint(m, -expr2/(1 + sum(y[i]^2 for i=1:5)) - y[6] <= 0)\n",
     "@objective(m, Min, y[6])\n",
     "\n",
-    "JuMP.optimize!(m)# retrieve solution info\n",
+    "JuMP.optimize!(m) # Retrieve solution info\n",
     "solution = JuMP.value.(y[1:5])\n",
-    "global_obj_value = JuMP.value.(y[6])"
+    "global_obj_value = JuMP.value.(y[6])\n",
+    "print(\"Global solution at y*=$solution with a value of f*=$global_obj_value\")"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "### Reference:\n",
-    "1. C. Jansson,Quasiconvex relaxations based on interval arithmetic, Linear Algebra andits Applications, 324 (2001), pp. 27–53.\n",
-    "2. S. Boyd and L. Vandenberghe,Convex optimization, Cambridge university press,2004."
+    "1. C. Jansson, Quasiconvex relaxations based on interval arithmetic, Linear Algebra and its Applications, 324 (2001), pp. 27–53.\n",
+    "2. S. Boyd and L. Vandenberghe, Convex optimization, Cambridge University Press, 2004."
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
@@ -243,15 +274,15 @@
    "lastKernelId": null
   },
   "kernelspec": {
-   "display_name": "Julia 1.3.1",
+   "display_name": "Julia 1.7.3",
    "language": "julia",
-   "name": "julia-1.3"
+   "name": "julia-1.7"
   },
   "language_info": {
    "file_extension": ".jl",
    "mimetype": "application/julia",
    "name": "julia",
-   "version": "1.3.1"
+   "version": "1.7.3"
   }
  },
  "nbformat": 4,