Wang Haihua
🚅 🚋😜 🚑 🚔
Goals of this notebook:
Learn some basic commands to solve an LP with Python.
These commands include LpProblem, LpVariable, and lpSum.
Python packages required: pulp
Additional resources: For more on PuLP, see https://pypi.org/project/PuLP/#description. See https://coin-or.github.io/pulp/CaseStudies/index.html for additional examples.
In its most general form, a linear program looks like
$$ \begin{array}{rcl} \max & c^\intercal x\\ \text{s.t.}& Ax &\le& b\\ & Bx &=& d\\ & Cx &\ge& f\\ & x & \in & \mathbb{R}^n. \end{array} $$There are three components to an LP: the variables $x$, the objective $\max ~ c^\intercal x$, and the constraints $Ax \le b$, $Bx = d$, and $Cx \ge f$. In order to code these using Python, we can follow the following six steps:
Step 1: Import Python's toolbox for solving LPs. This toolbox is called PuLP, and it contains many useful linear programming tools.
Step 2: Create an empty linear program. Intuitively, an empty linear program is Python's version of a sheet of paper on which we write the variables, the objetive, and the constraints.
Step 3: Add the variables $x$. The three components of the LP involve $x$, so we need to create these first.
Step 4: Add the objective $\max ~ c^\intercal x$.
Step 5: Add the constraints $Ax \le b, ~Bx = d$, and $Cx \ge f$.
Step 6: Solve the LP and print the results.
The following example will help us introduce the basic Python commands needed in Steps 1-6.
Alice wants to build a farm to produce corn. She can spend CHF 1000. It costs Alice CHF 3 to produce one kilogram of corn, and she can sell it for CHF 7. Alice can also buy additional farmland at a cost of CHF 100 per acre, and each acre can only grow 30 kilograms of corn. How many acres and kgs of corn should she buy to maximize profit?
Here is a model of Alice's problem.
$$ \begin{array}{rlcl} \max & 7 \times (\text{corn produced}) \\ \text{s.t.} & \text{corn produced} &\le& 30 \times (\text{acres purchased})\\ & 3 \times (\text{corn produced}) + 100 \times (\text{acres purchased}) &\le& 1000\\ &0 \le \text{corn produced}\\ &0 \le \text{acres purchased}. \end{array} $$Run the following line of code to import PuLP (you do not need to know how this line of code works).
Note: One way to run the code is to click in the box below and press the 'Run' button. Another way is to click in the box and press 'Shift + Enter'.
# Import the necessary LP toolbox from Python
from pulp import *
The function to create an empty linear program is LpProblem and it has three parts:
"Alice's Farm" - This is the name displayed when we print the linear program. You can choose any name you want.
LpMaximize - This makes the linear program a maximization problem.
my_LP - This is the name that code uses to identify our LP. You can choose any name you want.
Run the following line of code to create our linear program.
# Create an empty linear program
my_LP = LpProblem("Alice's_Farm", LpMaximize)
We need to create two variables: one for the amount of corn produced and one for the number of acres purchased.
Consider the variable for the amount of corn produced.
The function to create this variable is LpVariable and it has three components:
"Corn_produced" - This is the name displayed when we print the LP. You can choose any name you want.
"corn" - This is the name that the code uses to identify this variable. You can choose any name you want.
lowBound = 0 - This guarantees that the variable is lower bounded by 0 (i.e., that is it nonnegative).
Note: The command lowBound = 0 creates the inequality $\text{corn} \ge 0$.
Alternatively, we can add this as a constraint in Step 5.
However, this type of constraint is so common that PuLP has commands to add it at this step.
The LpVariable function can also be used to create a variable for the number of acres purchased.
Run the following line of code to add our variables.
# Create the variables
corn = LpVariable("Corn_produced", lowBound=0)
acres = LpVariable("Acres_purchased", lowBound=0)
In Step 2 that we created my_LP to be a maximization problem.
Therefore, we only need to add the objective function $c^\intercal x$.
The objective function for our problem is $7 \times \text{(corn produced)}$.
Using our variables from Step 3, the objective function becomes 7*corn.
We use the command += to add this to my_LP.
Run the following line of code to add this objective to my_LP.
# Add the objective
my_LP += 7*corn
The first constraint is $\text{corn produced} \le 30 \times (\text{acres purchased})$.
Using our variables this becomes corn <= 30* acres.
We use the command += to add this to my_LP.
Run the following line of code to create this constraint and add it to the my_LP.
# Add the first constraint
my_LP += corn <= 30*acres
The second constraint is $3 \times (\text{corn produced}) + 100 \times (\text{acres purchased}) \le 1000$.
Add this to our LP by running the following line of code.
# Add the second constraint
my_LP += 3*corn +100*acres <= 1000
Now our linear program is completely built! The following line of code will display your linear program so that you may check it.
# Display our linear program
print(my_LP)
Alice's_Farm: MAXIMIZE 7*Corn_produced + 0 SUBJECT TO _C1: - 30 Acres_purchased + Corn_produced <= 0 _C2: 100 Acres_purchased + 3 Corn_produced <= 1000 VARIABLES Acres_purchased Continuous Corn_produced Continuous
Note: After running the print(my_LP) command, the variables displayed will not explicitly state that the variables are lower bounded by zero.
Also, the variables will be displayed as Continuous. This can be ignored for now, but we will return to this when we consider discrete decision problems.
The following lines of code solves our linear program.
# Solve the linear program
my_LP.solve()
1
If everything went well, then the output should be 1.
This does not mean that the solution to my_LP equals 1, but rather the purpose of this output is just to tell us that everything went ok.
However, we (Alice and us) would like to know the optimal objective value and the optimal values of corn and acres.
Note : The possible output values are -3,-2,-1,0 or 1.
You can check what the different output values mean at
https://coin-or.github.io/pulp/technical/constants.html.
The optimal value of corn can be accessed using corn.value().
Similarly, the optimal value of acres can be accessed using acres.value().
The optimal objective value from my_LP is accessed using value(my_LP.objective).
Note :
The character %.2f is a formatting tool for rounding a decimal to two places.
It is not important for this tutorial to know how this formatting works.
For more see https://docs.python.org/2/library/string.html.
Run the following line of code to display the optimal values from my_LP.
# Print the optimal value and the variables values
opt_corn = corn.value()
print(f'Alice should produce {opt_corn:.2f} kilograms of corn.')
opt_acres = acres.value()
print(f'Alice should purchase {opt_acres:.2f} acres of land.')
opt_val = value(my_LP.objective)
print(f'Alice will have a profit of CHF {opt_val:.2f}.')
Alice should produce 157.89 kilograms of corn. Alice should purchase 5.26 acres of land. Alice will have a profit of CHF 1105.26.
Congratulations! We have successfully solved our first LP. If everything was run correctly, then the output should be
Alice should produce 157.89 kilograms of corn.
Alice should purchase 5.26 acres of land.
Alice will have a profit of CHF 1105.26.There are six main steps to solving a linear program. Moreover, the basic commands that we learned are already enough to solve many optimization problems. See https://coin-or.github.io/pulp/CaseStudies/a_blending_problem.html for another example.
In the next example, we show how to use the lpSum command in PuLP when dealing with multiple variables.
Alice's farm was a helpful example, but it was also simple with only two variables.
What if we had a more complicated LP with more variables?
Some LPs can have thousands of variables!
In Step 3 we could write the LpVariable command thousands of times on thousands of lines of code, but this would be tedious and prone to mistakes!
Similarly, writing down one constraint that involves one thousand variables would be painful.
In order to handle these complications, we can use for loops and the lpSum command in PuLP.
Let us illustrate this with another example.
Bob works for the electric company and must decide how to distribute electricity to the houses on his street. The street has seven houses, and each house requires at least 30 kWh per day. The company has two generators that can supply power to the houses. Generator 1 can supply at most 100 kWh per day and generator 2 can supply at most 150 kWh. The following price chart shows how much it costs (in CHF) each generator to supply 1 kWh to each house.
How should Bob supply power to the houses in order to minimize his daily cost?
The following is a mathematical model of Bob's problem. We use a variable $x_{g,h}$ to denote the number of kWh that generator $g$ sends to house $h$. For this model, we use the abbreviation $p_{g,h}$ to denote the price to send one kWh from generator $g$ to house $h$. For example $p_{2,6} = .23$.
$$ \begin{array}{rlll} \min & \displaystyle \sum_{g=1}^2 \sum_{h=1}^7 p_{g,h} \cdot x_{g,h}\\ \text{s.t.} & \displaystyle \sum_{h=1}^7 x_{1,h} \le 100 & \text{for each } g =1,2 &\text{[Generator 1 can produce at most 100 kWh]}\\ &\displaystyle \sum_{h=1}^7 x_{2,h} \le 150 & \text{for each } g =1,2 &\text{[Generator 2 can produce at most 150 kWh]}\\ & \displaystyle \sum_{g=1}^2 x_{g,h} \ge 30 & \text{for each } h =1,\dotsc, 7 &\text{[Each house needs at least 30 kWh]} \\ & x_{g,h} \ge 0 & \text{for each } g =1,2 \text{ and }h =1,\dotsc, 7. \end{array} $$Let us follow the steps outlined above.
Step 1 and Step 2 load PuLP and create an empty LP.
Note that this is a minimization problem, so we use LpMinimize in the function LpProblem.
We will also store the prices in a 2-dimensional Python list.
Run the following line of code to run these steps.
Note :
Lists are useful tools for storing information in Python, but we do not need to know how to manipulate them in this example.
All we need to know is that if we want to access a price, for instance the price from generator i to house j, then we use the code priceGen[i-1][j-1].
Recall that Python lists begin their indexing at 0 and not 1 (for a more precise indexing system, consider using Python dictionaries).
# Import the necessary LP toolbox from Python
from pulp import *
# Create an empty linear program
Electricity_LP = LpProblem("Bob's_electricity", LpMinimize)
# Load the price chart
# We use the notation gxhy to indicate generator x and house y
pricesGen = [[0.17, 0.25, 0.29, 0.16, 0.24, 0.20, 0.29],
[0.18, 0.21, 0.30, 0.20, 0.20, 0.23, 0.28]]
Step 3 is to create the variables.
There are 14 variables in total.
Here are two ways to define these variables.
One approach is to define each variable on a new line by using lines of code like variable_1 = LpVariable(), variable_2 = LpVariable(),....
However, this can be tedious and is prone to error.
Also, some problems have thousands of variables and adding one line for each variable takes too long.
A more convenient approach to defining these variables is to create the variables using a for loop. This allows us to turn 14 lines of code (or thousands of lines of code) into a few short lines.
We will take this second approach.
Run the following line of code to create the variables.
Note : We can store the LP variables in a Python list, tuple, or dictionary. Here, we store the variables in a 2-dimensional list.
# Create the variables
variables = [[LpVariable(f'g{gen}h{house}', lowBound=0)
for house in range(7)]
for gen in range(2)]
Step 4 is to add the objective function $\sum_{g=1}^2 \sum_{h=1}^7 p_{g,h} \cdot x_{g,h}$.
One approach to add the objective function is to write Electricity_LP += followed by all 14 terms in the objective function.
However, as was the case when we added variables, this can be tedious and lead to mistakes.
Fortunately, PuLP has a command called lpSum that allows us to add variables together using a for loop.
The command lpSum works like a for loop with the only difference being that the for statements come at the end.
Run the following line of code to add the objective.
# Add the objective
Electricity_LP += lpSum([pricesGen[gen][house]*variables[gen][house]
for gen in range(2)
for house in range(7)])
Step 5 is to add the constraints.
Luckily, we can use the lpSum command to help us again.
The lpSum command only adds the variables together.
It is left to us to add the inequalities <=100, <= 150, and >= 30 to the code.
Run the following line of code to add the constraints.
# Add the constraints
# The first set of constraints is for the generators
# Generator 1 can supply at most 100 kWh
Electricity_LP += lpSum(variables[0]) <= 100
# Generator 2 can supply at most 150 kWh
Electricity_LP += lpSum(variables[1]) <= 150
# The second set of constraints says each house must recieve at least 30 kWh
for house in range(7):
Electricity_LP += lpSum([variables[gen][house] for gen in range(2)]) >= 30
Great! Now we have everything we need to display and solve our LP.
# Display our linear program
print(Electricity_LP)
# Solve the linear program
Electricity_LP.solve()
# Print the optimal value and the variables values
opt_val = value(Electricity_LP.objective)
print(f'Bob must spend CHF {opt_val:.2f}.')
for gen in range(2):
for house in range(7):
opt_power = variables[gen][house].value()
print(f'Bob should send {opt_power:.2f} kWh from generator {gen+1} to house {house+1}.')
Bob's_electricity: MINIMIZE 0.17*g0h0 + 0.25*g0h1 + 0.29*g0h2 + 0.16*g0h3 + 0.24*g0h4 + 0.2*g0h5 + 0.29*g0h6 + 0.18*g1h0 + 0.21*g1h1 + 0.3*g1h2 + 0.2*g1h3 + 0.2*g1h4 + 0.23*g1h5 + 0.28*g1h6 + 0.0 SUBJECT TO _C1: g0h0 + g0h1 + g0h2 + g0h3 + g0h4 + g0h5 + g0h6 <= 100 _C2: g1h0 + g1h1 + g1h2 + g1h3 + g1h4 + g1h5 + g1h6 <= 150 _C3: g0h0 + g1h0 >= 30 _C4: g0h1 + g1h1 >= 30 _C5: g0h2 + g1h2 >= 30 _C6: g0h3 + g1h3 >= 30 _C7: g0h4 + g1h4 >= 30 _C8: g0h5 + g1h5 >= 30 _C9: g0h6 + g1h6 >= 30 VARIABLES g0h0 Continuous g0h1 Continuous g0h2 Continuous g0h3 Continuous g0h4 Continuous g0h5 Continuous g0h6 Continuous g1h0 Continuous g1h1 Continuous g1h2 Continuous g1h3 Continuous g1h4 Continuous g1h5 Continuous g1h6 Continuous Bob must spend CHF 45.50. Bob should send 10.00 kWh from generator 1 to house 1. Bob should send 0.00 kWh from generator 1 to house 2. Bob should send 30.00 kWh from generator 1 to house 3. Bob should send 30.00 kWh from generator 1 to house 4. Bob should send 0.00 kWh from generator 1 to house 5. Bob should send 30.00 kWh from generator 1 to house 6. Bob should send 0.00 kWh from generator 1 to house 7. Bob should send 20.00 kWh from generator 2 to house 1. Bob should send 30.00 kWh from generator 2 to house 2. Bob should send 0.00 kWh from generator 2 to house 3. Bob should send 0.00 kWh from generator 2 to house 4. Bob should send 30.00 kWh from generator 2 to house 5. Bob should send 0.00 kWh from generator 2 to house 6. Bob should send 30.00 kWh from generator 2 to house 7.
If everything ran correctly, then the output should say that Bob spends CHF 45.50.
A linear program can have multiple variables and constraints.
There are certain tools in Python such as dictionaries and the lpSum command to help us simplify our code when solving the LP.