notebooks/RCPSP tutorials/RCPSP-3 Local search.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1b2b937c",
   "metadata": {},
   "source": [
    "# Solving RCPSP with local search/metaheuristics/genetic algorithm\n",
    "\n",
    "One big family of combinatorial optimisation solver is local search. We include all metaheuristics and genetic algorithm into this simple terminology. \n",
    "\n",
    "In general, a local search algorithm explore the solution space by applying local changes to the current set of solutions.\n",
    "\n",
    "In the case of RCPSP, we have seen in [the first notebook](RCPSP%20%231%20Introduction.ipynb) that we can represent a solution with a priority list of tasks (equivalent to a permutation). \n",
    "Therefore local search algorithms are available on the search space being the ensemble of all permutation of tasks. We can imagine many kind of local changes possible to explore the permutation set.\n",
    "\n",
    "Local search (LS) algorithms are anytime algorithm, we have access to the current best solution whenever we want to stop the optimisation process. LS can't prove that a solution is optimal but it is rarely an issue in real world applications."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fb36d0f-a0c5-422c-b7d8-95e7435aac0c",
   "metadata": {},
   "source": [
    "## Prerequisites\n",
    "\n",
    "Concerning the python kernel to use for this notebook:\n",
    "- If running locally, be sure to use an environment with discrete-optimization and minizinc.\n",
    "- If running on colab, the next cell does it for you.\n",
    "- If running on binder, the environment should be ready.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cd42989-991e-4a40-8ca7-83ed9c84176c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# On Colab: install the library\n",
    "on_colab = \"google.colab\" in str(get_ipython())\n",
    "if on_colab:\n",
    "    import importlib\n",
    "    import os\n",
    "    import sys  # noqa: avoid having this import removed by pycln\n",
    "\n",
    "    !{sys.executable} -m pip install -U pip\n",
    "\n",
    "    # uninstall google protobuf conflicting with ray and sb3\n",
    "    ! pip uninstall -y protobuf\n",
    "\n",
    "    # install dev version for dev doc, or release version for release doc\n",
    "    !{sys.executable} -m pip install git+https://github.com/airbus/discrete-optimization@master#egg=discrete-optimization\n",
    "\n",
    "    # be sure to load the proper cffi (downgraded compared to the one initially on colab)\n",
    "    import cffi\n",
    "\n",
    "    importlib.reload(cffi)\n",
    "\n",
    "    # install and configure minizinc\n",
    "    !curl -o minizinc.AppImage -L https://github.com/MiniZinc/MiniZincIDE/releases/download/2.6.3/MiniZincIDE-2.6.3-x86_64.AppImage\n",
    "    !chmod +x minizinc.AppImage\n",
    "    !./minizinc.AppImage --appimage-extract\n",
    "    os.environ[\"PATH\"] = f\"{os.getcwd()}/squashfs-root/usr/bin/:{os.environ['PATH']}\"\n",
    "    os.environ[\"LD_LIBRARY_PATH\"] = (\n",
    "        f\"{os.getcwd()}/squashfs-root/usr/lib/:{os.environ['LD_LIBRARY_PATH']}\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19d8e0c6-ace8-45f8-ab8f-64b296d058be",
   "metadata": {},
   "source": [
    "### Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f451dd90-6e9e-4ae9-8749-40057c058a0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import nest_asyncio\n",
    "import numpy as np\n",
    "\n",
    "from discrete_optimization.datasets import fetch_data_from_psplib\n",
    "from discrete_optimization.generic_rcpsp_tools.ls_solver import (\n",
    "    LS_SOLVER,\n",
    "    LS_RCPSP_Solver,\n",
    ")\n",
    "from discrete_optimization.generic_tools.ea.ga import (\n",
    "    DeapCrossover,\n",
    "    DeapMutation,\n",
    "    Ga,\n",
    "    ObjectiveHandling,\n",
    ")\n",
    "from discrete_optimization.generic_tools.ls.hill_climber import (\n",
    "    HillClimber,\n",
    "    ModeMutation,\n",
    "    RestartHandler,\n",
    ")\n",
    "from discrete_optimization.generic_tools.ls.local_search import RestartHandlerLimit\n",
    "from discrete_optimization.generic_tools.ls.simulated_annealing import (\n",
    "    SimulatedAnnealing,\n",
    "    TemperatureSchedulingFactor,\n",
    ")\n",
    "from discrete_optimization.generic_tools.mutations.mixed_mutation import (\n",
    "    BasicPortfolioMutation,\n",
    ")\n",
    "from discrete_optimization.generic_tools.mutations.mutation_catalog import (\n",
    "    PermutationMutationRCPSP,\n",
    "    get_available_mutations,\n",
    ")\n",
    "from discrete_optimization.rcpsp.rcpsp_parser import get_data_available, parse_file\n",
    "from discrete_optimization.rcpsp.solver.rcpsp_pile import GreedyChoice, PileSolverRCPSP\n",
    "\n",
    "# patch asyncio so that applications using async functions can run in jupyter\n",
    "nest_asyncio.apply()\n",
    "\n",
    "# set logging level\n",
    "logging.basicConfig(level=logging.INFO)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c974812-e12c-4730-99d1-9b49d463badb",
   "metadata": {},
   "source": [
    "### Download datasets\n",
    "\n",
    "If not yet available, we import the datasets from [psplib](https://www.om-db.wi.tum.de/psplib/data.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a5255ec-2eb0-4674-9997-74de089be2dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "needed_datasets = [\"j301_1.sm\"]\n",
    "download_needed = False\n",
    "try:\n",
    "    files_available_paths = get_data_available()\n",
    "    for dataset in needed_datasets:\n",
    "        if len([f for f in files_available_paths if dataset in f]) == 0:\n",
    "            download_needed = True\n",
    "            break\n",
    "except:\n",
    "    download_needed = True\n",
    "\n",
    "if download_needed:\n",
    "    fetch_data_from_psplib()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6daae503",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parse some rcpsp file\n",
    "filepath = [f for f in get_data_available() if \"j1201_1.sm\" in f][0]\n",
    "rcpsp_problem = parse_file(filepath)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4787b053",
   "metadata": {},
   "source": [
    "## Hill climbing\n",
    "The hill climbing is the most basic local search algorithm. It explores the current best solution with a local move we can call `mutation`. If the new generated solution has a better objective value we overwrite the current best solution and repeat the process. The pseudocode is the following.\n",
    " \n",
    "1) Starts from a first solution $s$ and $obj=evaluation(s)$ the objective value of solution $s$ (to minimize)\n",
    "\n",
    "2) Compute $s'=mutate(s)$ and $obj'=evaluation(s')$\n",
    "\n",
    "3) if $obj'<obj$ then $(s,obj) \\leftarrow (s',obj')$\n",
    "\n",
    "4) Go back to 2) until some number of iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14c1b6ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "HillClimber?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "829ce2b1",
   "metadata": {},
   "source": [
    "The API of the HillClimber object is the following.\n",
    "```\n",
    "HillClimber(\n",
    "    problem: discrete_optimization.generic_tools.do_problem.Problem,\n",
    "    mutator: discrete_optimization.generic_tools.do_mutation.Mutation,\n",
    "    restart_handler: discrete_optimization.generic_tools.ls.local_search.RestartHandler,\n",
    "    mode_mutation: discrete_optimization.generic_tools.ls.local_search.ModeMutation,\n",
    "    params_objective_function: Union[discrete_optimization.generic_tools.do_problem.ParamsObjectiveFunction, NoneType] = None,\n",
    "    store_solution=False,\n",
    "    nb_solutions=1000)\n",
    "```\n",
    "We need to define a ```Mutation``` used in the algorithms.\n",
    "```RestartHandler``` is not used in HillClimbing, but will be described later in this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e99e46b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's compute the available mutation for our RCPSP Problem.\n",
    "_, mutations = get_available_mutations(rcpsp_problem)\n",
    "print(len(mutations), \" available for the problem \")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b01b2d2",
   "metadata": {},
   "source": [
    "We don't want to choose the mutation that we'll me using, we will build a \"meta\" mutation object that will in practice pick one of the mutation at each iteration (randomly). We call this meta mutation `BasicPortFolio` and can be initialized like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f263ab6",
   "metadata": {},
   "outputs": [],
   "source": [
    "dummy = rcpsp_problem.get_dummy_solution()\n",
    "list_mutation = [\n",
    "    mutate[0].build(rcpsp_problem, dummy, **mutate[1])\n",
    "    for mutate in mutations\n",
    "    if mutate[0] == PermutationMutationRCPSP\n",
    "]\n",
    "mixed_mutation = BasicPortfolioMutation(\n",
    "    list_mutation, weight_mutation=np.ones((len(list_mutation)))\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b58cc989",
   "metadata": {},
   "source": [
    "We can now initialize the HillClimbing solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca05e9b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "solver = HillClimber(\n",
    "    problem=rcpsp_problem,\n",
    "    mode_mutation=ModeMutation.MUTATE_AND_EVALUATE,\n",
    "    mutator=mixed_mutation,\n",
    "    restart_handler=RestartHandler(),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "388051af",
   "metadata": {},
   "source": [
    "Let's run the solver : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4201c17",
   "metadata": {},
   "outputs": [],
   "source": [
    "solver.solve??"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "093d25e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_hc = solver.solve(\n",
    "    initial_variable=dummy, nb_iteration_max=10000, max_time_seconds=20\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a63f20b",
   "metadata": {},
   "source": [
    "Let's plot the evolution of the solution found by the algorithm : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4863262",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1)\n",
    "ax.plot([x[1] for x in results_hc.list_solution_fits], marker=\"o\")\n",
    "ax.set_ylabel(\"- makespan\")\n",
    "ax.set_xlabel(\"# solution found\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f40b1af1",
   "metadata": {},
   "source": [
    "We improved a little bit the dummy solution :-)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6fc620a",
   "metadata": {},
   "source": [
    "## Simulated annealing\n",
    "[Simulated annealing](https://en.wikipedia.org/wiki/Simulated_annealing) is an algorithm from the class of metaheuristics. It is a famous method to escape local optima of a function because contrary to the Hill Climbing algorithm, the current solution can be different from the current best solution. \n",
    "The algorithm is very similar to the HC one, except for the 3) step.\n",
    "\n",
    "3) if $obj'<obj$ or $e^{\\frac{obj-obj'}{T}}\\geq random(0,1)$ then $(s, obj)\\leftarrow (s',obj')$ "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9eff86e",
   "metadata": {},
   "source": [
    "when $obj'$ is greater than $obj$, the exponent of the exponential is negative. \n",
    "When $T$, the temperature setting is high, it is very likely that the $e^{\\frac{obj-obj'}{T}}\\geq random(0,1)$ inequality is satisfied and that we update our current solution $s$ even though $obj'>obj$.\n",
    "To the contrary, with a low value $T$, we let less chance for $s'$ to be taken as new current solution. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45655ac9",
   "metadata": {},
   "source": [
    "- Simulated annealing relies heavily on initial temperature setting and the cooling schedule, which makes the temperature $T$ evolve through the algorithms : in the beginning we can allow high temperature to allow more exploration of the solution space, and then decrease the temperature to focus more on optimisation. in the library it will be done with `TemperatureScheduling` object, that can be custom.\n",
    "- For simulated annealing, the concept of `RestartHandler` makes sense : it consist in forcing the algorithm to go back to current best solution according to some criteria. Typically if there is no improvement for a given number of iteration, it can make sense to go back to current best solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73357bf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Come back to best solution every 300 iteration without improvement.\n",
    "restart_handler = RestartHandlerLimit(\n",
    "    nb_iteration_no_improvement=300,\n",
    ")\n",
    "# Multiply current temperature by a coefficient at every iteration. There might be better cooling schedule\n",
    "# in the litterature.\n",
    "temperature_scheduling = TemperatureSchedulingFactor(\n",
    "    temperature=10, restart_handler=restart_handler, coefficient=0.9999\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68047a88",
   "metadata": {},
   "source": [
    "We can now define the simulated annealing solver : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e40f760",
   "metadata": {},
   "outputs": [],
   "source": [
    "solver = SimulatedAnnealing(\n",
    "    problem=rcpsp_problem,\n",
    "    mutator=mixed_mutation,\n",
    "    temperature_handler=temperature_scheduling,\n",
    "    restart_handler=restart_handler,\n",
    "    mode_mutation=ModeMutation.MUTATE_AND_EVALUATE,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea316197",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_sa = solver.solve(\n",
    "    initial_variable=dummy, nb_iteration_max=10000, max_time_seconds=20\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a69fc0ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1)\n",
    "ax.plot([x[1] for x in results_sa.list_solution_fits], marker=\"o\")\n",
    "ax.set_ylabel(\"- makespan\")\n",
    "ax.set_xlabel(\"# solution found\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34bce2d8",
   "metadata": {},
   "source": [
    "Depending on your luck, you might have found a better solution than using the hillclimbing ;)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dacf1cfd",
   "metadata": {},
   "source": [
    "Since RCPSP is a widely developed use case in `discrete-optimisation` we developed a wrapper around Hill climbing and simulated annealing to ease their use : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea86fa0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "solver = LS_RCPSP_Solver(ls_solver=LS_SOLVER.HC, problem=rcpsp_problem)\n",
    "res = solver.solve(nb_iteration_max=10000, max_time_seconds=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f51d7dcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res.get_best_solution_fit())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "779b8e21",
   "metadata": {},
   "source": [
    "### Starting from a different initial solution\n",
    "We can use the results of a greedy solver as the initial solution for local search algorithm : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf317056",
   "metadata": {},
   "outputs": [],
   "source": [
    "solver = PileSolverRCPSP(problem=rcpsp_problem)\n",
    "res_greedy = solver.solve(greedy_choice=GreedyChoice.MOST_SUCCESSORS)\n",
    "print(\"res greedy : \", res_greedy.get_best_solution_fit())\n",
    "solver = LS_RCPSP_Solver(ls_solver=LS_SOLVER.SA, problem=rcpsp_problem)\n",
    "res = solver.solve(\n",
    "    nb_iteration_max=30000,\n",
    "    max_time_seconds=20,\n",
    "    starting_point=res_greedy.get_best_solution(),\n",
    ")\n",
    "print(\"local search : \", res.get_best_solution_fit())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b882a2d",
   "metadata": {},
   "source": [
    "## Genetic algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2278abab",
   "metadata": {},
   "source": [
    "`discrete-optimization` provides wrapper solver to ```deap``` library solver implement evolutionnary algorithms. Genetic algorithms can be seen as a generalisation of previous algorithms, the difference being that the algorithms are now handling a population of solutions instead of only the current one. \n",
    "To generate the new population for the current one the genetic algorithm is decomposed in several step : \n",
    "- selection : which individuals should be selected for the next generation\n",
    "- crossover : how to build new individuals from their parent (i.e the individual selected in the previous step)\n",
    "- mutate : how the new individual from previous step can acquire some mutation (i.e local changes)\n",
    "\n",
    "The genetic algorithm process is highly inspired by genetic selection and has been proven efficient for a wide range of optimisation problems."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32032e80",
   "metadata": {},
   "source": [
    "Currently the API of `Ga` is a bit different than others solvers, we need to specify the attribute of our object solution that we consider as a vector representation of the solution : in our case \"rcpsp_permutation\", we also have to specify manually what is the objective function we want to maximize."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5456137",
   "metadata": {},
   "outputs": [],
   "source": [
    "mutation = DeapMutation.MUT_SHUFFLE_INDEXES\n",
    "crossover = DeapCrossover.CX_UNIFORM_PARTIALY_MATCHED\n",
    "ga_solver = Ga(\n",
    "    rcpsp_problem,\n",
    "    encoding=\"rcpsp_permutation\",\n",
    "    objective_handling=ObjectiveHandling.AGGREGATE,\n",
    "    objectives=[\"makespan\"],\n",
    "    objective_weights=[-1],\n",
    "    pop_size=50,\n",
    "    max_evals=30000,\n",
    "    mut_rate=0.1,\n",
    "    crossover_rate=0.9,\n",
    "    crossover=crossover,\n",
    "    mutation=mutation,\n",
    ")\n",
    "results_ga = ga_solver.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b74fcb9",
   "metadata": {},
   "source": [
    "## Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ad8d1de",
   "metadata": {},
   "source": [
    "In this notebook we have shown 3 ways of optimizing the RCPSP problem using local search and evolutionnary algorithms. Most of those solvers are available for all implemented problem in `discrete-optimisation` as long as there is a vector representation of the solution that all those algorithms are able to optimize. The mentionned methods should not be discarded when you have a combinatorial optimisation problem, they're anytime solver and with the right parameters and mutations can become the most suited solution for your problem, especially when the problem is intractable with mathematical methods that we will discuss in next notebook."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.18"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "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
}