Lands_Validation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Validation of the optimal solution for LandS problem\n",
    "\n",
    "Disclaimer: I tried to keep it simple, using standard data structures and straigthforward formulations. Many steps can be done more efficiently. \n",
    "\n",
    "To validate the solution, we require and **exact** solver, that is, a solver that use rational numbers, so no rounding error is introduced at any step. In this case, we use QSopt."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import qsoptex\n",
    "import logging\n",
    "from fractions import Fraction\n",
    "\n",
    "#logging.basicConfig(level=logging.DEBUG)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We construct the data of the problem, using fractions instead of float"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "nscen = 1000000\n",
    "scen_ex = []\n",
    "prob_ex = []\n",
    "for i in range (1,101):\n",
    "    for j in range(1,101):\n",
    "        for k in range(1,101):\n",
    "            scen_ex.append([Fraction(4*(i-1),100),Fraction(4*(j-1),100),Fraction(4*(k-1),100)])\n",
    "            prob_ex.append(Fraction(1,int(nscen)))\n",
    "flowEx = [[Fraction(40),Fraction(24),Fraction(4)],\n",
    "        [Fraction(45),Fraction(27),Fraction(45,10)],\n",
    "        [Fraction(32),Fraction(192,10),Fraction(32,10)],\n",
    "        [Fraction(55),Fraction(33),Fraction(55,10)]]\n",
    "cost_ex = [Fraction(10),Fraction(7),Fraction(16),Fraction(6)]\n",
    "budget = Fraction(120)\n",
    "minCap = Fraction(12)\n",
    "nx = 4\n",
    "ny = 3\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We load the partition of the scenarios from the file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "partition_new = []\n",
    "for line in open('partition.txt'):\n",
    "    listWords = line.split(\" \")\n",
    "    partition_new.append([int(listWords[i]) for i in range(5,len(listWords))])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute the average and the probability of each subset of scenarios in the partition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "Aggscen_ex = []\n",
    "probs = []\n",
    "for subset in partition_new:\n",
    "    scen = [0,0,0]\n",
    "    for k in subset:\n",
    "        for j in range(ny):\n",
    "            scen[j] += scen_ex[k][j]\n",
    "    Aggscen_ex.append([Fraction(scen[j],len(subset)) for j in range(ny)])\n",
    "    probs.append(Fraction(len(subset),nscen))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we formulate the problem with the aggregated scenarios. The formulation is the following:\n",
    "$$[MP] := \\min_{x,y \\geq 0}   \\sum_{i=1}^4  cost_i x_i + \\sum_{s\\in S} p_s \\sum_{i=1}^3 \\sum_{j=1}^3flow_{ij}y_{ij}^s$$\n",
    "$$\\sum_{i=1}^4  x_i \\geq minCap$$\n",
    "$$\\sum_{i=1}^4 c_ix_i \\leq budget$$\n",
    "$$\\sum_{j=1}^3 y_{ij}^s \\leq x_i ~\\forall~ i\\in 1\\ldots 4, \\forall s\\in S$$\n",
    "$$\\sum_{i=1}^4 y_{ij}^s \\geq d_j^s  ~\\forall~ j\\in 1\\ldots 3,  \\forall s\\in S$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we are aggregating subsets of scenarios, which is equivalent to aggregate the last two constraints, so the solution is a lower bound of the problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = qsoptex.ExactProblem()\n",
    "for i in range(nx):\n",
    "    p.add_variable(name='x_'+str(i), objective = cost_ex[i], lower=0)\n",
    "for k in range(len(partition_new)):\n",
    "    for i in range(nx):\n",
    "        for j in range(ny):\n",
    "            objvar = flowEx[i][j]*probs[k]\n",
    "            p.add_variable(name='y_'+str(k)+'_'+str(i)+'_'+str(j), objective = objvar, lower=0)\n",
    "coef = {'x_'+str(i):cost_ex[i] for i in range(nx)}\n",
    "p.add_linear_constraint(qsoptex.ConstraintSense.LESS, coef, rhs=budget)    \n",
    "coef = {'x_'+str(i):1 for i in range(nx)}\n",
    "p.add_linear_constraint(qsoptex.ConstraintSense.GREATER, coef, rhs=minCap)\n",
    "for k in range(len(partition_new)):\n",
    "    for i in range(nx):\n",
    "        coef = {'x_'+str(i):-1}\n",
    "        for j in range(ny):\n",
    "            coef['y_'+str(k)+'_'+str(i)+'_'+str(j)] = Fraction(1)\n",
    "        p.add_linear_constraint(qsoptex.ConstraintSense.LESS, coef, rhs=0)\n",
    "for k in range(len(partition_new)):\n",
    "    for j in range(ny):\n",
    "        coef = {'y_'+str(k)+'_'+str(i)+'_'+str(j):1 for i in range(nx)} \n",
    "        p.add_linear_constraint(qsoptex.ConstraintSense.GREATER, coef, rhs=Aggscen_ex[k][j])\n",
    "p.set_objective_sense(qsoptex.ObjectiveSense.MINIMIZE)\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "status = p.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We solve and get a feasible solution, and a lower bound of the problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Fraction(2256294001, 10000000), 225.6294001]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LB = p.get_objective_value()\n",
    "[LB,float(LB)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Fraction(21, 25), Fraction(17, 5), Fraction(47, 25), Fraction(147, 25)]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol = [p.get_value('x_'+str(i)) for i in range(nx)]\n",
    "sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.84, 3.4, 1.88, 5.88]"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[float(p.get_value('x_'+str(i))) for i in range(nx)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we solve the second-stage problem for the solution $x=(0.84,3.4,1.88,5.88)$ for each original scenario."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: Python interface of QSopt does not provide the duals variables, so instead we formulate the dual problem, which is:\n",
    "\n",
    "$$\\max_{\\nu,\\mu \\geq 0}  \\sum_{j=1}^3 \\mu_j d_j^\\xi  -   \\sum_{i=1}^4 \\nu_i \\hat{x}_i$$\n",
    "$$\\mu_j - \\nu_i \\leq f_{ij}  \\qquad \\forall i\\in 1..4, j\\in 1..3$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "metadata": {},
   "outputs": [],
   "source": [
    "logging.basicConfig(level=0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = qsoptex.ExactProblem()\n",
    "for i in range(nx):\n",
    "    s.add_variable(name='v_'+str(i), lower=0)\n",
    "for j in range(ny):\n",
    "    s.add_variable(name='u_'+str(j), lower=0)\n",
    "for i in range(nx):\n",
    "    for j in range(ny):\n",
    "        coef = {'u_'+str(j):int(1), 'v_'+str(i):int(-1)}\n",
    "        s.add_linear_constraint(qsoptex.ConstraintSense.LESS, coef, rhs=flowEx[i][j])\n",
    "s.set_objective_sense(qsoptex.ObjectiveSense.MAXIMIZE)\n",
    "s.set_param(qsoptex.Parameter.SIMPLEX_DISPLAY, 0)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Warning: solving the subproblems can take a long time depending on your computer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "dualsSols = []\n",
    "dualsObjval = []\n",
    "for idx in range(len(partition_new)):\n",
    "    #if (idx%50 == 0 ):\n",
    "    #    print(idx)\n",
    "    #it = 0\n",
    "    for k in partition_new[idx]:\n",
    "        coef = {'v_'+str(i): -sol[i] for i in range(nx)}\n",
    "        for j in range(ny):\n",
    "            coef['u_'+str(j)] = scen_ex[k][j]\n",
    "        s.set_linear_objective(coef)\n",
    "        s.solve()\n",
    "        dual = [s.get_value('u_'+str(j)) for j in range(ny)]\n",
    "        #if (it > 0):\n",
    "        #    if (dual != dualsFin[-1:] ):\n",
    "        #        print(\"ERROR en idx=%d it=%d\" % (idx,it))\n",
    "        dualsSols.append(dual)\n",
    "        dualsObjval.append(s.get_objective_value())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we compute the upper bound, by suming all subproblems objectives and adding the cost of the first-stage solution $x$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Fraction(2256294001, 10000000), 225.6294001]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "UB = sum(dualsObjval)/1000000 + sum([sol[i]*cost_ex[i] for i in range(nx)])\n",
    "[UB, float(UB)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Fraction(2256294001, 10000000), 225.6294001]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[LB, float(LB)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the upper and lower bound are equal, this is the optimal solution of LandS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.7.7 64-bit",
   "language": "python",
   "name": "python37764bitbf82ff8733e44256b381353da7db86a3"
  },
  "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}