{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "title: Solvers and Solutions\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Originally Contributed by**: Arpit Bhatia"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The purpose of this part of the tutorial is to introduce you to solvers and how to use them with JuMP. We'll also learn\n",
    "what to do with a problem after the solver has finished optimizing it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is a Solver?\n",
    "A solver is a software package that incorporates algorithms for finding solutions to one or more classes of problem.\n",
    "For example, GLPK, which we used in the previous tutorials is a solver for linear programming (LP) and mixed integer\n",
    "programming (MIP) problems. It incorporates algorithms such as the simplex method, interior-point method etc. JuMP\n",
    "currently supports a number of open-source and commercial solvers which can be viewed\n",
    "[here](http://www.juliaopt.org/JuMP.jl/v0.19.1/installation/#Getting-Solvers-1)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is MathOptInterface?\n",
    "Each mathematical optimization solver API has its own concepts and data structures for representing optimization models\n",
    "and obtaining results. However, it is often desirable to represent an instance of an optimization problem at a higher\n",
    "level so that it is easy to try using different solvers. MathOptInterface (MOI) is an abstraction layer designed to provide\n",
    "an interface to mathematical optimization solvers so that users do not need to understand multiple solver-specific\n",
    "APIs. MOI can be used directly, or through a higher-level modeling interface like JuMP."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that JuMP reexports MathOptInterface and\n",
    "you can use the shortcut MOI to refer to MathOptInterface in your code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interacting with solvers\n",
    "JuMP models can be created in three different modes: `AUTOMATIC`, `MANUAL` and `DIRECT`.\n",
    "We'll use the following LP to illustrate them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "& \\max_{x,y} & x + 2y \\\\\n",
    "& \\;\\;\\text{s.t.} & x + y &\\leq 1 \\\\\n",
    "& & 0\\leq x, y &\\leq 1 \\\\\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using JuMP\n",
    "using GLPK"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `AUTOMATIC` Mode\n",
    "#### With Optimizer\n",
    "This is the easiest method to use a solver in JuMP. In order to do so, we simply set the solver inside the Model constructor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_auto = Model(GLPK.Optimizer)\n",
    "@variable(model_auto, 0 <= x <= 1)\n",
    "@variable(model_auto, 0 <= y <= 1)\n",
    "@constraint(model_auto, x + y <= 1)\n",
    "@objective(model_auto, Max, x + 2y)\n",
    "optimize!(model_auto)\n",
    "objective_value(model_auto)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### No Optimizer (at first)\n",
    "It is also possible to create a JuMP model with no optimizer attached. After the model object is initialized empty\n",
    "and all its variables, constraints and objective are set, then we can attach the solver at `optimize!` time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_auto_no = Model()\n",
    "@variable(model_auto_no, 0 <= x <= 1)\n",
    "@variable(model_auto_no, 0 <= y <= 1)\n",
    "@constraint(model_auto_no, x + y <= 1)\n",
    "@objective(model_auto_no, Max, x + 2y)\n",
    "set_optimizer(model_auto_no, GLPK.Optimizer)\n",
    "optimize!(model_auto_no)\n",
    "objective_value(model_auto_no)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we can also enforce the automatic mode by passing `caching_mode = MOIU.AUTOMATIC` in the Model function call.\n",
    "### `MANUAL` Mode\n",
    "This mode is similar to the `AUTOMATIC` mode, but there are less protections from the user getting errors from the solver\n",
    "API. On the other side, nothing happens silently, which might give the user more control. It requires attaching the solver\n",
    "before the solve step using the `MOIU.attach_optimizer()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_manual = Model(GLPK.Optimizer, caching_mode = MOIU.MANUAL)\n",
    "@variable(model_manual, 0 <= x <= 1)\n",
    "@variable(model_manual, 0 <= y <= 1)\n",
    "@constraint(model_manual, x + y <= 1)\n",
    "@objective(model_manual, Max, x + 2y)\n",
    "MOIU.attach_optimizer(model_manual)\n",
    "optimize!(model_manual)\n",
    "objective_value(model_manual)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `DIRECT` Mode\n",
    "Some solvers are able to handle the problem data directly. This is common for LP/MIP solver but not very common for\n",
    "open-source conic solvers. In this case we do not set a optimizer, we set a backend which is more generic and is able\n",
    "to hold data and not only solving a model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_direct = direct_model(GLPK.Optimizer())\n",
    "@variable(model_direct, 0 <= x <= 1)\n",
    "@variable(model_direct, 0 <= y <= 1)\n",
    "@constraint(model_direct, x + y <= 1)\n",
    "@objective(model_direct, Max, x + 2y)\n",
    "optimize!(model_direct)\n",
    "objective_value(model_direct)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solver Options\n",
    "Many of the solvers also allow options to be passed in. However, these options are solver-specific. To find out the various\n",
    "options available, please check out the individual solver packages. Some examples for the CBC solver are given below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Cbc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To turn off printing (i.e. silence the solver),"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = Model(optimizer_with_attributes(Cbc.Optimizer, \"logLevel\" => 0));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To increase the maximum number of iterations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = Model(optimizer_with_attributes(Cbc.Optimizer, \"max_iters\" => 10000));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To set the solution timeout limit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = Model(optimizer_with_attributes(Cbc.Optimizer, \"seconds\" => 5));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Querying Solutions\n",
    "So far we have seen all the elements and constructs related to writing a JuMP optimization model. In this section we\n",
    "reach the point of what to do with a solved problem. JuMP follows closely the concepts defined in MathOptInterface to\n",
    "answer user questions about a finished call to `optimize!(model)`. The three main steps in querying a solution are\n",
    "given below. We'll use the model we created in `AUTOMATIC` mode with an optimizer attached in this section.\n",
    "### Termination Statuses\n",
    "Termination statuses are meant to explain the reason why the optimizer stopped executing in the most recent call\n",
    "to `optimize!`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "OPTIMAL::TerminationStatusCode = 1"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "termination_status(model_auto)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can view the different termination status codes by referring to the docs or though checking the possible types using\n",
    "the below command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Enum MathOptInterface.TerminationStatusCode:\n",
       "OPTIMIZE_NOT_CALLED = 0\n",
       "OPTIMAL = 1\n",
       "INFEASIBLE = 2\n",
       "DUAL_INFEASIBLE = 3\n",
       "LOCALLY_SOLVED = 4\n",
       "LOCALLY_INFEASIBLE = 5\n",
       "INFEASIBLE_OR_UNBOUNDED = 6\n",
       "ALMOST_OPTIMAL = 7\n",
       "ALMOST_INFEASIBLE = 8\n",
       "ALMOST_DUAL_INFEASIBLE = 9\n",
       "ALMOST_LOCALLY_SOLVED = 10\n",
       "ITERATION_LIMIT = 11\n",
       "TIME_LIMIT = 12\n",
       "NODE_LIMIT = 13\n",
       "SOLUTION_LIMIT = 14\n",
       "MEMORY_LIMIT = 15\n",
       "OBJECTIVE_LIMIT = 16\n",
       "NORM_LIMIT = 17\n",
       "OTHER_LIMIT = 18\n",
       "SLOW_PROGRESS = 19\n",
       "NUMERICAL_ERROR = 20\n",
       "INVALID_MODEL = 21\n",
       "INVALID_OPTION = 22\n",
       "INTERRUPTED = 23\n",
       "OTHER_ERROR = 24"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(typeof(MOI.OPTIMAL))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution Statuses\n",
    "These statuses indicate what kind of result is available to be queried from the model. It's possible that no result is\n",
    "available to be queried. We shall discuss more on the dual status and solutions in the Duality tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "FEASIBLE_POINT::ResultStatusCode = 1"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "primal_status(model_auto)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "FEASIBLE_POINT::ResultStatusCode = 1"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dual_status(model_auto)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we saw before, the result (solution) status codes can be viewed directly from Julia."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Enum MathOptInterface.ResultStatusCode:\n",
       "NO_SOLUTION = 0\n",
       "FEASIBLE_POINT = 1\n",
       "NEARLY_FEASIBLE_POINT = 2\n",
       "INFEASIBLE_POINT = 3\n",
       "INFEASIBILITY_CERTIFICATE = 4\n",
       "NEARLY_INFEASIBILITY_CERTIFICATE = 5\n",
       "REDUCTION_CERTIFICATE = 6\n",
       "NEARLY_REDUCTION_CERTIFICATE = 7\n",
       "UNKNOWN_RESULT_STATUS = 8\n",
       "OTHER_RESULT_STATUS = 9"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(typeof(MOI.FEASIBLE_POINT))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Obtaining Solutions\n",
    "Provided the primal status is not `MOI.NO_SOLUTION`, we can inspect the solution values and optimal cost."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "value(x) = 0.0\n",
      "value(y) = 1.0\n",
      "objective_value(model_auto) = 2.0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "2.0"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@show value(x)\n",
    "@show value(y)\n",
    "@show objective_value(model_auto)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since it is possible that no solution is available to be queried from the model, calls to `value` may throw errors.\n",
    "Hence, it is recommended to check for the presence of solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "The model was not solved correctly.",
     "output_type": "error",
     "traceback": [
      "The model was not solved correctly.",
      "",
      "Stacktrace:",
      " [1] error(::String) at ./error.jl:33",
      " [2] top-level scope at In[16]:16",
      " [3] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1091",
      " [4] execute_code(::String, ::String) at /home/mbesancon/.julia/packages/IJulia/a1SNk/src/execute_request.jl:27",
      " [5] execute_request(::ZMQ.Socket, ::IJulia.Msg) at /home/mbesancon/.julia/packages/IJulia/a1SNk/src/execute_request.jl:86",
      " [6] #invokelatest#1 at ./essentials.jl:710 [inlined]",
      " [7] invokelatest at ./essentials.jl:709 [inlined]",
      " [8] eventloop(::ZMQ.Socket) at /home/mbesancon/.julia/packages/IJulia/a1SNk/src/eventloop.jl:8",
      " [9] (::IJulia.var\"#15#18\")() at ./task.jl:356"
     ]
    }
   ],
   "source": [
    "model_no_solution = Model(GLPK.Optimizer)\n",
    "@variable(model_no_solution, 0 <= x <= 1)\n",
    "@variable(model_no_solution, 0 <= y <= 1)\n",
    "@constraint(model_no_solution, x + y >= 3)\n",
    "@objective(model_no_solution, Max, x + 2y)\n",
    "\n",
    "optimize!(model_no_solution)\n",
    "\n",
    "if termination_status(model_no_solution) == MOI.OPTIMAL\n",
    "    optimal_solution = value(x)\n",
    "    optimal_objective = objective_value(model_no_solution)\n",
    "elseif termination_status(model_no_solution) == MOI.TIME_LIMIT && has_values(model_no_solution)\n",
    "    suboptimal_solution = value(x)\n",
    "    suboptimal_objective = objective_value(model_no_solution)\n",
    "else\n",
    "    error(\"The model was not solved correctly.\")\n",
    "end"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.5.1",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}