solver.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Benchmark the solver"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Package versions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gurobipy as gb\n",
    "from gurobipy import GRB\n",
    "\n",
    "import scipy.sparse as sp\n",
    "from scipy.optimize import linprog, milp, Bounds, LinearConstraint\n",
    "import scipy\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cvxpy version:  1.2.1\n",
      "pyomo version:  6.4.2\n",
      "numpy version:  1.23.3\n",
      "pandas version:  1.5.0\n"
     ]
    }
   ],
   "source": [
    "print(\"numpy version: \", np.__version__)\n",
    "print(\"pandas version: \", pd.__version__)\n",
    "print(\"scipy version: \", scipy.__version__)\n",
    "# print(\"gekko version: \", m.version)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Benchmark case"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "min. $ c^Tx + d $\n",
    "\n",
    "s.t.\n",
    "\n",
    "$ Gx \\leq h $\n",
    "\n",
    "$ Ax = b $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def data_prep(x_len, G_len, A_len):\n",
    "    np.random.seed(2022)\n",
    "    c = np.random.randn(x_len)\n",
    "\n",
    "    np.random.seed(2021)\n",
    "    d = np.random.randn(x_len)\n",
    "\n",
    "    np.random.seed(2020)\n",
    "    G = np.random.randn(G_len, x_len)\n",
    "\n",
    "    np.random.seed(2019)\n",
    "    h = np.random.randn(G_len)\n",
    "\n",
    "    np.random.seed(2018)\n",
    "    A = np.random.randn(A_len, x_len)\n",
    "\n",
    "    np.random.seed(2017)\n",
    "    b = np.random.randn(A_len)\n",
    "    return c, d, G, h, A, b\n",
    "\n",
    "x_len = 2000\n",
    "G_len = 4000\n",
    "A_len = 2000\n",
    "c, d, G, h, A, b = data_prep(2*x_len, G_len, A_len)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gbm1(c, d, G, h, A, b, x_len):\n",
    "    # Create a new model\n",
    "    m1 = gb.Model(\"m1\")\n",
    "    # Create variables\n",
    "    x = m1.addMVar(shape=2*x_len, vtype=GRB.BINARY, name=\"x\")\n",
    "    x[0:x_len].setAttr(\"vtype\", GRB.CONTINUOUS)\n",
    "    m1.setObjective(c @ x, GRB.MINIMIZE)\n",
    "\n",
    "    # Add constraints\n",
    "    m1.addConstr(G @ x <= h, name=\"ic\")\n",
    "    m1.addConstr(A @ x == b, name=\"ec\")\n",
    "    m1.update()\n",
    "    return m1\n",
    "\n",
    "def gbm2(c, d, G, h, A, b, x_len):\n",
    "    # Create a new model\n",
    "    m2 = gb.Model(\"m2\")\n",
    "    x1 = m2.addVars(range(x_len), vtype=GRB.CONTINUOUS, name='x1')\n",
    "    x2 = m2.addVars(range(x_len), vtype=GRB.BINARY, name='x2')\n",
    "    obj1 = sum(x1[gen] * c[gen] for gen in range(x_len))\n",
    "    obj2 = sum(x2[gen] * c[gen + x_len] for gen in range(x_len))\n",
    "    m2.setObjective(expr=obj1 + obj2, sense=gb.GRB.MINIMIZE)\n",
    "\n",
    "    # Add constraints\n",
    "    ec1 = sum(A[0, gen] * x1[gen] for gen in range(x_len))\n",
    "    ec2 = sum(A[1, gen] * x2[gen] for gen in range(x_len))\n",
    "\n",
    "    for i in range(A_len):\n",
    "        ec = sum(A[i, gen] * x1[gen] for gen in range(x_len))\n",
    "        m2.addConstr(ec == b[i], name=f'ec{i}')\n",
    "\n",
    "    for i in range(G_len):\n",
    "        ic = sum(G[i, gen] * x1[gen] for gen in range(x_len))\n",
    "        m2.addConstr(ic == h[i], name=f'ic{i}')\n",
    "    m2.update()\n",
    "    return m2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set parameter Username\n",
      "Academic license - for non-commercial use only - expires 2023-04-08\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<gurobi.Model MIP instance m1: 6000 constrs, 4000 vars, No parameter changes>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gbm1(c, d, G, h, A, b, x_len)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.25 s ± 17 ms per loop (mean ± std. dev. of 5 runs, 5 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit  -r5 -n5\n",
    "gbm1(c, d, G, h, A, b, x_len)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit  -r5 -n5\n",
    "gbm2(c, d, G, h, A, b, x_len)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve with gurobi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## gurobi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Generate a random problem\n",
    "# np.random.seed(0)\n",
    "# m, n= 40, 25\n",
    "\n",
    "# A = np.random.rand(m, n)\n",
    "# b = np.random.randn(m)\n",
    "\n",
    "# # Construct a CVXPY problem\n",
    "# x = cp.Variable(n, integer=True)\n",
    "# objective = cp.Minimize(cp.sum_squares(A @ x - b))\n",
    "# prob = cp.Problem(objective)\n",
    "# prob.solve(solver='CVXOPT')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.13 ('ams')",
   "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.13"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "d9202bbf0491b3c6145e36d5ba22687a8b95e024b0a1e787962165edd6294d6d"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}