{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Multi-commodity Min-Cost Flow Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Linear Multi-commodity min-cost flow (MMCF) problems can be modeled in a number of ways depending on how the term “commodity” is defined. \n",
    "There are three major options: a commodity may originate at a subset of nodes in the network and be destined for another subset of nodes, or it may originate from a single node and be destined for a subset of nodes, or it may originate from a single node and be destined for a single node.\n",
    "Let us consider a model of the first case.\n",
    "<br>In the MMCF problem, a directed graph $G(N,A)$, where $N$ a set of nodes is and $A$ is a set of Arcs is given. $A$ set of $K$ commodities has to be routed on $G$ at minimal total cost while satisfying the usual flow-conservation constraints at the nodes.<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<strong>Parameters</strong> <br>$c^k_{ij}$ : cost of assigning commodity $k$ to arc $ij$ <br>$d_{ij}$ : capacity of arc $ij$ <br>$b^k_{i}$ : demand/supply of commodity $k$ at node $i$ <br><strong>Decision variable</strong> <br>$x^k_{ij}$ : amount of flow of commodity $k$ to arc $ij$ <br><br>The arc-formulation according to [Cappanera and Frangioni](https://pubsonline.informs.org/doi/10.1287/ijoc.15.4.369.24887) is given as:<br>$$  \n",
    "    \\begin{array}{ll}\n",
    "    \\mbox{minimize}   & \\sum_{k\\in \\mathcal{K}}\\sum_{{ij}\\in \\mathcal{A}} c^k_{ij} x^k_{ij} \\quad\\quad\\quad(1)\\\\\n",
    "    \\mbox{subject to} & \\sum_{{ij}\\in \\mathcal{A}} x^k_{ij} - \\sum_{{ji}\\in \\mathcal{A}} x^k_{ji} = b^k_{i}\\quad  \\forall{{i}\\in \\mathcal{N}}, \\forall{{k}\\in \\mathcal{k}} \\quad(2)\\\\\n",
    "    & 0\\leq x^k_{ij} \\leq d^k_{ij} \\quad  \\forall{{ij}\\in \\mathcal{A}}, \\forall{{k}\\in \\mathcal{k}} \\quad(3)\\\\\n",
    "    & \\sum_{{k}\\in \\mathcal{K}} x^k_{ij} \\leq d_{ij} \\quad  \\forall{{ij}\\in \\mathcal{A}}, \\quad(4)\\\\\n",
    "    \\end{array}\n",
    "$$<br>Equation (1) defines the objective function and equation (2) is the flow-conservation constraint. Equations (3) and (4) are respectively the individual and mutual capacity constraints.<br>The goal here is to solve instances from [MNETGEN](http://groups.di.unipi.it/optimize/Data/MMCF.html#MNetGen) (a popular MMCF benchmark instance set) using the formulation above. Please check the link for full details on the instance generation and file formats."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following code, we introduce the various varaibles that hold our instance data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import packages\n",
    "import cvxpy as cp\n",
    "import numpy as np\n",
    "\n",
    "problemName = 'instance/64-4-1'\n",
    "nodeFile = problemName+'.nod'\n",
    "arcFile = problemName+'.arc'\n",
    "mutFile = problemName+'.mut'\n",
    "supFile = problemName+'.sup'\n",
    "instanceInfo = [] # [commodity, nodes, arcs, cap_arcs]\n",
    "numNodes = 0\n",
    "numCommodity = 0\n",
    "numArcs = 0\n",
    "numCapacitatedArcs = 0\n",
    "bigM = 99999999"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above code assumes the instance files are stored in a folder named instances. Solving a single instance requires 4 files with details of the file formats at [MNETGEN](http://groups.di.unipi.it/optimize/Data/MMCF.html#MNetGen).<br>In the next code, we read the each of the 4 files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Read the node file\n",
    "nodeFileData = open(nodeFile, \"r\")\n",
    "\n",
    "for eachLine in nodeFileData:\n",
    "    instanceInfo.append(int(eachLine))\n",
    "numCommodity = instanceInfo[0]\n",
    "numNodes = instanceInfo[1]\n",
    "numArcs = instanceInfo[2]\n",
    "nodeFileData.close()\n",
    "\n",
    "#Read mutual capacity pointer file\n",
    "mutualCapacityPointer = np.loadtxt(mutFile)\n",
    "\n",
    "mcpFileData.close()\n",
    "\n",
    "#Read arc file\n",
    "uniqueArcs = {}\n",
    "arcCapacity = {}\n",
    "mutualCapacity = {}\n",
    "cost = {}\n",
    "temp_arc = [0]*numCommodity\n",
    "temp_cost = [0]*numCommodity\n",
    "arcFileData = open(arcFile, \"r\")\n",
    "\n",
    "for eachLine in arcFileData:\n",
    "    data = eachLine.split(\"\t\")\n",
    "    uniqueArcs[int(data[0])] = [int(data[1]),int(data[2])]\n",
    "    mutualCapacity[int(data[0])] = int(data[6])\n",
    "    inner_index = int(data[3])-1\n",
    "    outer_index = int(data[0])-1\n",
    "    if temp_arc[inner_index]==0:\n",
    "        temp_arc[inner_index] = [0]*numArcs\n",
    "        temp_arc[inner_index][outer_index]=int(data[5])\n",
    "    else:\n",
    "        temp_arc[inner_index][outer_index]=int(data[5])\n",
    "    if temp_cost[inner_index]==0:\n",
    "        temp_cost[inner_index] = [bigM]*numArcs\n",
    "        temp_cost[inner_index][outer_index] = float(data[4])\n",
    "    else:\n",
    "        temp_cost[inner_index][outer_index] = float(data[4])\n",
    "    \n",
    "\n",
    "for i in range(numCommodity):\n",
    "    arcCapacity[i+1] = temp_arc[i]    \n",
    "    cost[i+1] = temp_cost[i] \n",
    "    \n",
    "arcFileData.close()\n",
    "\n",
    "\n",
    "#Read sup file\n",
    "supplyDemand = {}\n",
    "temp_supDem = [0]*numCommodity\n",
    "supFileData = open(supFile, \"r\")\n",
    "\n",
    "\n",
    "for eachLine in supFileData:\n",
    "    data = list(map(int, eachLine.split(\"\t\")))\n",
    "    inner_index = data[1]-1\n",
    "    outer_index = data[0]-1\n",
    "    if temp_supDem[inner_index]==0:\n",
    "        temp_supDem[inner_index] = [0]*numNodes\n",
    "        temp_supDem[inner_index][outer_index]= data[2]\n",
    "    else:\n",
    "        temp_supDem[inner_index][outer_index]= data[2]\n",
    "\n",
    "for i in range(numCommodity):\n",
    "    supplyDemand[i+1] = temp_supDem[i]  \n",
    "    \n",
    "supFileData.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We skip the details of the above code to read the problem data since there are other ways for it to be accomplished.<br>In the final section, we introduce the model and solve the MMCF problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MMCF objective function value is 290806.3029314588\n"
     ]
    }
   ],
   "source": [
    "#---------------------------Solve the math model------------------------------\n",
    "#\n",
    "C = [(a,k) for a in uniqueArcs.keys() for k in range(numCommodity)] #decision variable indexer\n",
    "X_key = {}\n",
    "my_ind = 0\n",
    "for (a,k) in C:\n",
    "    X_key[a,k] = my_ind\n",
    "    my_ind+=1\n",
    "X = cp.Variable(shape=len(C), nonneg=True, name=\"X\") \n",
    "\n",
    "allConstraints = []\n",
    "#constraint 1\n",
    "for i in range(numNodes):\n",
    "    for k in range(numCommodity):\n",
    "        sumIn = 0\n",
    "        sumOut = 0\n",
    "        for akey,value in uniqueArcs.items():\n",
    "            if value[0]==(i+1):\n",
    "                sumIn +=X[X_key[akey,k]]\n",
    "            elif value[1]==(i+1):\n",
    "                sumOut +=X[X_key[akey,k]]\n",
    "        const1=sumIn-sumOut==supplyDemand[k+1][i]\n",
    "        allConstraints.append(const1)\n",
    " \n",
    "\n",
    "#constraint 2\n",
    "for (a,k) in C:\n",
    "    maxVal = bigM#sys.float_info.max\n",
    "    if arcCapacity[k+1][a-1]==-1:\n",
    "        const2=X[X_key[a,k]]<=maxVal\n",
    "        allConstraints.append(const2)\n",
    "    else:\n",
    "        const2=X[X_key[a,k]]<=arcCapacity[k+1][a-1]\n",
    "        allConstraints.append(const2)\n",
    "        \n",
    "#constraint 3\n",
    "for a in uniqueArcs.keys():\n",
    "    if mutualCapacity[a]!=0:\n",
    "        totalVars = [X[X_key[a,k]] for k in range(numCommodity)]\n",
    "        lhs = cp.sum(totalVars)\n",
    "        const3=lhs<=mutualCapacityPointer[mutualCapacity[a]-1][1]\n",
    "        allConstraints.append(const3)\n",
    "\n",
    "actualCost = np.asarray([cost[k+1][a-1] for (a,k) in C])\n",
    "objExpr = cp.sum(actualCost@X) #objective function\n",
    "objFunc = cp.Minimize(objExpr)\n",
    "prob = cp.Problem(objFunc, allConstraints)\n",
    "prob.solve()\n",
    "\n",
    "print(\"MMCF objective function value is {}\".format(prob.value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is evident from the math model code that the $x^k_{ij}$ was first deflated to $x^k_{a}$ then to $x$. To emphasize, the change to 1 dimension is to demonstrate a way to get pass the 2 dimension limitation in CVXPY as noted in the manual.To achieve this, we build an indexer for $x$ and then store the mapping of the 2-d format in a dictionary. This way the code is still readable as in the math model and it shows how already coded math models in other platforms can be ported to CVXPY."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
