02_OptimizationTesting.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    " %matplotlib inline\n",
    "import seaborn\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "from BorealWeights import BorealWeightedProblem\n",
    "from gradutil import *\n",
    "from pyomo.opt import SolverFactory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialize data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "revenue, carbon, deadwood, HA = init_boreal()\n",
    "opt = SolverFactory('glpk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "carbon_clean = carbon.dropna(axis=0, how='any')\n",
    "HA_clean = HA.dropna(axis=0, how='any')\n",
    "deadwood_clean = deadwood.dropna(axis=0, how='any')\n",
    "revenue_clean = revenue.dropna(axis=0, how='any')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Problem formulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Carbon without Nans"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets solve the problem with just one objective and using data without Nan-values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "min(carbon_clean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "bproblem = BorealWeightedProblem(carbon_clean.values)\n",
    "opt.solve(bproblem.model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_solution(problem,data):\n",
    "    res_dict = dict()\n",
    "    for i in problem.model.I:\n",
    "        for j in problem.model.J:\n",
    "            res_dict[j] = res_dict.get(j,0) + int(problem.model.x[i,j].value)\n",
    "    print('Handling, # of stands')        \n",
    "    for key in res_dict:\n",
    "        print(\"{:8} {}\".format(list(data)[key], res_dict[key]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_solution(bproblem, carbon_clean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So it looks like we really are able to solve this problem using the original data! (Without Nans)\n",
    "And also the distribution of variables makes sense (Set aside storages carbon the most, so it is the most common management regime)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Carbon where Nan:s replaced with BAUs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets replace Nan:s with the corresponding BAU values and lets try to solve the problem. All the single optimization tasks are maximizing, so in the final solution there should be no zeros anyway."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_revenue = nan_to_bau(revenue)\n",
    "n_carbon = nan_to_bau(carbon)\n",
    "n_deadwood = nan_to_bau(deadwood)\n",
    "n_ha = nan_to_bau(HA)\n",
    "\n",
    "revenue_norm = new_normalize(n_revenue.values)\n",
    "carbon_norm = new_normalize(n_carbon.values)\n",
    "deadwood_norm = new_normalize(n_deadwood.values)\n",
    "ha_norm = new_normalize(n_ha.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "carbon_zeros = carbon_norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "zero_bproblem = BorealWeightedProblem(carbon_zeros)\n",
    "opt.solve(zero_bproblem.model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_solution(zero_bproblem, carbon)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution are still quite the same than with the previous try. Good."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Verifying that solution looks reasonable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "lst = values_to_list(zero_bproblem.model, n_carbon.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "min(lst)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no zeros at all, so at least by that aspect the result looks rational."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(lst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.min(carbon.dropna(axis=0, how='any').values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.max(carbon.dropna(axis=0, how='any').values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving all single objective optimization tasks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solving carbon storage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "carbon_bproblem = zero_bproblem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "carbon_values = values_to_list(carbon_bproblem.model, n_carbon.values)\n",
    "print(\"Minimum: {}, maximum: {}\".format(min(carbon_values), max(carbon_values)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.max(carbon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_solution(carbon_bproblem, n_carbon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "this_bproblem = carbon_bproblem\n",
    "sum([this_bproblem.model.x[i,j].value for i in this_bproblem.model.I for j in this_bproblem.model.J])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(carbon_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solving HA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "HA_zeros = ha_norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "HA_bproblem = BorealWeightedProblem(HA_zeros)\n",
    "opt.solve(HA_bproblem.model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "HA_values = values_to_list(HA_bproblem.model, n_ha.values)\n",
    "print(\"Minimum: {}, maximum: {}\".format(min(HA_values), max(HA_values)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_solution(HA_bproblem, n_ha)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(HA_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These results make sense also. The values match to the ones attained in the previous notebook, which is great."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solving deadwood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "deadwood_zeros = deadwood_norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "deadwood_bproblem = BorealWeightedProblem(deadwood_zeros)\n",
    "opt.solve(deadwood_bproblem.model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "deadwood_values = values_to_list(deadwood_bproblem.model, n_deadwood.values)\n",
    "print(\"Minimum: {}, maximum: {}\".format(min(deadwood_values), max(deadwood_values)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_solution(deadwood_bproblem, n_deadwood)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(deadwood_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matching."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solving Timber revenue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "revenue_zeros = revenue_norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "revenue_bproblem = BorealWeightedProblem(revenue_zeros)\n",
    "opt.solve(revenue_bproblem.model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "revenue_values = values_to_list(revenue_bproblem.model, n_revenue.values)\n",
    "print(\"Minimum: {}, maximum: {}\".format(min(revenue_values), max(revenue_values)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_solution(revenue_bproblem, n_revenue)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First of all, it would make sense that the BAU scheme should produce the most profit and so there should be more BAU handlings. Apparently this is not the case. The NT\\*R schemes are prolonged versions of BAU, so probably they actually are able to produce bigger timbers and so more revenue. I think this still is correct result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There were some stands with value 0.0, so lets check them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "rv = np.array(revenue_values)\n",
    "revenue.iloc[rv == 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apparently on some stands there is no way to make any profit, so it is ok that those are only zeros."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(revenue_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final objective value is still the right one."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing optimizations to ones in the papers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The single objective results are documented in http://onlinelibrary.wiley.com/doi/10.1111/1365-2664.12790/full\n",
    "so it is meaningful to compare our results to that one."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Results in the paper:\n",
    "\"The maximum capacity of the landscape \n",
    "- (i) to provide harvest revenues (NPV) was 250 M\u20ac (average 5800 \u20ac ha\u22121),\n",
    "- (ii) to store carbon was 4459 \u00d7 10\u00b3 MgC (average 10\u00b3 MgC ha\u22121), \n",
    "- (iii) for deadwood index was 218 150 m\u00b3 (average 5\u00b71 m\u00b3 ha\u22121) and \n",
    "- (iv) for the combined habitat availability was 20 211 (no units) (average 0\u00b747 ha\u22121).\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For us the correspondig values are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"(i) Harvest revenues {:.0f} M\u20ac\".format(sum(revenue_values)/1000000))\n",
    "print(\"(ii) Carbon storage {:.0f} x 10\u00b3 MgC\".format(sum(carbon_values)/1e+3))\n",
    "print(\"(iii) Deadwood index {:.0f} m3\".format(sum(deadwood_values)))\n",
    "print(\"(iv) Combined Habitat {:.0f}\".format(sum(HA_values)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is still something weird with the data values given in the paper:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Total ha-1 calculated according to the values given in paper:')\n",
    "print('-'*62)\n",
    "print('Revenue/(average timber revenue/ha-1) = {}'.format(250000000/5800))\n",
    "print('Carbon/(average carbon storage /ha-1) = {}'.format(4459*1000/103))\n",
    "print('Deadwood/(average deadwood index/ha-1) = {}'.format(218150/5.1))\n",
    "print('Combined habitat/(average habitat/ha-1) = {}'.format(20211/0.47))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now all these values indicate that there should be ~43000 hectars in total. Paper still states that there were 68 700 hectars. I don't really know if that is a real problem regarding the optimization task, but it's still odd."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets draw the same pictures than in the paper, so we can see if there are some deviations there also."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "pylab.rcParams['figure.figsize'] = (15,12)\n",
    "\n",
    "ind = list(revenue)\n",
    "ind.append('\u00d6ptim.')\n",
    "\n",
    "fig, ax = plt.subplots(2,2)\n",
    "val = revenue.sum().values\n",
    "val = np.append(val, sum(revenue_values))\n",
    "ax[0,0].bar(ind, val)\n",
    "ax[0,0].set_title('Timber Harvest Revenues')\n",
    "\n",
    "val = carbon.sum().values\n",
    "val = np.append(val, sum(carbon_values))\n",
    "ax[0,1].bar(ind, val)\n",
    "ax[0,1].set_title('Carbon storage')\n",
    "\n",
    "val = deadwood.sum().values\n",
    "val = np.append(val, sum(deadwood_values))\n",
    "ax[1,0].bar(ind, val)\n",
    "ax[1,0].set_title('Deadwood')\n",
    "\n",
    "val = HA.sum().values\n",
    "val = np.append(val, sum(HA_values))\n",
    "ax[1,1].bar(ind, val)\n",
    "ax[1,1].set_title('Habitat availability')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These really doesn't look like the same as in the paper. Something must be wrong with the way I am summing values without optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "Image('../images/articleOptims.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pay-off Table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.dstack((n_revenue, n_carbon, n_deadwood, n_ha))\n",
    "problems = [revenue_bproblem, carbon_bproblem, deadwood_bproblem, HA_bproblem]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "payoff = [[np.sum(values_to_list(p.model, data[:, :, i]))\n",
    "           for i in range(np.shape(data)[-1])]\n",
    "          for p in problems]\n",
    "z_ideal = np.max(payoff, axis=0)\n",
    "z_nadir = np.min(payoff, axis=0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "for p in payoff:\n",
    "    for f in p:\n",
    "        print('{:11.1f}'.format(f), end=' ')\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_ideal, z_nadir,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These look correct, so everything is fine"
   ]
  }
 ],
 "metadata": {
  "git": {
   "suppress_outputs": true
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}