Multi Scenario MORDM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multi-Scenario MORDM\n",
    "\n",
    "Multi-scenario MORMD is an extension of normal MORDM to better include robustness considerations within the search phase. It starts from the scenario discovery results resulting from MORDM. Next, from the experiments within this box, a set of scenarios is selected. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Import all modules and set up the model and the formulation\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import networkx as nx\n",
    "\n",
    "from ema_workbench import (Model, CategoricalParameter,\n",
    "                           ScalarOutcome, IntegerParameter, RealParameter,\n",
    "                           MultiprocessingEvaluator, Policy, Scenario)\n",
    "\n",
    "from dike_model_function import DikeNetwork  # @UnresolvedImport\n",
    "\n",
    "\n",
    "def sum_over(*args):\n",
    "    return sum(args)\n",
    "\n",
    "from ema_workbench import (Model, )\n",
    "\n",
    "from ema_workbench.em_framework.evaluators import perform_experiments\n",
    "from ema_workbench.em_framework.samplers import sample_uncertainties\n",
    "from ema_workbench.util import ema_logging\n",
    "import time\n",
    "\n",
    "from problem_formulation import get_model_for_problem_formulation\n",
    "\n",
    "\n",
    "ema_logging.log_to_stderr(level = \"NOTSET\")\n",
    "\n",
    "dike_model, planning_steps = get_model_for_problem_formulation(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Set uncertainties and levers paramaters to be used later\n",
    "import copy\n",
    "uncertainties = copy.deepcopy(dike_model.uncertainties)\n",
    "levers = copy.deepcopy(dike_model.levers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Import previously done experiments\n",
    "experiments = pd.read_csv('outcomes/Experiments_from_Exploration_1000scenarios_100Policies_pf5.csv')\n",
    "outcomes_100000 = pd.read_csv('outcomes/Outcomes_totals_from_Exploration_1000scenarios_100Policies_pf5.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "experiments_DF = pd.DataFrame(experiments)\n",
    "experiments_DF = experiments_DF.iloc[:, 1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "outcomes_DF = pd.DataFrame(outcomes_100000)\n",
    "outcomes_DF = outcomes_DF.iloc[:, 1:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "there are a number of ways in which you can make an informed selection. In this notebook the five worst performing scenarios are selected in terms of deaths in the A5 Dikering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Look into the problem formulation of comparing A5 with A4\n",
    "outcomes_A45 = outcomes_DF.loc[:,['A.4_Total Dike Investment Costs','A.4_Total Expected Number of Deaths', \n",
    "                           'A.5_Expected Total Damage', 'A.5_Total Dike Investment Costs','A.5_Total Expected Number of Deaths', 'RfR Total Costs']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "# take worst values for total deaths in A5\n",
    "row_high_death = np.array(outcomes_A45.sort_values('A.5_Total Expected Number of Deaths', ascending=False)[0:5].index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "# also all we need are the uncertainty columns --> for maxrows\n",
    "experiments_uncert = experiments_DF.loc[:, ['A.0_ID flood wave shape', 'A.1_Bmax', 'A.1_Brate', 'A.1_pfail',\n",
    "       'A.2_Bmax', 'A.2_Brate', 'A.2_pfail', 'A.3_Bmax', 'A.3_Brate',\n",
    "       'A.3_pfail', 'A.4_Bmax', 'A.4_Brate', 'A.4_pfail', 'A.5_Bmax',\n",
    "       'A.5_Brate', 'A.5_pfail']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#Find the uncertainties of the scenarios with the worst outcome for deaths in A5\n",
    "selected = experiments_uncert.loc[row_high_death, experiments_uncert.columns[:]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Search for each scenario\n",
    "\n",
    "For each of the five selected scenarios, many-objective optimization is used to find a pareto approximate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we have to do the same thing for each scenario, it is convenient to wrap this in a function we can call with the scenario under which we want to optimize. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from ema_workbench import Scenario\n",
    "\n",
    "scenarios = [Scenario(f\"{index}\", **row) for index, row in selected.iterrows()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running optimisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Load packages for running the MORDM analyses\n",
    "from ema_workbench import MultiprocessingEvaluator, ema_logging, SequentialEvaluator\n",
    "from ema_workbench.em_framework.evaluators import BaseEvaluator\n",
    "\n",
    "from ema_workbench.em_framework.optimization import (HyperVolume,\n",
    "                                                     EpsilonProgress)\n",
    "\n",
    "from __future__ import (unicode_literals, print_function, absolute_import,\n",
    "                        division)\n",
    "\n",
    "\n",
    "from ema_workbench import (Model, MultiprocessingEvaluator,\n",
    "                           ScalarOutcome, IntegerParameter, optimize, Scenario, Constraint)\n",
    "\n",
    "from ema_workbench.em_framework.optimization import EpsilonProgress\n",
    "from ema_workbench.util import ema_logging\n",
    "\n",
    "from problem_formulation import get_model_for_problem_formulation\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "if __name__ == '__main__':\n",
    "    ema_logging.log_to_stderr(ema_logging.INFO)\n",
    "\n",
    "    model, steps = get_model_for_problem_formulation(3)\n",
    "    \n",
    "    #Record the epsilon convergence\n",
    "    convergence_metrics = [EpsilonProgress()]\n",
    "    \n",
    "    #Set epsilon values\n",
    "    espilon = [1e5, 0.001,1e5, 0.001,1e5, 0.001,1e5, 0.001,1e5, 0.001, 1e5, 1e5]\n",
    "    \n",
    "    #Set number of function evaluations\n",
    "    nfe = 15000\n",
    "    \n",
    "    #Loop over the 5 selection scenarios\n",
    "    for i in range(0,5):\n",
    "        #Select scenario\n",
    "        ref_scenario = scenarios[i]\n",
    "        \n",
    "        #Run model with said scenario\n",
    "        with MultiprocessingEvaluator(model) as evaluator:\n",
    "            results, convergence = evaluator.optimize(nfe=nfe, searchover='levers',\n",
    "                                                      epsilons=espilon,\n",
    "                                                      convergence=convergence_metrics,\n",
    "                                                      reference=ref_scenario)\n",
    "        \n",
    "        #Save the results to own files\n",
    "        results.to_csv(f\"outcomes/results_optimisation_scenario34_{i}.csv\")\n",
    "        convergence.to_csv(f'outcomes/convergence_optimisation_scenario34{i}.csv')\n",
    "        \n",
    "        #Plot and save convergence figures\n",
    "        fig, (ax1, ax2) = plt.subplots(ncols=2, sharex=True)\n",
    "        fig, ax1 = plt.subplots(ncols=1)\n",
    "        ax1.plot(convergence.epsilon_progress)\n",
    "        ax1.set_xlabel('nr. of generations')\n",
    "        ax1.set_ylabel('$\\epsilon$ progress')\n",
    "        sns.despine()\n",
    "        plt.savefig(f'figures/convergence_scenario34_{i}.png')\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.13"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}