Simplex-Matrix.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation of The Simplex Algorithm via Matrix Operations\n",
    "\n",
    "The following code implements the Simplex method with matrix operations, as opposed to the tableau method.  \n",
    "\n",
    "We begin by writing out a constrained optimization problem in *standard form* below.  The matrix $A$ holds the coefficients of the inequality constraints, the vector $b$ is the vector of solutions, and the vector $c$ holds the coefficients of the variables of the objective function that is being optimized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# User Defined Input\n",
    "\n",
    "# Example Input\n",
    "A = np.array([[-2, 1, 1, 0, 0],\n",
    "             [-1, 2, 0, 1, 0],\n",
    "             [1, 0, 0, 0, 1]])\n",
    "\n",
    "b = np.array([2, 7, 3])\n",
    "\n",
    "c = np.array([-1, -2, 0, 0, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we continue to establish the function ``Simplex`` that solves a linear constrained optimization problem using a matrix method implementation of the Simplex Algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([-2, -1,  0]),\n",
       " [1, 0, 2],\n",
       " array([0, 0]),\n",
       " [3, 4],\n",
       " array([5., 3., 3.]),\n",
       " array([1., 2.]))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def Simplex(A, b, c):\n",
    "    '''Takes input vars, computs corresponding values,\n",
    "    then uses while loop to iterate until a basic optimal solution is reached.\n",
    "    RETURNS: cbT, cbIndx, cnT, cnIndx, bHat, cnHat.\n",
    "    cbT, cbIndex is final basic variable values, and indices\n",
    "    cnT, cnIndex is final nonbasic variable values and indices\n",
    "    bHat is final solution values, \n",
    "    cnHat is optimality condition'''\n",
    "    \n",
    "    #sizes of basic and nonbasic vectors\n",
    "    basicSize = A.shape[0] # number of constraints, m\n",
    "    nonbasicSize = A.shape[1] - basicSize #n-m, number of variables\n",
    "        \n",
    "    # global index tracker of variables of basic and nonbasic variables (objective)\n",
    "    # that is, index 0 corresponds with x_0, 1 with x_1 and so on.  So each index corresponds with a variable\n",
    "    cindx = [i for i in range(0, len(c))]\n",
    "    \n",
    "    #basic variable coefficients\n",
    "    cbT = np.array(c[nonbasicSize:])\n",
    "\n",
    "    #nonbasic variable coefficients\n",
    "    cnT = np.array(c[:nonbasicSize])\n",
    "    \n",
    "    # run core simplex method until reach the optimal solution\n",
    "    while True:\n",
    "        \n",
    "        # keep track of current indices of basic and non-basic variables\n",
    "        cbIndx = cindx[nonbasicSize:]\n",
    "        cnIndx = cindx[:nonbasicSize]\n",
    "        \n",
    "        # basis matrix\n",
    "        B = A[:, cbIndx]\n",
    "        Binv = np.linalg.inv(B)\n",
    "        \n",
    "        # nonbasic variable matrix\n",
    "        N = A[:, cnIndx]\n",
    "        \n",
    "        # bHat, the values of the basic variables\n",
    "        # recall that at the start the basic variables are the slack variables, and \n",
    "        # have values equal the vector b (as primary variables are set to 0 at the start)\n",
    "        bHat = Binv @ b\n",
    "        yT = cbT @ Binv\n",
    "        \n",
    "        # use to check for optimality, determine variable to enter basis\n",
    "        cnHat = cnT - (yT @ N)\n",
    "        \n",
    "        # find indx of minimum value of cnhat, this is the variable to enter the basis\n",
    "        cnMinIndx = np.argmin(cnHat)\n",
    "\n",
    "        # break out of loop, returning values if all values of cnhat are above 0\n",
    "        if(all(i>=0 for i in cnHat)):\n",
    "            # use cbIndx to get index values of variables in bHat, and the corresponding index\n",
    "            # values in bHat are the final solution values for each of the corresponding variables\n",
    "            # ie value 0 in dbIndx corresponds with first variable, so whatever the index for the 0 is\n",
    "            # is the index in bHat that has the solution value for that variable.\n",
    "            return cbT, cbIndx, cnT, cnIndx, bHat, cnHat\n",
    "        \n",
    "        # this is the index for the column of coeffs in a for the given variable\n",
    "        indx = cindx[cnMinIndx]\n",
    "\n",
    "        Ahat = Binv @ A[:, indx]\n",
    "        \n",
    "        # now we want to iterate through Ahat and bHat and pick the minimum ratios\n",
    "        # only take ratios of variables with Ahat_i values greater than 0\n",
    "        # pick smallest ratio to get variable that will become nonbasic.\n",
    "        ratios = []\n",
    "        for i in range(0, len(bHat)):\n",
    "            Aval = Ahat[i]\n",
    "            Bval = bHat[i]\n",
    "\n",
    "            # don't look at ratios with val less then or eqaul to 0, append to keep index\n",
    "            if(Aval <= 0):\n",
    "                ratios.append(10000000)\n",
    "                continue\n",
    "            ratios.append(Bval / Aval)\n",
    "\n",
    "        ratioMinIndx = np.argmin(ratios)\n",
    "\n",
    "        #switch basic and nonbasic variables using the indices.\n",
    "        cnT[cnMinIndx], cbT[ratioMinIndx] = cbT[ratioMinIndx], cnT[cnMinIndx]\n",
    "        # switch global index tracker indices\n",
    "        cindx[cnMinIndx], cindx[ratioMinIndx + nonbasicSize] = cindx[ratioMinIndx + nonbasicSize], cindx[cnMinIndx]\n",
    "        # now repeat the loop\n",
    "        \n",
    "\n",
    "Simplex(A, b, c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following we proceed to test the function with different constrained optimization problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([-2,  0, -3]),\n",
       " [1, 4, 0],\n",
       " array([0, 0]),\n",
       " [3, 2],\n",
       " array([12.,  3.,  3.]),\n",
       " array([0.25, 1.25]))"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# example test\n",
    "A = np.array([[2, 1, 1, 0, 0],\n",
    "             [2, 3, 0, 1, 0],\n",
    "             [3, 1, 0, 0, 1]])\n",
    "c = np.array([-3, -2, 0, 0, 0])\n",
    "b = np.array([18, 42, 24])\n",
    "\n",
    "Simplex(A, b, c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([-1, -2,  0]),\n",
       " [0, 1, 5],\n",
       " array([0, 0, 1]),\n",
       " [3, 4, 2],\n",
       " array([0.66666667, 3.33333333, 0.33333333]),\n",
       " array([1.33333333, 0.33333333, 1.66666667]))"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# another example test\n",
    "A = np.array([[1, 1, 1, 1, 0, 0],\n",
    "            [-1, 2, -2, 0, 1, 0],\n",
    "            [2, 1, 0, 0, 0, 1]])\n",
    "\n",
    "b = np.array([4, 6, 5])\n",
    "c = np.array([-1, -2, 1, 0, 0, 0])\n",
    "\n",
    "Simplex(A, b, c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As seen above, the function ``Simplex`` outputs the correct values.  ``Simplex`` returns more information than necessary (it does not just return the solution), but it can be useful to see the final values of all the key matrices it uses in the algorithm, so we may gain an intuition into what is going on."
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}