SAA for stochastic knapsack problem.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sample Average Approximation for the stochastic knapsack problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "df3 = pd.read_excel (r'knapsack.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "w = pd.DataFrame(df3.iloc[1:6, 2:12].values)\n",
    "r = df3.iloc[10:16, 2].values\n",
    "q = df3.iloc[7,2]\n",
    "c = df3.iloc[8,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyomo.environ import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create model\n",
    "m = ConcreteModel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SETS ######################\n",
    "\n",
    "m.i = Set(initialize = np.arange(w.shape[0]))\n",
    "#m.k.pprint()\n",
    "\n",
    "m.j = Set(initialize = np.arange(w.shape[1]))\n",
    "#m.n.pprint()\n",
    "\n",
    "## PARAMETERS ###############\n",
    "\n",
    "m.N = Param(initialize = w.shape[1])\n",
    "\n",
    "# Transform matrix to dictionary\n",
    "w.to_dict()\n",
    "w = w.stack().to_dict()\n",
    "r = dict(enumerate(r))\n",
    "\n",
    "m.w = Param(m.i, m.j, initialize = w)\n",
    "#m.w.pprint()\n",
    "\n",
    "m.r = Param(m.i, initialize = r)\n",
    "#m.r.pprint()\n",
    "\n",
    "m.q = Param(initialize = q)\n",
    "#m.q.pprint()\n",
    "\n",
    "m.c = Param(initialize = c)\n",
    "#m.c.pprint()\n",
    "\n",
    "## VARIABLES ################\n",
    "\n",
    "m.x = Var(m.i, within = Binary, initialize = 0)\n",
    "#m.x.pprint()\n",
    "\n",
    "m.z = Var(m.j, within =NonNegativeReals, bounds = (0, None), initialize = 0)\n",
    "#m.z.pprint()\n",
    "\n",
    "## CONSTRAINTS ##############\n",
    "\n",
    "def penalties(model, j):\n",
    "    return(m.z[j] >= (sum(m.w[i,j] * m.x[i] for i in m.i) - m.q))\n",
    "m.pen = Constraint(m.j, rule = penalties)\n",
    "#m.pen.pprint()\n",
    "\n",
    "## OBJECTIVE ################\n",
    "\n",
    "def objective(model):\n",
    "    return(sum(m.r[i] * m.x[i] for i in m.i) - m.c/m.N * sum(m.z[j] for j in m.j))\n",
    "m.obj = Objective(rule = objective, sense = maximize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using license file /Users/fietekrutein/gurobi.lic\n",
      "Academic license - for non-commercial use only\n",
      "Read LP format model from file /var/folders/3f/ty6j5m0x59l2j7rp1ktlmbx80000gn/T/tmpzz_9j9xq.pyomo.lp\n",
      "Reading time = 0.00 seconds\n",
      "x16: 11 rows, 16 columns, 61 nonzeros\n",
      "Changed value of parameter IntFeasTol to 1e-09\n",
      "   Prev: 1e-05  Min: 1e-09  Max: 0.1  Default: 1e-05\n",
      "Changed value of parameter MIPGap to 0.0\n",
      "   Prev: 0.0001  Min: 0.0  Max: inf  Default: 0.0001\n",
      "Gurobi Optimizer version 9.0.0 build v9.0.0rc2 (mac64)\n",
      "Optimize a model with 11 rows, 16 columns and 61 nonzeros\n",
      "Model fingerprint: 0x2bd6ab04\n",
      "Variable types: 11 continuous, 5 integer (5 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 7e+00]\n",
      "  Objective range  [2e-01, 9e-01]\n",
      "  Bounds range     [1e+00, 1e+00]\n",
      "  RHS range        [1e+00, 2e+01]\n",
      "Presolve removed 1 rows and 1 columns\n",
      "Presolve time: 0.00s\n",
      "Presolved: 10 rows, 15 columns, 60 nonzeros\n",
      "Variable types: 10 continuous, 5 integer (5 binary)\n",
      "\n",
      "Root relaxation: objective 2.266291e+00, 4 iterations, 0.00 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0    2.26629    0    3          -    2.26629      -     -    0s\n",
      "H    0     0                       1.8751400    2.26629  20.9%     -    0s\n",
      "H    0     0                       1.9323400    2.26629  17.3%     -    0s\n",
      "     0     0    2.14923    0    3    1.93234    2.14923  11.2%     -    0s\n",
      "     0     0    2.12025    0    4    1.93234    2.12025  9.72%     -    0s\n",
      "     0     0    2.08974    0    4    1.93234    2.08974  8.15%     -    0s\n",
      "     0     0    2.03291    0    3    1.93234    2.03291  5.20%     -    0s\n",
      "     0     0    2.00279    0    3    1.93234    2.00279  3.65%     -    0s\n",
      "     0     0    1.99839    0    4    1.93234    1.99839  3.42%     -    0s\n",
      "     0     0    1.94320    0    2    1.93234    1.94320  0.56%     -    0s\n",
      "     0     0     cutoff    0         1.93234    1.93234  0.00%     -    0s\n",
      "\n",
      "Cutting planes:\n",
      "  Gomory: 1\n",
      "  Cover: 3\n",
      "  MIR: 9\n",
      "\n",
      "Explored 1 nodes (25 simplex iterations) in 0.01 seconds\n",
      "Thread count was 4 (of 4 available processors)\n",
      "\n",
      "Solution count 2: 1.93234 1.87514 \n",
      "\n",
      "Optimal solution found (tolerance 0.00e+00)\n",
      "Best objective 1.932340000000e+00, best bound 1.932340000000e+00, gap 0.0000%\n"
     ]
    }
   ],
   "source": [
    "from pyomo.opt import SolverFactory\n",
    "import pyomo.environ\n",
    "\n",
    "opt = SolverFactory('gurobi')\n",
    "opt.options['IntFeasTol']= 10e-10\n",
    "opt.options['MIPGap'] = 0\n",
    "\n",
    "results = opt.solve(m, load_solutions=True, tee = True)"
   ]
  }
 ],
 "metadata": {
  "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}