two-stage-problem.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A two stage stochastic program\n",
    "\n",
    "Assume that you are warming your house with fuel. Currently fuel is at 750€/1000 liters. \n",
    "You can store up to 1000l for the winter, but you don't know yet how much you will use, and how much it will cost to buy fuel during winter.\n",
    "\n",
    "We assume that there are $4$ scenarios, with probabilities $p$, futur fuel cost $c$ and demand $d$ given afterwards.\n",
    "We want to minimize the expected costs of warming yourself during winter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1x4 Array{Int64,2}:\n",
       " 500  700  800  1000"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using JuMP, Clp\n",
    "\n",
    "S = 4 \n",
    "p = [0.1 0.4 0.4 0.1]\n",
    "c0 = 750\n",
    "c1 = [700 770 900 900]\n",
    "d = [500 700 800 1000]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 1\n",
    "\n",
    "Solve the corresponding 2 stage program. Give the optimal value, and optimal controls."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "588000.0"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m_2S = Model(solver = ClpSolver())\n",
    "@variable(m_2S, u1)\n",
    "@variable(m_2S, u2[1:4])\n",
    "@variable(m_2S, cout[1:4])\n",
    "@constraint(m_2S, u1 <=1000)\n",
    "for i = 1:4\n",
    "    @constraint(m_2S, u1+u2[i] >= d[i])\n",
    "end\n",
    "for i = 1:4\n",
    "    @constraint(m_2S, u2[i] >= 0)\n",
    "end\n",
    "@constraint(m_2S, u1 >= 0)\n",
    "for i=1:4\n",
    "    @constraint(m_2S, cout[i] >= u1*c0 + c1[i]*u2[i])\n",
    "end\n",
    "@objective(m_2S, Min, sum(p*cout))\n",
    "solve(m_2S)\n",
    "value_2S = getvalue(sum(p*cout))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Resolution gives us optimal value for u1 = 700 L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "700.0"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "getvalue(u1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 2\n",
    "\n",
    "Solve the open-loop version of this problem. Give the optimal value, and optimal controls.\n",
    "Check the inequality between the value of the open-loop problem and 2-stage problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "750000.0"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m_OL = Model(solver = ClpSolver())\n",
    "@variable(m_OL, u0)\n",
    "@variable(m_OL, u1)\n",
    "@variable(m_OL, cout[1:4])\n",
    "@constraint(m_OL, u0 <=1000)\n",
    "for i =1:4\n",
    "    @constraint(m_OL, u0+u1 >= d[i])\n",
    "end\n",
    "\n",
    "@constraint(m_OL, u0 >= 0)\n",
    "@constraint(m_OL, u1 >= 0)\n",
    "\n",
    "for i=1:4\n",
    "    @constraint(m_OL, cout[i] >= u0*c0 + c1[i]*u1)\n",
    "end\n",
    "@objective(m_OL, Min, sum(p*cout))\n",
    "solve(m_OL)\n",
    "value_OL = getvalue(sum(p*cout))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, as we don't allow u1 to be chosen according to the need in fuel, we decide at stage 0 how much fuel we need to get. As the fuel price is exepected to rise above c0, the solution is to buy as much fuel as we can at stage 0, i.e. u0 = 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1000.0"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "getvalue(u0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 3\n",
    "\n",
    "Solve the anticipative version of this problem. Give the optimal value, and optimal controls. Check the inequality between the value of the anticipative problem and 2-stage problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "560000.0"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m_a = Model(solver = ClpSolver())\n",
    "@variable(m_a, u0[1:4])\n",
    "@variable(m_a, u1[1:4])\n",
    "@variable(m_a, cout[1:4])\n",
    "for i =1:4\n",
    "    @constraint(m_a, u0[i] <=1000)\n",
    "end\n",
    "for i = 1:4\n",
    "    @constraint(m_a, u0[i]+u1[i] >= d[i])\n",
    "end\n",
    "for i = 1:4\n",
    "    @constraint(m_a, u0[i] >= 0)\n",
    "    @constraint(m_a, u1[i] >= 0)\n",
    "    \n",
    "end\n",
    "\n",
    "for i=1:4\n",
    "    @constraint(m_a, cout[i] >= u0[i]*c0 + c1[i]*u1[i])\n",
    "end\n",
    "@objective(m_a, Min, sum(p*cout))\n",
    "solve(m_a)\n",
    "value_a = getvalue(sum(p*cout))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1-element Array{Float64,1}:\n",
       " 560000.0"
      ]
     },
     "execution_count": 110,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "getvalue(p*cout)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected we each time buy exactly the quantity of fuel needed in each sceanario, at the best price. We have the inequality :\n",
    "\n",
    "$$value_a \\leq value_{2S} \\leq value_{OL} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 4\n",
    "\n",
    "We call (P_mean) the problem where cost and demand are replaced by their expectation. Solve this problem, giving value and first and second stage optimal control.\n",
    "\n",
    "Evaluate the value of this first stage control, that is the expected cost of using this first stage control with adapted recourse.\n",
    "\n",
    "Compare both values to the precedents problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "DimensionMismatch(\"matrix A has dimensions (1,4), matrix B has dimensions (1,4)\")",
     "output_type": "error",
     "traceback": [
      "DimensionMismatch(\"matrix A has dimensions (1,4), matrix B has dimensions (1,4)\")",
      "",
      " in _generic_matmatmul! at ./linalg/matmul.jl:454"
     ]
    }
   ],
   "source": [
    "m_mean = Model(solver = ClpSolver())\n",
    "d_mean = sum(p*d)\n",
    "c1_mean = mean(c1)\n",
    "@variable(m_mean, u0)\n",
    "@variable(m_mean, u1)\n",
    "@constraint(m_mean,u0 <=1000)\n",
    "@constraint(m_mean, u0 >= 0 )\n",
    "@constraint(m_mean, u1 >= 0)\n",
    "@constraint(m_mean, u0 + u1 >= d_mean)\n",
    "@objective(m_mean, Min, c0*u0 + c1_mean * u1)\n",
    "solve(m_mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "750.0"
      ]
     },
     "execution_count": 122,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "getvalue(u0)\n",
    "#Expected value u0 = 750"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now evaluating the expected cost of using such a method to estimate the quantity of fuel needed :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "603000.0"
     ]
    }
   ],
   "source": [
    "u1 = max(0,d - getvalue(u0))\n",
    "average_cost =0\n",
    "for i= 1:4\n",
    "    average_cost += p[i] * (u1[i]*c1[i] + getvalue(u0)*c0)\n",
    "end\n",
    "print(average_cost)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then have this inequality, reflecting the smart use of information it is possible to make\n",
    "\n",
    "$$value_a \\leq value_{2S} \\leq value_{average} \\leq value_{OL} $$"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 0.4.5",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}