How to: Linear Programming in Python

Optimizing a real world problem

Posted by Daniel Newman on October 11, 2016

Introduction

In this post, I'd like to show you how to setup and solve a linear program using Python. In my example, I'll use Python and the PuLP library to model and solve a linear program1 for assigning students to classes at a science fair. I chose this example since it is easy to understand but difficult to do by hand. I also like using real-world examples and this is an actual problem that I helped an elementary school solve.

Who This is Intended For

This post is primarily tailored towards people who already know the basics of linear programming and want to learn how to model and solve them using Python (I'll also assume you know the basics of Python).

For those who don't have any experience with linear programming, you may find parts of this post hard to understand, particularly the equations, but I hope that if you do read this post, you'll come away with a better understanding of how linear programming can be used to solve an optimization problem.

Many people first learned how to setup and solve Linear Programs with Excel (using the Solver add-on). If this is the case for you, the Excel file in the Github repo for this project has the model set up and might help you follow along here2.

The Problem (in plain words)

A friend who works at an elementary school came to me for help on an interesting problem. She was setting up a science fair for students in her grade and she wanted help placing students into activities. The rules for placements were simple:

  • The fair had 3 periods/sessions, one after another
  • Each student had to do one activity during each period
  • There were 11 activities (each one offered three times in a row)
  • Students couldn't do the same activity more than once
  • There couldn't be too many (>9) or too few kids (<5) in any given activity at once.

Students got to rank their top 4 choices of activities. Our goal was to honor these choices as much as possible so that students got to do activities that they found most interesting, while still fulfilling all of the constraints above.

What You'll Need

To get started, you will need to install:

  • Python (I use 2.7 here)
  • A solver such as GLPK
  • The PuLP library

Install a Solver

A solver is separate program that you install that your Python code is going to call into. You can think of the solver as the engine for solving linear programs. Many solvers are fairly complex programs made with the help of academics knowledgeable in the field of Operations Research. Because these programs need to run very fast/efficiently, they are usually written in a compiled language such as C.

I recommend starting with GLPK since it is free to use and open-source (common alternatives include CPLEX and Gurobi). On a linux machine (debian based) GLPK can be installed with:


$ sudo apt-get install glpk

On a mac, the easiest way is to use Homebrew:


$ brew install glpk

(sorry Windows users... I don't use Windows much, so you are on your own here. BSD users and non-Debian Linux users: I'm not worried about you figuring this out on your own :)

Install PuLP

Solvers such as GLPK can't be called directly from Python (at least not easily), but there are a few libraries that make this easy to do. Personally, I prefer to use PuLP.

PuLP translates your Python code into a mathematical model that can be read into the solver. See here for installation instructions, though if you use pip, you can just run:


$ pip install pulp

The Data

Each student ranked their first four choices of activities. The data can be seen here (names have been changed for the students' privacy). The input file is a CSV with the students' first and last names as well as the top 4 (in ranked order) activities that they wanted to do during the science fair.

The Code

Here's the code for setting up this problem using Python and PuLP. I'll go over the constraints and the objective function bits in more detail in the next section. The code below will:

  • Load in the data (the CSV file containing student choices)
  • Setup the linear model
  • Run the solver
  • Write out the activity assignments to a CSV file for the teacher's use


# Copyright (c) 2016 Daniel Newman (see License file)

import csv
import pulp # library used for linear programming

# PREPARATION: SETUP FUNCTIONS FOR READING/WRITING CSV
#------------------------------------------------------

def read_in_csv_to_list(filename, headers=False):
    """
    :param filename: A string indicating the file name (or relative path)
    :param headers: Whether the first row of the CSV contains headers (default False)
    :return: A list where each item in the list corresponds to a row in the
             CSV (each item in turn is a list corresponding to columns)
    """
    start_row = 0
    if headers:
        start_row = 1
    with open(filename, 'r') as f:
        data = [row for row in csv.reader(f.read().splitlines())]
        return data[start_row:]

def write_list_to_csv(filename, data, headers=None):
    """
    :param filename: A string indicating the file name (or relative path)
    :param headers: A list of headers to write to the first row of the file
    :param data: A list of lists to write to the CSV file
    :return: A list where each item in the list corresponds to a row in the
             CSV (each item in turn is a list corresponding to columns)
    """
    with open(filename, 'wb') as csvfile:
        csv_writer = csv.writer(csvfile, delimiter=',', quotechar='"', quoting=csv.QUOTE_ALL)
        if headers:
            csv_writer.writerow(headers)
        for row in data:
            csv_writer.writerow(row)

# STEP 1: LOAD IN THE RAW DATA
#------------------------------------------------------
student_choice_data = read_in_csv_to_list('student_choices.csv', headers=True)

# STEP 2: SET CONSTANT VARIABLES
#------------------------------------------------------
ACTIVITIES_LIST = ['100ChartPicture', 'BubbleMania', 'Closeto20',
                   'Measuring', 'MysteryNumber', 'PicoPhonyZilch',
                   'PuzzleinaBag', 'RacetoaFlat', 'Salute',
                   'ShapeCodes', 'SteppingStones']

NUM_ACTIVITIES = len(ACTIVITIES_LIST)
NUM_PERIODS = 3 # there are three sessions
NUM_STUDENTS = len(student_choice_data) # this can change based on the input data, but remains constant once here
MIN_NUM_IN_ACTIVITY = 5
MAX_NUM_IN_ACTIVITY = 9

# STEP 3: LOAD STUDENT PREFERENCES INTO A USEFUL FORMAT
#------------------------------------------------------

# assign each activity a number as an activity_id
activities_dict = {activity: n for n, activity in enumerate(ACTIVITIES_LIST)}
activities_dict_reverse = {n: activity for n, activity in enumerate(ACTIVITIES_LIST)}

# assign each student a number as an ID and setup a dictionary with the number as key
# and a dictionary with choices (ordered list of activity_ids by preference) and student name
student_info_dicts = {}
for i, row in enumerate(student_choice_data):
    first_name, last_name, choice1, choice2, choice3, choice4 = row
    name =  last_name + ', ' + first_name
    choices = [choice1, choice2, choice3, choice4]
    student_info_dicts[i] = {'choices': [activities_dict[c] for c in choices],
                             'choices_full_names': choices,
                             'name': name,
                             'assignments': []}


# STEP 4: DEFINE THE PROBLEM TYPE
#------------------------------------------------------

# define the problem as an optimization to maximize the objective function
prob = pulp.LpProblem("StudentClass", pulp.LpMaximize)

# STEP 5: SETUP DECISION VARIABLES
#------------------------------------------------------

# start off by making a simple list of coordinates for decision variables, in the format (i, j)
# where "i" indicates the  student number and "j" indicates the activity number:
decision_var_matrix = []
for i in range(NUM_STUDENTS):
    for j in range(NUM_ACTIVITIES):
        decision_var_matrix.append((i, j))
# at this point, decision_var_matrix = [(0, 0), (0, 1) ... (0, 11), (1, 0), (1, 1), etc...]

# now setup a set of decision variables, one for each period
decision_vars_list = []
for i in range(NUM_PERIODS):
    variable_type_name = 'period_%s_decision_variable' % (i + 1)
    decision_vars_list.append(pulp.LpVariable.dicts(name=variable_type_name,
                                                    indexs=decision_var_matrix,
                                                    lowBound=0,
                                                    upBound=1,
                                                    cat=pulp.LpInteger))
    # note: For continuous variables, set "cat" (i.e. category) to pulp.LpContinuous instead

# We now have one set of decision variables for each period. Think of each set of decision
# variables as an m-by-n matrix of binary (1 or 0) values where 'm' is the number of students and
# 'n' is the number of activities. For the given period, a value of 1 indicates that the student
# is in that activity in that period and a value of 0 means they are not.

# STEP 6: SETUP CONSTRAINTS
#------------------------------------------------------

# CONSTRAINT 1: each student must be in one and only one activity per period
for i in range(NUM_STUDENTS):
    for k in range(NUM_PERIODS):
        vars_to_sum = [decision_vars_list[k][(i, j)] for j in range(NUM_ACTIVITIES)]
        prob += pulp.lpSum(vars_to_sum) == 1

# CONSTRAINT 2: a student cannot repeat an activity
for i in range(NUM_STUDENTS):
    for j in range(NUM_ACTIVITIES):
        vars_to_sum = [decision_vars_list[k][(i, j)] for k in range(NUM_PERIODS)]
        prob += pulp.lpSum(vars_to_sum) <= 1

# CONSTRAINT 3: each activity must have a minimum number of students
for j in range(NUM_ACTIVITIES):
    for k in range(NUM_PERIODS):
        vars_to_sum = [decision_vars_list[k][(i, j)] for i in range(NUM_STUDENTS)]
        prob += pulp.lpSum(vars_to_sum) >= MIN_NUM_IN_ACTIVITY

# CONSTRAINT 4: each activity can only have up to a maximum number of students
for j in range(NUM_ACTIVITIES):
    for k in range(NUM_PERIODS):
        vars_to_sum = [decision_vars_list[k][(i, j)] for i in range(NUM_STUDENTS)]
        prob += pulp.lpSum(vars_to_sum) <= MAX_NUM_IN_ACTIVITY

# CONSTRAINT 5: each student must get 3 of their 4 choices
for i in range(NUM_STUDENTS):
    # remember, student_info_dicts[i]['choices'] looks like [x1, x2, x3, x4] where
    # x1...x4 are the activity numbers for the student's first through fourth
    # choices respectively
    vars_to_sum = [decision_vars_list[k][(i, j)] for j in student_info_dicts[i]['choices']
                                                 for k in range(NUM_PERIODS)]
    prob += pulp.lpSum(vars_to_sum) == 3

# STEP 6: SETUP OBJECTIVE FUNCTION
#------------------------------------------------------

# The objective function is gonna give 1000 "points" for each 1st choice match made,
# 100 for each second choice, 10 for each third choice and 1 for each fourth choice.
# To do this, for each set of decision variables, we iterate through each student
# and get their choices. For each choice, the corresponding decision variables are
# multiplied by the appropriate number of points (1000, 100, 10 or 1) and added to
# the objective_function_parts list. When the list is done being assembled, we set
# the objective function by running lpSum on the list. Note that with this objective
# function, there is no preference on _when_ a student does a preferred activity.
# For example, the same number of points are rewarded whether the student does their
# first choice activity in the first period or the last period.
objective_function_parts = []
for i in range(NUM_STUDENTS):
    for k in range(NUM_PERIODS):
        choice1, choice2, choice3, choice4 = student_info_dicts[i]['choices']
        objective_function_parts.append([decision_vars_list[k][(i, choice1)]*1000])
        objective_function_parts.append([decision_vars_list[k][(i, choice2)]*100])
        objective_function_parts.append([decision_vars_list[k][(i, choice3)]*10])
        objective_function_parts.append([decision_vars_list[k][(i, choice4)]*1])

prob += pulp.lpSum(objective_function_parts)

# STEP 7: SOLVE THE LP
#------------------------------------------------------
solution_found = prob.solve()

# check that a solution was found (anything other than solution_found == 1
# indicates some issue with finding a solution that fits constraints)
assert solution_found == 1, "solution not found"
print("Objective value:", pulp.value(prob.objective))

# STEP 8: EXTRACT THE SOLUTION INTO A USEFUL FORMAT
#------------------------------------------------------

# In this case, let's try to take the solution and put it into a CSV output file. While we're at
# it, we'll collect some data to see how we did

for period in range(NUM_PERIODS):
    for i in range(NUM_STUDENTS):
        for j in range(NUM_ACTIVITIES):
            if decision_vars_list[period][(i,j)].varValue == 1:
                student_info_dicts[i]['assignments'].append(activities_dict_reverse[j])

results = []
for i in range(NUM_STUDENTS):
    student_dict = student_info_dicts[i]
    name = student_dict['name']
    choices = student_dict['choices_full_names']
    choice1, choice2, choice3, choice4 = choices
    assignments = student_dict['assignments']
    assignment1, assignment2, assignment3 = assignments
    row = [name, choice1, choice2, choice3, choice4, assignment1, assignment2, assignment3]
    results.append(row)

headers = ['name', 'choice1', 'choice2', 'choice3', 'choice4',
           'assignment1', 'assignment2' ,'assignment3', 'assignment4']
write_list_to_csv('assignment_outputs.csv', data=results, headers=headers)


# STEP 9: SANITY CHECK OUR RESULTS
#------------------------------------------------------

total_got_choice1 = 0
total_got_choice2 = 0
total_got_choice3 = 0
total_got_choice4 = 0
# due to constraint 5, we know that every student gets 3 of 4 choices, but it doesn't
# hurt to check that over here
total_got_3_of_4 = 0

for i in range(NUM_STUDENTS):
    choices = student_info_dicts[i]['choices_full_names']
    choice1, choice2, choice3, choice4 = choices
    assignments = student_info_dicts[i]['assignments']

    got_choice1, got_choice2, got_choice3, got_choice4 = 0, 0, 0, 0
    if choice1 in assignments: got_choice1 = 1
    if choice2 in assignments: got_choice2 = 1
    if choice3 in assignments: got_choice3 = 1
    if choice4 in assignments: got_choice4 = 1
    if (got_choice1 + got_choice2 + got_choice3 + got_choice4) == 3:
        total_got_3_of_4 += 1
    total_got_choice1 += got_choice1
    total_got_choice2 += got_choice2
    total_got_choice3 += got_choice3
    total_got_choice4 += got_choice4

print("percent who got choice 1:", float(total_got_choice1) / NUM_STUDENTS * 100)
print("percent who got choice 2:", float(total_got_choice2) / NUM_STUDENTS * 100)
print("percent who got choice 3:", float(total_got_choice3) / NUM_STUDENTS * 100)
print("percent who got choice 4:", float(total_got_choice4) / NUM_STUDENTS * 100)
print("percent who got 3 of 4 choices:", float(total_got_3_of_4) / NUM_STUDENTS * 100)

The Problem (in mathematical terms)

Let's start by defining some variables:

  • `M` denotes the total number of students (85 in this example)
  • `N` denotes the number of activities (11)
  • `O` denotes the number of periods (3)
  • Let `x_{ijk}` be a 0/1 value where `x_{ijk}=1` indicates that student `i` is assigned to activity `j` during period `k`, where `i=1,...,M`, `j=1,...,N` and `k=1,...,O`
  • Let `C_{1ij}` be a 0/1 value where `C_{1ij}=1` indicates that student `i` ranked activity `j` as their first choice, where `i=1,...,M`, `j=1,...,N`
  • Let `C_{2ij}`, `C_{3ij}` and `C_{4ij}` be the same thing, but for 2nd through 4th ranked choices

The objective function (i.e. the equation that we are trying to maximize in order to get the best student-to-activity assignments) is then defined as:

Maximize `\sum_{i=1}^{M}\sum_{j=1}^{N}\sum_{k=1}^{O} 1000x_{ijk}C_{1ij}` `+ \sum_{i=1}^{M}\sum_{j=1}^{N}\sum_{k=1}^{O} 100x_{ijk}C_{2ij}` `+ \sum_{i=1}^{M}\sum_{j=1}^{N}\sum_{k=1}^{O} 10x_{ijk}C_{3ij}` `+ \sum_{i=1}^{M}\sum_{j=1}^{N}\sum_{k=1}^{O} x_{ijk}C_{4ij}`

In other words, we'll give 1000 points for each student who gets their first choice, 100 for each second, 10 for each third and 1 for each fourth choice so that fulfilling higher ranked choices is preferred (the goal is to maximize the total number of "points").

Here's the relevant code snippet:


objective_function_parts = []
for i in range(NUM_STUDENTS):
    for k in range(NUM_PERIODS):
        choice1, choice2, choice3, choice4 = student_info_dicts[i]['choices']
        objective_function_parts.append([decision_vars_list[k][(i, choice1)]*1000])
        objective_function_parts.append([decision_vars_list[k][(i, choice2)]*100])
        objective_function_parts.append([decision_vars_list[k][(i, choice3)]*10])
        objective_function_parts.append([decision_vars_list[k][(i, choice4)]*1])

prob += pulp.lpSum(objective_function_parts)

Below are the mathematical representations of the constraints as well as the relevant code snippets for each. You can see that the mathematical notation translates fairly well to the code implementations:

1) Each student is in one and only one activity per period:

`\sum_{j=1}^{N} x_{ijk} = 1`, for `i=1,...,M` and `k=1,...,O`


for i in range(NUM_STUDENTS):
    for k in range(NUM_PERIODS):
        vars_to_sum = [decision_vars_list[k][(i, j)] for j in range(NUM_ACTIVITIES)]
        prob += pulp.lpSum(vars_to_sum) == 1

2) A student cannot repeat an activity:

`\sum_{k=1}^{P} x_{ijk} <= 1`, for `i=1,...,M` and `j=1,...,N`


for i in range(NUM_STUDENTS):
    for j in range(NUM_ACTIVITIES):
        vars_to_sum = [decision_vars_list[k][(i, j)] for k in range(NUM_PERIODS)]
        prob += pulp.lpSum(vars_to_sum) <= 1

3) Each activity must have a minimum number of students per period (let `\alpha` represent some minimum number which in our example is 5):

`\sum_{i=1}^{M} x_{ijk} >= \alpha`, for `j=1,...,N` and `k=1,...,O`


for j in range(NUM_ACTIVITIES):
    for k in range(NUM_PERIODS):
        vars_to_sum = [decision_vars_list[k][(i, j)] for i in range(NUM_STUDENTS)]
        prob += pulp.lpSum(vars_to_sum) >= MIN_NUM_IN_ACTIVITY

4) Each activity can have up to a maximum number of students (let `\beta` represent some maximum number which in our example is 9):

`\sum_{i=1}^{M} x_{ijk} <= \beta`, for `j=1,...,N` and `k=1,...,O`


for j in range(NUM_ACTIVITIES):
    for k in range(NUM_PERIODS):
        vars_to_sum = [decision_vars_list[k][(i, j)] for i in range(NUM_STUDENTS)]
        prob += pulp.lpSum(vars_to_sum) <= MAX_NUM_IN_ACTIVITY

5) Each student must get three out of their four choices3:

`\sum_{j=1}^{N}\sum_{k=1}^{O} x_{ijk} = 3`, for `i=1,...,M`


for i in range(NUM_STUDENTS):
    # remember, student_info_dicts[i]['choices'] looks like [x1, x2, x3, x4] where
    # x1...x4 are the activity numbers for the student's first through fourth
    # choices respectively
    vars_to_sum = [decision_vars_list[k][(i, j)] for j in student_info_dicts[i]['choices']
                                                 for k in range(NUM_PERIODS)]
    prob += pulp.lpSum(vars_to_sum) == 3

Conclusions

The example I gave here shows how I solved a particular problem using Python and PuLP. I hope that this example gives you some idea as to how to setup your own linear, integer and mixed-integer programs using these tools. All the code here is available on the github repo that I setup for this tutorial. My code uses the MIT License, so feel free to adapt it for your own use, whether for personal, academic or commercial use.

If you have any questions or feedback, please reach out to me here.

Footnotes

[1] The example here is technically an Integer Linear Program, but the process for setting up a "regular" linear program is very similar and I'll point out how to do this in the code example.

[2] In the Excel file, I setup the entire model, which can be viewed by going to "Tools>solver...". However, because there are too many decision variables, the standard solver that comes with Excel won't run when you click "solve". To do so, you'll need to install OpenSolver. For large problems, you will find that a solver such as GLPK is much faster and gets better solutions than the solver engine that Excel uses.

[3] With the given data, it won't be possible for everyone to get all 3 of their top choices (though it will be possible for everyone to get their first two choices). This constraint ensures that no one gets choices that they didn't rank in their top 4. In other words, all else equal, it's better for both Sally and John to get their 1st, 2nd and 4th choices than for Sally to get her top 3 choices and John to get his 1st, 2nd choices plus an activity that he didn't rank in his top 4.

Back to the top