simplex.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The simplex method "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook executes and shows each step of the simplex method using Dantzig's rule as described in:\n",
    "    \n",
    "&nbsp;&nbsp;&nbsp;&nbsp;D. J. Rader (2010). _Deterministic Operations Research: Models and Algorithms_. Wiley."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Setup code.__ You only need to run this code cell once."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import symbols, solve, latex \n",
    "from sympy.matrices import Matrix, zeros\n",
    "from sympy.interactive import printing\n",
    "from IPython.display import Markdown, display\n",
    "printing.init_printing(use_latex=\"mathjax\")\n",
    "\n",
    "def printmd(string):\n",
    "    display(Markdown(string))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Input the variable names of the canonical LP you want to solve.__ Make sure you change the variable names in all three places."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_1, x_2, x_3, x_4, x_5 = symbols('x_1, x_2, x_3, x_4, x_5')\n",
    "variables = [x_1, x_2, x_3, x_4, x_5] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Input the constraint matrix $\\mathbf{A}$, the constraint vector $\\mathbf{b}$, and the objective function vector $\\mathbf{c}$ of the canonical LP you want to solve.__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = Matrix([[1, 5, 0, 1, 0],\n",
    "            [-2, 3, 4, 0, 1]])\n",
    "b = Matrix([1, 2])\n",
    "c = Matrix([3, 2, 1, 0, 0.25])\n",
    "maximize = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Input the initial basis.__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "basis = {x_2, x_3}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Run the simplex method!__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Determine m and n \n",
    "(m, n) = A.shape\n",
    "\n",
    "# Compute initial BFS\n",
    "x = zeros(n, 1)\n",
    "for i, var in enumerate(variables):\n",
    "    if var in basis:\n",
    "        x[i] = var\n",
    "\n",
    "soln = solve(A*x - b)\n",
    "bfs = zeros(n, 1)\n",
    "for i, var in enumerate(variables):\n",
    "    if var in basis:\n",
    "        bfs[i] = soln[var]\n",
    "\n",
    "# Print the initial BFS and basis\n",
    "printmd(\"The intial BFS and basis are\")\n",
    "printmd(\"&nbsp;&nbsp;&nbsp;&nbsp;$\\mathbf{{x}}^{{0}} = {0}$&nbsp;&nbsp;&nbsp;&nbsp;$\\mathcal{{B}}^{{0}} = {1}$\".format(latex(bfs.T), latex(basis)))\n",
    "\n",
    "# Start main loop\n",
    "iteration = 0\n",
    "optimal_solution_found = False\n",
    "unbounded = False\n",
    "\n",
    "while True:\n",
    "    iteration += 1\n",
    "    printmd(\"**Iteration {0}**\".format(iteration))\n",
    "    \n",
    "    # Compute the nonbasic variables\n",
    "    # Sort the nonbasic variables according to the original variable order\n",
    "    #   for more sensible printing\n",
    "    nonbasic_variables = list(set(variables) - basis)\n",
    "    nonbasic_variables.sort(key=lambda x: variables.index(x))\n",
    "    \n",
    "    # Compute simplex directions and reduced costs\n",
    "    direction = {}\n",
    "    reduced_cost = {}\n",
    "    for nb in nonbasic_variables:\n",
    "        # Set up simplex direction template\n",
    "        d = zeros(n, 1)\n",
    "        for i, var in enumerate(variables):\n",
    "            if var in basis:\n",
    "                d[i] = var\n",
    "            elif var == nb:\n",
    "                d[i] = 1\n",
    "        \n",
    "        # Solve for simplex direction\n",
    "        soln = solve(A*d)\n",
    "        direction[nb] = zeros(n,1)\n",
    "        for i, var in enumerate(variables):\n",
    "            if var in basis:\n",
    "                direction[nb][i] = soln[var]\n",
    "            elif var == nb:\n",
    "                direction[nb][i] = 1\n",
    "        \n",
    "        # Compute reduced cost\n",
    "        reduced_cost[nb] = c.dot(direction[nb])\n",
    "        \n",
    "    # Print simplex directions and associated reduced costs\n",
    "    printmd(\"&nbsp;&nbsp;&nbsp;&nbsp;The simplex directions are:\")\n",
    "    for nb in nonbasic_variables:\n",
    "        printmd(\"&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\\mathbf{{d}}^{{{0}}} = {1}$&nbsp;&nbsp;&nbsp;&nbsp;$\\\\bar{{c}}_{{{0}}} = {2}$\".format(nb, latex(direction[nb].T), reduced_cost[nb]))\n",
    "        \n",
    "    # Choose entering variable according to Dantzig's rule\n",
    "    if maximize:\n",
    "        (entering, max_reduced_cost) = max(reduced_cost.items(), key=(lambda x: x[1]))\n",
    "        if max_reduced_cost <= 0:\n",
    "            optimal_solution_found = True\n",
    "            break\n",
    "    else:\n",
    "        (entering, min_reduced_cost) = min(reduced_cost.items(), key=(lambda x: x[1]))\n",
    "        if min_reduced_cost >= 0:\n",
    "            optimal_solution_found = True\n",
    "            break\n",
    "    \n",
    "    # Print entering variable\n",
    "    printmd(\"&nbsp;&nbsp;&nbsp;&nbsp;Using Dantzig's rule, the entering variable is ${0}$.\".format(latex(entering)))\n",
    "            \n",
    "    # Check for unboundedness\n",
    "    if all(direction[entering][j] >= 0 for j, var in enumerate(variables)):\n",
    "        unbounded = True\n",
    "        break\n",
    "    \n",
    "    # Compute step size and leaving variable using MRT\n",
    "    ratio = {}\n",
    "    for j, var in enumerate(variables):\n",
    "        if direction[entering][j] < 0:\n",
    "            ratio[var] = bfs[j] / -direction[entering][j]\n",
    "    (leaving, lambda_max) = min(ratio.items(), key=(lambda x: x[1]))\n",
    "    \n",
    "    # Print the step size and leaving variable\n",
    "    printmd(\"&nbsp;&nbsp;&nbsp;&nbsp;Using the minimum ratio test, the step size is $\\lambda_{{\\max}} = {0}$ and the leaving variable is ${1}$.\".format(lambda_max, latex(leaving)))\n",
    "        \n",
    "    # Update solution and basis\n",
    "    bfs = bfs + lambda_max * direction[entering]\n",
    "    basis.add(entering)\n",
    "    basis.remove(leaving)\n",
    "    \n",
    "    # Print the next BFS and basis\n",
    "    printmd(\"&nbsp;&nbsp;&nbsp;&nbsp;$\\mathbf{{x}}^{{{0}}} = {1}$&nbsp;&nbsp;&nbsp;&nbsp;$\\mathcal{{B}}^{{{0}}} = {2}$\".format(iteration, latex(bfs.T), latex(basis)))\n",
    "\n",
    "# Print the reason for termination\n",
    "if optimal_solution_found:\n",
    "    printmd(\"&nbsp;&nbsp;&nbsp;&nbsp;Since there are no improving simplex directions, the BFS $\\mathbf{{x}}^{{{0}}} = {1}$ is optimal with value {2}.\".format(iteration - 1, latex(bfs.T), c.dot(bfs)))\n",
    "elif unbounded:\n",
    "    printmd(\"&nbsp;&nbsp;&nbsp;&nbsp;Since $\\mathbf{{d}}^{{{0}}}$ is an improving simplex direction with all nonnegative components, the LP is unbounded.\".format(latex(entering)))\n",
    "else:\n",
    "    print(\"ERROR. Something went wrong: algorithm terminated, but optimal solution or certificate of unboundedness were not found\")"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}