Introduction to Linear Programming with PuLP

Try me

Requirements

Install in your environment

Pip Installation

The simplest way to install PuLP in your environment is using pip. If you have installed Python and pip in your environment, just open a terminal and try:

pip install pulp

Conda Installation

If you use Conda, open a Conda Terminal and try:

conda install –c conda-forge pulp

Google Colabs installation

Run the following code cell to try this notebook in Google Colabs:

!pip install pulp

Binder installation

Run the following code cell to try this notebook in Binder:

!pip install pulp
!pip install pandas
!pip install numpy

The PuLP library

In this tutorial, we will learn to model and solve Linear Programming Problems using the Python open source Linear Programming library PuLP.

To guide this example, we will use a simple LPP formulated in class:

maximise \(z = 300x + 250y\)

Subject to:

\(2x + y \leq 40\)

\(x + 3y \leq 45\)

\(x \leq 12\)

In Pulp this problem can be solved with the following code cell:

# Let´s start importing the library PuLP to solve linear programs
import pulp

# Create an instance of the problem class using LpProblem
model = pulp.LpProblem("Production_Mix_example", pulp.LpMaximize) #this will create an instance of an LP Maximise problem

# Create the variables of the problem
x = pulp.LpVariable('x', lowBound=0, cat='Continuous')
y = pulp.LpVariable('y', lowBound=0, cat='Continuous')

# Add the objective function to the model
model += 300 * x + 250 * y, "Profit"

# And the constraints
model += 2 * x + y <= 40, "Man Power"
model += x + 3 * y <= 45, "Machine Operating Time"
model += x <=12, "Marketing"

# solve the problem
model.solve()
pulp.LpStatus[model.status]

# Print our decision variable values
print("Production of Product A = {}".format(x.varValue))
print("Production of Product B = {}".format(y.varValue))

Ok, let us see step by step the main functions used to model and solve this problem with PuLP.

Code Explanation

Problem Class LpProblem

PuLP uses classes providing different methods to model and solve LPPs. The class that will contain our model is the LpProblem class. To create a new LpProblem we use the pulp LpProblem function:

• LpProblem(name=’None’, sense=1): Creates a new Linear Programming Problem. The parameter name (default ‘None’) assigns a name to the problem. The parameter sense (either pulp.LpMinimise or pulp.LpMaximize) sets the type of objective function. The default is minimise.

So, the line code:

model = pulp.LpProblem("Production_Mix_example", pulp.LpMaximize)

Created a problem called “Production_Mix_example” and set the objective function to maximise.

Variable class LpVariable

The definition of a LPP program with PuLP is very similar to the standard procedure used to model a problem. First, we need to define the unknown variables in our problem. For this purpose we use the class LpVariable. The function LpVariable allows us to create a variable:

• LpVariable(name, lowBound=None, upBound=None, cat=’Continuous’, e=None): Creates an instance of variable with the following properties:

  • Name: The name of the variable to be used in the solution.

  • lowBoud: The lower bound of the variable, the default is unsrestricted (-Inf).

  • upBound: The upper bound of the variable. The default is unrestricted (Inf).

  • cat: Either ‘Continuous’ for continuous variables, ‘Binary’ for binary variables or ‘Integer’ for Integer variables. We will see in detail binary and integer variables in the course unit for Mixed Integer Programming, but now you know that you will be able to model and solve this type of problems with PuLP. The default is ‘Continuous’.

  • e: This parameter is outside the scope of this course and can be neglected for now.

We can define the variables of our problem using the LpVariable function:

# Create the variables of the problem
x = pulp.LpVariable('x', lowBound=0, cat='Continuous')
y = pulp.LpVariable('y', lowBound=0, cat='Continuous')

Note that we have created two variables, x and y, with lower bound 0 and type continuous. We could have also use upBound to set the upper bound of the variables. For instance, if we wanted to set the upper bound of x to 12, we could have used this parameter.

Adding expressions

In PuLP, both objective function and constraints are expressions (algebraic expressions containing variables) that have to be added to the instance problem using the standard operand ‘+=’. For instance, to add the objective function in this example, we could write:

model += 300 * x + 250 * y, "Profit"

With this line of code, we have added a new expression with name “Profit” that multiplies the technological coefficients to the variables x and y (as defined in the code snippet in the previous section).

Note that we have used the operand ‘+=’ to add the expression to the model variable. We can use the same procedure to add the constraints:

# And the constraints
model += 2 * x + y <= 40, "Man Power"
model += x + 3 * y <= 45, "Machine Operating Time"
model += x <=12, "Marketing"

Note that we added names to the constraints. This is not mandatory, but it is a good practice to do so. Note also that we needed one line of code per constraint. This is not a problem for this simple example, but it can be a problem for more complex problems. We will see how to solve this issue in the next section.

Solving the problem

Once we have defined the problem, we can solve it using the solve method of the LpProblem class:

# solve the problem
model.solve()

Note that we have not specified any solver. PuLP will use the default solver (CBC) to solve the problem. We will see in next chapters how to use other solvers.

Checking the status of the solution

Once the problem has been solved, we can check the status of the solution using the LpStatus method of the LpProblem class:

pulp.LpStatus[model.status]

The status of the solution can be:

• Optimal: An optimal solution has been found.

• Not Solved: The problem has not been solved yet.

• Infeasible: There are no feasible solutions (e.g. if you set the constraints x <= 1 and x >=2).

• Unbounded: The constraints are not bounded, maximising the solution will tend towards infinity (e.g. if the only constraint was x >= 3).

• Undefined: The optimal solution may exist but may not have been found.

• Error: An error occurred (e.g. if you try to retrieve the infeasible or unbounded status from a model which is not infeasible or unbounded).

• None: No solution status exists (e.g. if you create a model but do not solve it).

Retrieving the solution

Once the problem has been solved, we can retrieve the solution using the varValue method of the LpVariable class:

print(x.varValue)

Making PuLP programs more scalable

In this section, we will see how to make PuLP programs more scalable. We will use some Python basic concepts related to iterables. You can check these tutorials if you need a refresher:

• Basic Python Iterators

• Intermediate Python Iterators

Using dictionary variables

Note however that if we use the script above, the problem just does not scale up because we need a line of code for every unknown. What if we have hundreds of unknowns? Luckily for us, PuLP provides a convenient method to write more efficient codes for our program, the LpVariable.dicts method, which basically allows us to create a set of variables with the same category, upper bounds and lower bounds at once:

• LpVariable.dicts(name, index, lowBound=None, upBound=None, cat=’Continuous’): Creates a dictionary containing variables of type cat (default ‘Continuous’), indexed with the keys contained in the iterable index and bounded by lowBound (default -Inf) and upBound (default Inf).

For instance, we can define our variables using the LpVariable.dicts method like this:

# First we define a tuple with the variable names x and y
variable_names = (1,2)
# Then we create a variable from a dictionary, using the variable names as keys
X = pulp.LpVariable.dicts("vars", variable_names, lowBound=0, cat='Continuous')

Note that we used a tuple with the variable names to create both variables in a single line of code. If we had 20 variables, we could have used a tuple with the 20 variable names and created all the variables in a single line of code. This is a much more scalable way to create variables in PuLP!

Using lpSum

We can also use the lpSum function to create linear expressions from a vector of variables, just as we can use summations to express the objective function and constraints in a more compact mathematical form. For instance, we can write the objective function like this:

\(C = [c_1, c_2] = [300, 250]\)

\(X = [x_1, x_2]\)

\(\max z = 300x_1 + 250x_2 = \sum_{j=1}{2}{c_j*x_j}\)

PuLP provides a convenience function named lpSum to achieve the same result efficiency. lpSum takes an array of expressions and returns the summation of the elements in the array. Let us see it action:

# We define the objective function coefficients
C = [300, 250]
# Then we add the objective function to the model like
model += (
    pulp.lpSum([
        C[i] * X[variable_names[i]]
        for i in range(len(X))])
), "Profit"

Notice that we have used list comprehension to create the array where every element is the product of a coefficient times the corresponding decision variable, and passed this array to the lpSum function using an index array. We have used range and length to go through all decision variables.

Using matrices and arrays to define problem constraints

We can also use matrices and arrays to define the constraints of the problem. For instance, we can write the constraints of the problem like this:

\(A = \begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \\ a_{31} & a_{32} \end{bmatrix} = \begin{bmatrix} 2 & 1 \\ 1 & 3 \\ 1 & 0 \end{bmatrix}\)

\(b = \begin{bmatrix} b_1 \\ b_2 \\ b_3 \end{bmatrix} = \begin{bmatrix} 40 \\ 45 \\ 12 \end{bmatrix}\)

\(\sum_{j=1}^{2}{a_{ij}*x_j} \leq b_i \forall i\)

In Python, we can achieve the same compactness defining the coefficients of the constraints in of a matrix, the LHS coefficients in a vector, and use a for loop to add the constraints to the model:

# need We also define the name for the constraints
constraint_names = ['Man Power', 'Machine Operating Time', 'Marketing']


# Define the RHS constraint coefficients in a matrix
A=[[2, 1], #Coefficients of the first constraint
   [1, 3], #Coefficients of the second constraint
   [1, 0]] #Coefficients of the third constraint

# And vector b
b = [40, 45, 12] #limits of the three constraints

# need We also define the name for the constraints
constraint_names = ['Man Power', 'Machine Operating Time', 'Marketing']
# Now we add the constraints using
# model += expression, name
# eg model += 2*X + y <= 40
# We add all constraints in a loop, using a vector and the function lpSum to generate the linear expression:
for i in range(len(A)):
    model += pulp.lpSum([
        A[i][j] * X[variable_names[j]]
        for j in range(len(variable_names))]) <= b[i] , constraint_names[i]

Here we just put all data in iterables, and used a for loop to add the constraints to the model. This is a much more scalable way to add constraints to a model!

Let’s put everything together and solve the problem:

import pulp
# First we define a tuple with the variable names x and y
variable_names = (1,2)
# Then we create a variable from a dictionary, using the variable names as keys
X = pulp.LpVariable.dicts("x",
                                     (i for i in variable_names),
                                     lowBound=0,
                                     cat='Continuous')

# We define the objective function coefficients
coefficients = [300, 250]

# We create the model
model = pulp.LpProblem("Profit maximising problem", pulp.LpMaximize)

# We use the function lpSum to generate the linear expression from the coeeficients and variables
model += (
    pulp.lpSum([
        coefficients[i] * X[variable_names[i]]
        for i in range(len(X))])
), "Profit"

# Define the RHS constraint coefficients in a matrix
A=[[2, 1], #Coefficients of the first constraint
   [1, 3], #Coefficients of the second constraint
   [1, 0]] #Coefficients of the third constraint

# And vector b
b = [40, 45, 12] #limits of the three constraints

# need We also define the name for the constraints
constraint_names = ['Man Power', 'Machine Operating Time', 'Marketing']
# Now we add the constraints using
# model += expression, name
# eg model += 2*X + y <= 40
# We add all constraints in a loop, using a vector and the function lpSum to generate the linear expression:
for i in range(len(A)):
    model += pulp.lpSum([
        A[i][j] * X[variable_names[j]]
        for j in range(len(variable_names))]) <= b[i] , constraint_names[i]

# Solve our problem
model.solve()
pulp.LpStatus[model.status]

# Print our decision variable values
for variable in model.variables():
    print ("{} = {}".format(variable.name, variable.varValue))

x_1 = 12.0
x_2 = 11.0

C:\Users\ffraile\PycharmProjects\operations-research-notebooks\venv\lib\site-packages\pulp\pulp.py:1352: UserWarning: Spaces are not permitted in the name. Converted to '_'
  warnings.warn("Spaces are not permitted in the name. Converted to '_'")

Now that we have created our model, we can get the solution just by calling the method solve(). The status of the solution can be read in the LpStatus attribute:

Displaying the solution using Pandas

Pandas is a great library to manipulate tables. You can find a basic tutorial here.

We can use Pandas to display the solution in a nice table. First we need to import Pandas:

import pandas as pd

We are also going to use the library display to show styled output using Markdown notation.

Now, let us display the solution in a nice table using Pandas. We are going to first display the solution value using markdown and then we will use Pandas to create a table with the results.

# We are going to use panda to display the results as tables using Panda
import pandas as pd
#And we will use numpy to perform array operations
import numpy as np
#We will use display and Markdown to format the output of code cells as Markdown
from IPython.display import display, Markdown


# Solution
max_z = pulp.value(model.objective)

#We use display and Mardown to show the value using markdown
display(Markdown("The value of the objective function is **%.2f**"%max_z))


# Print our decision variable values
display(Markdown("The following tables show the values obtained: "))
# First we create a dataframe from the dictionary of the solution. We want to use the variable indexes to present the results and
# place the different values provided by the solver in the data frame.
var_df = pd.DataFrame.from_dict(X, orient="index",
                                columns = ["Variables"])
# First we add the solution. We apply a lambda function to get only two decimals:
var_df["Solution"] = var_df["Variables"].apply(lambda item: "{:.2f}".format(float(item.varValue)))
# We do the same for the reduced cost:
var_df["Reduced cost"] = var_df["Variables"].apply(lambda item: "{:.2f}".format(float(item.dj)))


# We use the display function to represent the results:
display(var_df)


# we define a dictionary with the constraints:
const_dict = dict(model.constraints)
#We create a list of records from the dictionary and exclude the Expression to have a more compact solution.
con_df = pd.DataFrame.from_records(list(const_dict.items()), exclude=["Expression"], columns=["Constraint", "Expression"])

#Now we add columns for the solution, the slack and shadow price

con_df["Right Hand Side"] = con_df["Constraint"].apply(lambda item: "{:.2f}".format(-const_dict[item].constant))
con_df["Slack"] = con_df["Constraint"].apply(lambda item: "{:.2f}".format(const_dict[item].slack))
con_df["Shadow Price"] = con_df["Constraint"].apply(lambda item: "{:.2f}".format(const_dict[item].pi))

# And we display the results
display(con_df)

The value of the objective function is 6350.00

The following tables show the values obtained:

	Variables	Solution	Reduced cost
x	vars_x	12.00	0.00
y	vars_y	11.00	0.00

	Constraint	Right Hand Side	Slack	Shadow Price
0	Man_Power	40.00	5.00	-0.00
1	Machine_Operating_Time	45.00	-0.00	83.33
2	Marketing	12.00	-0.00	216.67

Analysis questions

1. Search for the exercise “The Good Carpenter” and try to model it.

2. Ask an AI assistant to model the problem and to provide a Python script to solve it. You can use this sample prompt:

Provide the problem model and a Python script using PuLP to solve the problem definition below assuming that decision variables are continuous: [Copy the problem definition from the exercise “The Good Carpenter” here]

3. Try the solution provided by the AI assistant. Does it work?

4. Is the solution provided by the AI assistant scalable? If the answer is not, ask it to provide a solution that is more scalable and explain its solution, using the following prompt template:

The Python script provided does not scale up well if the number of decision variables or the number of constraints increase. Can you provide a more scalable solution using PulP? Please describe the changes you made and motivate why

Did the AI assistant use the same methods explained in this tutorial?

5. Ask the assistant to use the problem you provided as an example to demonstrate how indices, summations, and “for all” expressions can make mathematical models more compact and scalable. Also, ask the assistant if these mathematical expressions are conceptually similar to the PuLP functions used in the scalable solution.

Solved exercises

The following notebooks include exercises solved with PuLP, using different solvers, do not forget to check them out!

• Making Chappie solved with CBC

• Blending problem solved with CBC

• Blending Craft Beer solved with CBC

• Blending Craft Beer solved with GRB

• Chappie II solved with CBC

• Manufacturing solar cell panels solved with CBC

• Manufacturing solar cell panels solved with GRB

• Petroleum Blending Solved with CBC