benders.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "587a65e5-6be9-4b24-8b39-c7123d576e21",
   "metadata": {},
   "outputs": [],
   "source": [
    "# using problem is here: https://homepages.rpi.edu/~mitchj/matp6620/benders_eg/benders_eg_slides.html"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89eb5cf6-dc30-4948-abfa-39037ed473ac",
   "metadata": {},
   "source": [
    "# Original Problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "88c78ad0-7d24-4fee-9b0a-42c715031010",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Restricted license - for non-production use only - expires 2024-10-28\n",
      "Gurobi Optimizer version 10.0.2 build v10.0.2rc0 (mac64[x86])\n",
      "\n",
      "CPU model: Intel(R) Core(TM) i7-8559U CPU @ 2.70GHz\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "\n",
      "Optimize a model with 9 rows, 6 columns and 18 nonzeros\n",
      "Model fingerprint: 0xac069ec6\n",
      "Variable types: 4 continuous, 2 integer (2 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+00]\n",
      "  Objective range  [5e+00, 2e+01]\n",
      "  Bounds range     [1e+00, 1e+00]\n",
      "  RHS range        [1e+00, 1e+00]\n",
      "Found heuristic solution: objective 2.0000000\n",
      "Presolve removed 3 rows and 1 columns\n",
      "Presolve time: 0.00s\n",
      "Presolved: 6 rows, 5 columns, 12 nonzeros\n",
      "Variable types: 4 continuous, 1 integer (1 binary)\n",
      "\n",
      "Root relaxation: cutoff, 0 iterations, 0.00 seconds (0.00 work units)\n",
      "\n",
      "Explored 1 nodes (0 simplex iterations) in 0.01 seconds (0.00 work units)\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 1: 2 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 2.000000000000e+00, best bound 2.000000000000e+00, gap 0.0000%\n",
      "x[0] = 1.0\n",
      "x[1] = 1.0\n",
      "x[2] = 0.0\n",
      "x[3] = 0.0\n",
      "y[0] = 1.0\n",
      "y[1] = 0.0\n",
      "Obj: 2.0\n"
     ]
    }
   ],
   "source": [
    "from gurobipy import Model, GRB\n",
    "\n",
    "# Create a new model\n",
    "m = Model(\"mip_problem\")\n",
    "\n",
    "# Create variables\n",
    "x = m.addVars(4, vtype=GRB.CONTINUOUS, name=\"x\")\n",
    "y = m.addVars(2, vtype=GRB.BINARY, name=\"y\")\n",
    "\n",
    "# Set objective\n",
    "m.setObjective(8*x[0] + 9*x[1] + 5*x[2] + 6*x[3] - 15*y[0] - 10*y[1], GRB.MAXIMIZE)\n",
    "\n",
    "# Add constraints\n",
    "m.addConstr(x[0] + x[2] <= 1)\n",
    "m.addConstr(x[0] + x[3] <= 1)\n",
    "m.addConstr(x[1] + x[2] <= 1)\n",
    "m.addConstr(x[1] + x[3] <= 1)\n",
    "m.addConstr(x[0] - y[0] <= 0)\n",
    "m.addConstr(x[1] - y[0] <= 0)\n",
    "m.addConstr(x[2] - y[1] <= 0)\n",
    "m.addConstr(x[3] - y[1] <= 0)\n",
    "m.addConstr(-y[0] - y[1] <= -1)\n",
    "\n",
    "# Model update\n",
    "m.update()\n",
    "\n",
    "# Optimize the model\n",
    "m.optimize()\n",
    "\n",
    "# Print the solution\n",
    "if m.status == GRB.Status.OPTIMAL:\n",
    "    for v in m.getVars():\n",
    "        print(f'{v.varName} = {v.x}')\n",
    "    print(f'Obj: {m.objVal}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2062c18d-b012-4bd3-aec1-bff0bc5ebed6",
   "metadata": {},
   "source": [
    "# Benders Decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "19e9b616-9fd3-42d9-88c8-80495737f002",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Round 0\n",
      "  opt_y    = {0: 0, 1: 0}\n",
      "  main_obj = -inf\n",
      "  sub_obj  = 0.0\n",
      "  lambda   = [-2.0, -1.0]\n",
      "  add <gurobi.TempConstr: t >= -2.0 y[0] + -1.0 y[1]>\n",
      "Round 1\n",
      "  opt_y    = {0: 1.0, 1: 1.0}\n",
      "  main_obj = -3.0\n",
      "  sub_obj  = 8.0\n",
      "  lambda   = [15.0, 10.0]\n",
      "  add <gurobi.TempConstr: t >= -17.0 + 15.0 y[0] + 10.0 y[1]>\n",
      "Round 2\n",
      "  opt_y    = {0: 1.0, 1: 0.0}\n",
      "  main_obj = -2.0\n",
      "  sub_obj  = -2.0\n",
      "  lambda   = [0.0, -1.0]\n",
      "  add <gurobi.TempConstr: t >= -2.0 + 0.0 y[0] + -1.0 y[1]>\n",
      "\n",
      "Result\n",
      "  x[0] = 1.0\n",
      "  x[1] = 1.0\n",
      "  x[2] = 0.0\n",
      "  x[3] = 0.0\n",
      "  y[0] = 1.0\n",
      "  y[1] = 0.0\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "from gurobipy import Model, GRB, quicksum\n",
    "\n",
    "# Creat main problem\n",
    "main = Model(\"main\")\n",
    "main_y = main.addVars(2, vtype=GRB.BINARY, name=\"y\")\n",
    "main_t = main.addVar(lb=-GRB.INFINITY, vtype=GRB.CONTINUOUS, name=\"t\")\n",
    "main.addConstr(- main_y[0] - main_y[1] <= -1)\n",
    "main.setObjective(main_t)\n",
    "main.Params.OutputFlag = 0 # silent/verbose mode\n",
    "\n",
    "# Create sub problem\n",
    "sub = Model(\"sub\")\n",
    "x = sub.addVars(4, vtype=GRB.CONTINUOUS, name=\"x\")\n",
    "y = sub.addVars(2, lb=0, ub=1, vtype=GRB.CONTINUOUS, name=\"y\")\n",
    "sub.addConstr(x[0] + x[2] <= 1)\n",
    "sub.addConstr(x[0] + x[3] <= 1)\n",
    "sub.addConstr(x[1] + x[2] <= 1)\n",
    "sub.addConstr(x[1] + x[3] <= 1)\n",
    "sub.addConstr(x[0] - y[0] <= 0)\n",
    "sub.addConstr(x[1] - y[0] <= 0)\n",
    "sub.addConstr(x[2] - y[1] <= 0)\n",
    "sub.addConstr(x[3] - y[1] <= 0)\n",
    "d0 = sub.addConstr(0 <= 0)  # dummy\n",
    "d1 = sub.addConstr(0 <= 0)  # dummy\n",
    "sub.setObjective(-8*x[0] - 9*x[1] - 5*x[2] - 6*x[3] + 15*y[0] + 10*y[1])\n",
    "sub.Params.OutputFlag = 0 # silent/verbose mode\n",
    "\n",
    "# Set initial solution randomly\n",
    "opt_y = {i: 0 for i in range(2)}\n",
    "main_obj = float(\"-inf\")\n",
    "\n",
    "# Benders decomposition\n",
    "index = 0\n",
    "while True:\n",
    "    # Update sub problem\n",
    "    sub.remove([d0, d1])\n",
    "    d0 = sub.addConstr(y[0] == opt_y[0])\n",
    "    d1 = sub.addConstr(y[1] == opt_y[1])\n",
    "    sub.update()\n",
    "\n",
    "    # Solve sub problem\n",
    "    sub.optimize()\n",
    "    sub_obj = sub.ObjVal\n",
    "    l = [d0.Pi, d1.Pi]  # get dual variables\n",
    "\n",
    "    # Output logging\n",
    "    print(f\"Round {index}\")\n",
    "    print(f\"  opt_y    = {opt_y}\")\n",
    "    print(f\"  main_obj = {main_obj}\")\n",
    "    print(f\"  sub_obj  = {sub_obj}\")\n",
    "    print(f\"  lambda   = {l}\")\n",
    "    \n",
    "    # Add constraint to the main problem\n",
    "    c = main.addConstr(main_t >= sub_obj + l[0] * (main_y[0] - opt_y[0]) + l[1] * (main_y[1] - opt_y[1]))\n",
    "    main.update()\n",
    "    print(f\"  add {main_t >= sub_obj + quicksum(l[i] * (main_y[i] - opt_y[i]) for i in range(2))}\")\n",
    "\n",
    "    # Solve main problem\n",
    "    main.optimize()\n",
    "    main_obj = main.ObjVal\n",
    "    opt_y = {i: main_y[i].X for i in range(2)}\n",
    "\n",
    "    # Check termination condition\n",
    "    if math.isclose(main_obj, sub_obj, abs_tol=1e-6):\n",
    "        break\n",
    "    \n",
    "    index += 1\n",
    "\n",
    "print()\n",
    "print(\"Result\")\n",
    "for v in sub.getVars():\n",
    "    print(f'  {v.varName} = {v.x}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb6b9031-5de5-414b-8792-a38a9c35a51f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}