mip/cutting_stock/cutting_stock.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cutting stock problem with column generation"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Restricted master problem\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    \\text{min} \\quad & \\sum_{p \\in P} c_{p} x_{p} \\\\\n",
    "    \\text{s.t} \\quad & d_{i} \\leq \\sum_{p \\in P} a_{i, p} x_{p} & \\forall \\; i \\in I \\\\\n",
    "    & x_{p} \\geq 0 & \\forall \\; p \\in P \\\\\n",
    "    & x_{p} \\in \\mathbb{Z} & \\forall \\; p \\in P\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pricing problem\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    \\text{max} \\quad & \\sum_{i \\in I} \\pi_{i} y_{i} \\\\\n",
    "    \\text{s.t} \\quad & \\sum_{i \\in I} w_{i} y_{i} \\leq W\\\\\n",
    "    & y_{i} \\geq 0 & \\forall \\; i \\in I \\\\\n",
    "    & y_{i} \\in \\mathbb{Z} & \\forall \\; i \\in I\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.optimize import linprog\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      w   d\n",
      "0  75.0  38\n",
      "1  71.3  44\n",
      "2  67.0  30\n",
      "3  61.4  41\n",
      "4  58.0  36\n"
     ]
    }
   ],
   "source": [
    "dataset = pd.read_csv(\"data.txt\", sep=\" \")\n",
    "print(dataset.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_knapsack(W, w, duals):\n",
    "    return linprog(-duals, A_ub=np.atleast_2d(w), b_ub=np.atleast_1d(W), bounds=(0, np.inf), integrality=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compare duality difference: 0.0\n",
      "Iteration: 0; Reduced cost: -0.591\n",
      "Iteration: 1; Reduced cost: -0.535\n",
      "Iteration: 2; Reduced cost: -0.533\n",
      "Iteration: 3; Reduced cost: -0.500\n",
      "Iteration: 4; Reduced cost: -0.515\n",
      "Iteration: 5; Reduced cost: -0.500\n",
      "Iteration: 6; Reduced cost: -0.511\n",
      "Iteration: 7; Reduced cost: -0.500\n",
      "Iteration: 8; Reduced cost: -0.400\n",
      "Iteration: 9; Reduced cost: -0.413\n",
      "Iteration: 10; Reduced cost: -0.433\n",
      "Iteration: 11; Reduced cost: -0.333\n",
      "Rounding solution 339.0\n",
      "Integer solution: 334.0\n"
     ]
    }
   ],
   "source": [
    "# Total width\n",
    "W = 100.0\n",
    "\n",
    "# Width and amount associated with each demand\n",
    "w = dataset.w.values\n",
    "d = dataset.d.values\n",
    "\n",
    "# LP parameters\n",
    "A = np.eye(dataset.shape[0]) * (W // w)\n",
    "c = np.ones_like(w)\n",
    "\n",
    "# Initial solution\n",
    "sol = linprog(c, A_ub=-A, b_ub=-d, bounds=(0, None))\n",
    "sol_dual = linprog(-d, A_ub=A.T, b_ub=c, bounds=(0, None))\n",
    "diff = np.abs(sol_dual.x + sol.ineqlin.marginals).sum()\n",
    "print(f\"Compare duality difference: {diff}\")\n",
    "\n",
    "# Iterate\n",
    "for _ in range(1000):\n",
    "    duals = -sol.ineqlin.marginals\n",
    "    price_sol = solve_knapsack(W, w, duals)\n",
    "    y = price_sol.x\n",
    "    if 1 + price_sol.fun < -1e-4:\n",
    "        print(f\"Iteration: {_}; Reduced cost: {(1 + price_sol.fun):.3f}\")\n",
    "        A = np.hstack((A, y.reshape((-1, 1))))\n",
    "        c = np.append(c, 1)\n",
    "        sol = linprog(c, A_ub=-A, b_ub=-d, bounds=(0, None))\n",
    "    else:\n",
    "        break\n",
    "\n",
    "sol_round = linprog(c, A_ub=-A, b_ub=-d, bounds=(0, np.inf), integrality=0)\n",
    "print(f\"Rounding solution {np.ceil(sol_round.x).sum()}\")\n",
    "sol = linprog(c, A_ub=-A, b_ub=-d, bounds=(0, np.inf), integrality=1)\n",
    "print(f\"Integer solution: {sol.x.sum()}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcwAAAEiCAYAAAB0leJ0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAFpElEQVR4nO3dPXKbYABFUZxxoZRZSorsKsvJrlJ4KSntwjOkSRlZ1wgJPnFOjQeEfu7Q+D3N8zxPAMCHvmx9AQAwAsEEgEAwASAQTAAIBBMAAsEEgEAwASAQTAAIBBMAgud64NfvPxed4M/vX4v+DgDu5RRq6AkTAALBBIBAMAEgEEwACAQTAALBBIBAMAEgEEwACAQTAALBBIBAMAEgEEwACAQTAIKneZ7ncuDb+7ITfPth5QSAfbNWAgArEUwACAQTAALBBIBAMAEgEEwACAQTAALBBIBAMAEgEEwACAQTAALBBIBAMAEguPlayVJLV06mydIJAJ9jrQQAViKYABAIJgAEggkAgWACQCCYABAIJgAEggkAgWACQCCYABAIJgAEggkAgWACQCCYABDsdt7rGkunwcyCARyTeS8AWIlgAkAgmAAQCCYABIIJAIFgAkAgmAAQCCYABIIJAIFgAkAgmAAQCCYABIIJAMFDrpUsZeUEYHtLf4uv8fpy+XfcEyYABIIJAIFgAkAgmAAQCCYABIIJAIFgAkAgmAAQCCYABIIJAIFgAkAgmAAQCCYABNZKVmDlBBiB36rzTs+Xj/GECQCBYAJAIJgAEAgmAASCCQCBYAJAIJgAEAgmAASCCQCBYAJAIJgAEAgmAASCCQCBtZINLV0OmKZjrAfACLZYABlldWSU65wmayUAsBrBBIBAMAEgEEwACAQTAALBBIBAMAEgEEwACAQTAALBBIBAMAEgEEwACAQTAAJrJYMaaQUAPmuLJR/fqfOOcG+slQDASgQTAALBBIBAMAEgEEwACAQTAALBBIBAMAEgEEwACAQTAALBBIBAMAEgEEwACAQTAALzXgdzhJke9mOLz9son/EtJsw4z7wXAKxEMAEgEEwACAQTAALBBIBAMAEgEEwACAQTAALBBIBAMAEgEEwACAQTAALBBIDAWgnJKAsQwH6MtFbz+nL5nJ4wASAQTAAIBBMAAsEEgEAwASAQTAAIBBMAAsEEgEAwASAQTAAIBBMAAsEEgEAwASCwVsJNWTmB/1v63Zim+38/jvA9Pj1fPsYTJgAEggkAgWACQCCYABAIJgAEggkAgWACQCCYABAIJgAEggkAgWACQCCYABAIJgAE1krYpZGWHNiXaz4797TF5/QIqyNLWSsBgJUIJgAEggkAgWACQCCYABAIJgAEggkAgWACQCCYABAIJgAEggkAgWACQCCYABBYK+HhWGTYD+/Fbbiv67NWAgArEUwACAQTAALBBIBAMAEgEEwACAQTAALBBIBAMAEgEEwACAQTAALBBIBAMAEgEEwACMx7wT8mk47t3u//0vNdc85RbHFvzHsBwEoEEwACwQSAQDABIBBMAAgEEwACwQSAQDABIBBMAAgEEwACwQSAQDABIBBMAAislcCVrJysz5IH92atBABWIpgAEAgmAASCCQCBYAJAIJgAEAgmAASCCQCBYAJAIJgAEAgmAASCCQCBYAJAEP4/O/CRpesYFjl4dI+25OMJEwACwQSAQDABIBBMAAgEEwACwQSAQDABIBBMAAgEEwACwQSAQDABIBBMAAgEEwCCp3me53Lg2/utLwWoRlmBuGaRZalRXuNeFzn2YIsln1PY7vKECQCBYAJAIJgAEAgmAASCCQCBYAJAIJgAEAgmAASCCQCBYAJAIJgAEAgmAASCCQCBYAJAYN4LDsQU1bF5/88z7wUAKxFMAAgEEwACwQSAQDABIBBMAAgEEwACwQSAQDABIBBMAAgEEwACwQSAQDABILBWAlw00srFSNe61L1f4xHuqbUSAFiJYAJAIJgAEAgmAASCCQCBYAJAIJgAEAgmAASCCQCBYAJAIJgAEAgmAASCCQBBXisBgCPzhAkAgWACQCCYABAIJgAEggkAgWACQCCYABAIJgAEggkAwV+AzQhHonEGzAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 700x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=[7, 3], dpi=100)\n",
    "hmap = ax.imshow(A > 1e-6, cmap=\"Blues\")\n",
    "plt.axis('off')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv",
   "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.9.7"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}