Skip to content

Instantly share code, notes, and snippets.

@tuf22191

tuf22191/Linear Programming Solver 5 10 2018

Created May 17, 2018 23:31
Show Gist options
  • Save tuf22191/65e3a26deb91ca95daa96baa4c80893f to your computer and use it in GitHub Desktop.
Save tuf22191/65e3a26deb91ca95daa96baa4c80893f to your computer and use it in GitHub Desktop.
Download ZIP
Linear Programming Simplex Method!
Raw
Linear Programming Solver 5 10 2018
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this tutorial we are going to build our very own practical linear programming solver using the Simplex Method. We are going to do the standard MAXIMIZATION (not minimization) problem discussed in the video above (here is a link just in case: https://www.youtube.com/watch?v=XK26I9eoSl8&t=300s) . Another great article to understand linear programming some more is here https://www.analyticsvidhya.com/blog/2017/02/lintroductory-guide-on-linear-programming-explained-in-simple-english/. as well."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"#let's import our libraries! \n",
"import numpy as np\n",
"import pandas as pd"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our input vectors look like this from the video: \n",
"\n",
"$$ \\left[\n",
" \\begin{matrix}\n",
" P = 6x+5y+4z \\\\\n",
" 2x+z+y\\leq180 \\\\\n",
" 3y+2z+x\\leq300 \\\\\n",
" 2z+y+2x\\leq240 \\\\\n",
" \\end{matrix}\n",
"\\right] $$\n",
"\n",
"To convert it we have to alphabetically put the independent variables on the left hand side and the result values on the right hand side: \n",
"\n",
"(matrix A below)\n",
"$$ \\left[\n",
"\\begin{array}{ccc|c}\n",
" 6x &5y &4z & P \\\\\n",
" 2x &y &z & \\leq180 \\\\\n",
" x & 3y & 2z & \\leq300 \\\\\n",
" 2x & y &2z & \\leq240 \\\\\n",
"\\end{array}\n",
"\\right] $$\n",
"\n",
"We need to get rid of the less than or equal to signs by introducing the slack variables. Also the objective function we are trying to maximize is the top equation and we want to maximize the value of P. According to the algorithm, we can only have constants on the right hand side, as in the P variable is actually considered a slack variable. We simple subtract it to the other side:\n",
"(matrix B below)\n",
"$$ \\left[\n",
"\\begin{array}{cccc|c}\n",
" -6x &-5y &-4z & P & 0 \\\\\n",
" 2x &y &z & s1 & =180 \\\\\n",
" x & 3y & 2z & s2 & =300 \\\\\n",
" 2x & y &2z & s3 &=240 \\\\\n",
"\\end{array}\n",
"\\right] $$\n",
"\n",
"Let's generate this in cleaner format with the columns added for each slack variable: "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"df = pd.DataFrame({1:['-6x','2x', 'x', '2x'], 2:['-5y','y', '3y', 'y'],\n",
" 3:['-4z','z', '2z', '2z'], 4:['P','s1', 's2', 's3'],\n",
" 5:['0','180', '300', '240']})\n",
"df # show the dataframe"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's clean up the dataframe a bit. The 4th column is misleading. Each slack variable should have it's own column: \n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style>\n",
" .dataframe thead tr:only-child th {\n",
" text-align: right;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: left;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>constants</th>\n",
" <th>p</th>\n",
" <th>s1</th>\n",
" <th>s2</th>\n",
" <th>s3</th>\n",
" <th>x</th>\n",
" <th>y</th>\n",
" <th>z</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0</td>\n",
" <td>P</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>-6x</td>\n",
" <td>-5y</td>\n",
" <td>-4z</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>180</td>\n",
" <td>0</td>\n",
" <td>s1</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>2x</td>\n",
" <td>1y</td>\n",
" <td>z</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>300</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>s2</td>\n",
" <td>0</td>\n",
" <td>1x</td>\n",
" <td>3y</td>\n",
" <td>2z</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>240</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>s3</td>\n",
" <td>2x</td>\n",
" <td>1y</td>\n",
" <td>2z</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" constants p s1 s2 s3 x y z\n",
"0 0 P 0 0 0 -6x -5y -4z\n",
"1 180 0 s1 0 0 2x 1y z\n",
"2 300 0 0 s2 0 1x 3y 2z\n",
"3 240 0 0 0 s3 2x 1y 2z"
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df2 = pd.DataFrame({'x':['-6x','2x', '1x', '2x'], 'y':['-5y','1y', '3y', '1y'],\n",
" 'z':['-4z','z', '2z', '2z'], 'p':['P', '0', '0', '0'],\n",
" 's1':['0','s1', '0', '0'], 's2':['0','0', 's2', '0'], \n",
" 's3':['0','0', '0', 's3'], \n",
" 'constants':['0','180', '300', '240']})\n",
"df2 # show the dataframe"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Firstly, we have to put the objective function (the function we are maximizing) to the bottom. "
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"swapping rows in place\n",
" constants p s1 s2 s3 x y z\n",
"1 180 0 s1 0 0 2x 1y z\n",
"2 300 0 0 s2 0 1x 3y 2z\n",
"3 240 0 0 0 s3 2x 1y 2z\n",
"0 0 P 0 0 0 -6x -5y -4z\n",
"after resetting the index from 0 to 3\n",
" constants p s1 s2 s3 x y z\n",
"0 180 0 s1 0 0 2x 1y z\n",
"1 300 0 0 s2 0 1x 3y 2z\n",
"2 240 0 0 0 s3 2x 1y 2z\n",
"3 0 P 0 0 0 -6x -5y -4z\n"
]
}
],
"source": [
"# to swap rows in place in dataframes https://stackoverflow.com/questions/32547440/python-pandas-how-to-move-one-row-to-the-first-row-of-a-dataframe\n",
"indexes_of_rows_at_top = list(range(1,len(df2)))\n",
"indexes_of_top_row = [0]\n",
"df3 = df2.reindex(indexes_of_rows_at_top+indexes_of_top_row)\n",
"print('swapping rows in place')\n",
"print(df3)\n",
"df3.index = range(len(df3))\n",
"print('after resetting the index from 0 to 3')\n",
"print(df3)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now , let's keep only the coefficients in the matrix and we can start the Simplex Method!"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style>\n",
" .dataframe thead tr:only-child th {\n",
" text-align: right;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: left;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>x</th>\n",
" <th>y</th>\n",
" <th>z</th>\n",
" <th>s1</th>\n",
" <th>s2</th>\n",
" <th>s3</th>\n",
" <th>p</th>\n",
" <th>constants</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>2</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>180</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1</td>\n",
" <td>3</td>\n",
" <td>2</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>300</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>2</td>\n",
" <td>1</td>\n",
" <td>2</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>240</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>-6</td>\n",
" <td>-5</td>\n",
" <td>-4</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" x y z s1 s2 s3 p constants\n",
"0 2 1 1 1 0 0 0 180\n",
"1 1 3 2 0 1 0 0 300\n",
"2 2 1 2 0 0 1 0 240\n",
"3 -6 -5 -4 0 0 0 1 0"
]
},
"execution_count": 60,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df4 = pd.DataFrame({'x':[2, 1, 2,-6], 'y':[1, 3, 1, -5],\n",
" 'z':[1, 2, 2, -4], 'p':[0, 0, 0, 1],\n",
" 's1':[1, 0, 0,0], 's2':[0, 1, 0,0], \n",
" 's3':[0, 0, 1,0], \n",
" 'constants':[180, 300, 240,0]})\n",
"columns_order = ['x', 'y','z', 's1', 's2', 's3', 'p', 'constants']\n",
"df5 = df4[columns_order] # we need to order the variables from independent -> to slack -> to constants\n",
"df5 # show the dataframe"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"Now we are ready to manipulate the above matrix using the simplex method!!! \n",
"Recall that we have to find the smallest value . \n",
"The key thing to realize is we can convert the dataframe to a numpy array. First let's convert each column in the dataframe to np.float64 type so we can do accurate calculations and then we will convert to a numpy 2d array so we can do matrix manipulations"
]
},
{
"cell_type": "code",
"execution_count": 61,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x float64\n",
"y float64\n",
"z float64\n",
"s1 float64\n",
"s2 float64\n",
"s3 float64\n",
"p float64\n",
"constants float64\n",
"dtype: object\n",
"[[ 2. 1. 1. 1. 0. 0. 0. 180.]\n",
" [ 1. 3. 2. 0. 1. 0. 0. 300.]\n",
" [ 2. 1. 2. 0. 0. 1. 0. 240.]\n",
" [ -6. -5. -4. 0. 0. 0. 1. 0.]]\n",
"(4L, 8L)\n"
]
}
],
"source": [
"#step 1. transform to np.float64 type\n",
"for column_name in df5.columns.tolist():\n",
" df5[column_name] = df5[column_name].astype(np.float64)\n",
"\n",
"print(df5.dtypes) # shows the data types of each of the columns \n",
"#step 2. get the numpy array\n",
"simplex_matrix = df5.values\n",
"print(simplex_matrix)\n",
"print(simplex_matrix.shape)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Great that is exactly what we wanted! Let's follow the algorithm (keep df5 around: we will still need df5, our dataframe to see which column is which!):"
]
},
{
"cell_type": "code",
"execution_count": 80,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"#helper function to find index of number in array\n",
"def indexOf(number, nparr):\n",
" for x in range(0, len(nparr)):\n",
" if(nparr[x] == number):\n",
" return x\n",
"\n",
"#helper function to find the pivot row given you found the pivot column already\n",
"def findPivotRow(mat, pivot_col):\n",
" mat_copy = np.array(mat) #does a \"deep copy\"\n",
" constant_column = mat[:, -1] # gets every row with column index = -1 (which corresponds to last column)\n",
" pivot_column = mat[:, pivot_col]\n",
" #the division step to find pivot_row\n",
" division_step = np.array(pivot_column) # does a \"deep copy\", we will\n",
" for x in range(0, len(pivot_column)):\n",
" if((pivot_column[x] == 0)):# we will be dividing by zero and hence put an np.inf instead\n",
" division_step[x] = 100000000000.0\n",
" elif(pivot_column[x]<0): #pivot must be a positive number \n",
" division_step[x] = 100000000000.0 # set to infinity as well\n",
" else: #we are not dividing by zero or by negative number -> just do simple division\n",
" division_step[x] = constant_column[x]/pivot_column[x]\n",
" \n",
" # now just find the minimum division step\n",
" minval = min(division_step) #this is why np.infinities help us to ignore the minimum values\n",
" row_col = indexOf(minval, division_step)\n",
" return row_col\n",
"#step 3 simplex algorihtm\n",
"def simplex_method(mat): \n",
" \"\"\"\n",
" Assumes that mat is in standard form of: independent variables-> to slack -> to constants column at the end\n",
" \"\"\"\n",
" #mat = mat*1.0\n",
" a = mat.astype(np.float64)\n",
" mat = a\n",
" print(\"matrix is \\n\\n\")\n",
" print(mat)\n",
" counter=0\n",
" while(min(mat[-1])<0.0): # while we have negative in bottom row\n",
" if(counter == 500):\n",
" print(\"problem\")\n",
" break\n",
" counter+=1 \n",
" #find the pivot column \n",
" a = min(mat[-1])\n",
" pivot_col =indexOf(a, mat[-1])\n",
" #print(\"pivot_col is \"+str(pivot_col))\n",
" #find the pivot row\n",
" pivot_row = findPivotRow(mat, pivot_col)\n",
" #print(\"pivot_row is \"+str(pivot_row))\n",
" #set pivot to one \n",
" #print(\"bef first row op, mat is \"+str(mat.shape)+\"\\n\\n\")\n",
" #print(mat)\n",
" mat[pivot_row] = mat[pivot_row]/mat[pivot_row, pivot_col]\n",
" #print(\"bef rest of row ops, mat is \"+str(mat.shape)+\"\\n\\n\")\n",
" #print(mat)\n",
" #set other values to zero\n",
" other_row_indices = list(range(0, len(mat)))\n",
" other_row_indices.remove(pivot_row) #we have already updated the pivot row\n",
" for x in other_row_indices:\n",
" #print(mat[pivot_row, pivot_col])\n",
" divider = mat[x, pivot_col]/mat[pivot_row, pivot_col]\n",
" mat[x] = mat[x] - divider*mat[pivot_row]\n",
" #print(mat[x])\n",
"# print(\"after row ops, mat is \"+str(mat.shape)+\"\\n\\n\")\n",
"# print(mat)\n",
" #done with loop, let's identify the values of each variable \n",
" # we can store the value of the variables according to the index of the variables\n",
" \n",
" #determine the basic and non basic columns! \n",
" basic_cols = np.array(mat[0])\n",
" #set basic_cols[i] to 1 if i'th column is basic.\n",
" for i in range(0, len(mat[0])):\n",
" basic_cols[i] = 1 if sum(mat[:, i]) == 1 else 0\n",
" \n",
" #determine the values of each of the columns\n",
" column_values_left_to_right = []\n",
" for i in range(0, len(basic_cols)):\n",
" if(basic_cols[i]==1): #is basic column\n",
" row_val_where_1 = indexOf(1, mat[:, i]) #find where we have a 1\n",
" column_values_left_to_right.append(mat[row_val_where_1,-1])\n",
" else: #is not basic column\n",
" column_values_left_to_right.append(0) #according to video\n",
" \n",
" return column_values_left_to_right\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's call our simplex method! We will also map the output to the column names since we just get a list of numbers back. "
]
},
{
"cell_type": "code",
"execution_count": 78,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"matrix is \n",
"\n",
"\n",
"[[ 2. 1. 1. 1. 0. 0. 0. 180.]\n",
" [ 1. 3. 2. 0. 1. 0. 0. 300.]\n",
" [ 2. 1. 2. 0. 0. 1. 0. 240.]\n",
" [ -6. -5. -4. 0. 0. 0. 1. 0.]]\n",
"\n",
"\n",
"\n",
"##########################################\n",
"##########################################\n",
"##########################################\n",
"\n",
"x is 48.0\n",
"y is 84.0\n",
"z is 0\n",
"s1 is 0\n",
"s2 is [[ 0.00000000e+00 0.00000000e+00 2.00000000e-01 2.60000000e+00\n",
" 8.00000000e-01 0.00000000e+00 1.00000000e+00 7.08000000e+02]]\n",
"s3 is 60.0\n",
"p is 708.0\n",
"constants is 0\n"
]
}
],
"source": [
"column_values = simplex_method(simplex_matrix)\n",
"#use df5 to get the column names:\n",
"print('\\n'*3+\"##########################################\\n\"*3)\n",
"for index, column in enumerate(df5.columns.tolist()):\n",
" print(str(column) + \" is \" + str(column_values[index]))"
]
},
{
"cell_type": "code",
"execution_count": 81,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"matrix is \n",
"\n",
"\n",
"[[ 2. 3. 2. 1. 0. 0. 1000.]\n",
" [ 1. 1. 2. 0. 1. 0. 800.]\n",
" [ -7. -8. -10. 0. 0. 1. 0.]]\n",
"\n",
"\n",
"\n",
"##########################################\n",
"##########################################\n",
"##########################################\n",
"\n",
"x1 is 200.0\n",
"x2 is 0\n",
"x3 is 300.0\n",
"s1 is 0\n",
"s2 is 0\n",
"P is 4400.0\n"
]
}
],
"source": [
"#there is a slight bug at the s2 so please use the code with caution!!\n",
"#let's try another problem found here: https://www.youtube.com/watch?v=WeK4JjNLSgw\n",
"simplex_mat2 = np.array([\n",
" [2.0,3.0,2.0,1.0,0.0,0.0,1000.0],\n",
" [1.0,1.0,2.0,0.0,1.0,0.0,800.0],\n",
" [-7.0,-8.0,-10.0,0.0,0.0,1.0,0.0]\n",
"])\n",
"column_values = simplex_method(simplex_mat2)\n",
"#use df5 to get the column names:\n",
"print('\\n'*3+\"##########################################\\n\"*3)\n",
"for index, column in enumerate(['x1', 'x2','x3','s1','s2','P']):\n",
" print(str(column) + \" is \" + str(column_values[index]))\n"
]
},
{
"cell_type": "code",
"execution_count": 83,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"matrix is \n",
"\n",
"\n",
"[[ 1. 3. 2. 1. 0. 0. 10.]\n",
" [ 1. 5. 1. 0. 1. 0. 8.]\n",
" [ -8. -10. -7. 0. 0. 1. 0.]]\n",
"\n",
"\n",
"\n",
"##########################################\n",
"##########################################\n",
"##########################################\n",
"\n",
"x1 is 8.0\n",
"x2 is 0\n",
"x3 is 0\n",
"s1 is 2.0\n",
"s2 is 0\n",
"P is 64.0\n"
]
}
],
"source": [
"#yet another problem! found here: https://www.youtube.com/watch?v=iwDiG2mR6FM\n",
"simplex_mat3 = np.array([\n",
" [1,3,2,1,0,0,10],\n",
" [1,5,1,0,1,0,8],\n",
" [-8,-10,-7,0,0,1,0]\n",
"])\n",
"column_values = simplex_method(simplex_mat3)\n",
"#use df5 to get the column names:\n",
"print('\\n'*3+\"##########################################\\n\"*3)\n",
"for index, column in enumerate(['x1', 'x2','x3','s1','s2','P']):\n",
" print(str(column) + \" is \" + str(column_values[index]))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"That's all for now folks! Next we write a Quadratic programming solver!"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.14"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment