Make me look good. Click on the cell below and press Ctrl+Enter.
from IPython.core.display import HTML
HTML(open('css/custom.css', 'r').read())
We learned the basics of using Pyomo to model and solve linear programs in Project 3.
Many of you have also used Pyomo in SA305, your course on linear programming.
In SA305, you may have learned about the benefits of using symbolic input parameters or set-based notation when formulating linear programs.
In this lesson, we'll go over an example of how to use set-based notation in Pyomo.
Example. (POMP Chapter 3) Consider the following warehouse location problem. Let $N$ be a set of candidate warehouse locations, and let $M$ be a set of customer locations. For each warehouse $n \in N$, the cost of delivering product to customer $m \in M$ is given by $d_{n,m}$. We want to determine the warehouse locations that will minimize the total cost of product delivery.
Below is a formulation of this problem. The binary variables $y_n$ are used to define whether or not a warehouse should be built: $y_n$ is 1 if warehouse $n$ is built and 0 otherwise. The variables $x_{n,m}$ indicate the fraction of demand for customer $m$ that is served by warehouse $n$.
Quick check. Is this a linear program?
Write your answer here. Double-click to edit.
No! Although the objective function is linear, and constraints (WL.2) - (W.5) are also linear, the variables $y_n$ are required to be binary. This requirement makes the above formulation an integer linear program. You'll learn more about these in SA405.
Nevertheless, we can still use this model to learn about how to use set-based notation in Pyomo.
In the code cell below, we begin by importing pyomo.environ
as pyo
, defining all the required data for the example, and defining a concrete model.
Note the use of tuples as keys for the dictionary d
. For example, in the code below, we have
d[('Harlingen', 'NYC')] = 1956
The idea is that d[(n, m)]
represents $d_{n,m}$ in the formulation above.
In Python, a tuple is a collection of items in round parentheses ()
, separated by commas.
For our purposes here, you can think of a tuple as a list-like object that you can use as a dictionary key.
For more details on lists vs. tuples, take a look at this article from Real Python.
# First import Pyomo
import pyomo.environ as pyo
# Define data
N = ['Harlingen', 'Memphis', 'Ashland']
M = ['NYC', 'LA', 'Chicago', 'Houston']
d = {
('Harlingen', 'NYC'): 1956,
('Harlingen', 'LA'): 1606,
('Harlingen', 'Chicago'): 1410,
('Harlingen', 'Houston'): 330,
('Memphis', 'NYC'): 1096,
('Memphis', 'LA'): 1792,
('Memphis', 'Chicago'): 531,
('Memphis', 'Houston'): 567,
('Ashland', 'NYC'): 485,
('Ashland', 'LA'): 2322,
('Ashland', 'Chicago'): 324,
('Ashland', 'Houston'): 1236
}
P = 2
# Define the model as a Concrete Model in Pyomo
model = pyo.ConcreteModel()
N
and M
defined above.# Define sets
model.N = pyo.Set(initialize=N, doc='Set of candidate warehouse locations')
model.M = pyo.Set(initialize=M, doc='Set of customer locations')
d
defined in the code above.# Define parameters
model.d = pyo.Param(model.N, model.M, initialize=d, doc='cost of delivering product to customer m from warehouse n')
Now, we define the decision variables $x_{n,m}$ for $n \in N$ and $m \in M$, and $y_n$ for $n \in N$. Note that
we use the keyword argument bounds
to restrict $x_{n,m}$ to values between 0 and 1, and
we use the keyword argument within
to specify the domain of the decision variable $y_n$ to be pyo.Binary
.
# Define decision variables
model.x = pyo.Var(model.N, model.M, bounds=(0,1), doc='fraction of demand for customer m served by warehouse n')
model.y = pyo.Var(model.N, within=pyo.Binary, doc='1 if warehouse n is open, 0 otherwise')
Next we use construction rules to define the objective function and constraints.
We can add the objective function (WL.1)
$$ \min_{x,y} \quad \sum_{n\in N} \sum_{m \in M} d_{n,m} x_{n,m} \qquad (\rm{WL.1}) $$
like this:
# Define objective (WL.1)
def obj_rule(model):
return sum(model.d[n,m] * model.x[n,m] for n in model.N for m in model.M)
model.obj = pyo.Objective(rule=obj_rule, sense=pyo.minimize)
Next, we can add the constraint (WL.2) that ensures each customer's demand is exactly satisfied
$$ \sum_{n\in N} x_{n,m} = 1 \quad \forall m \in M \qquad (\rm{WL.2}) $$
like this:
# Define "100% of each customer's demand satisfied" constraint (WL.2)
def one_per_cust_rule(model, m):
return sum(model.x[n,m] for n in model.N) == 1
model.one_per_cust = pyo.Constraint(model.M, rule=one_per_cust_rule, doc="Each customer's demand is fully met")
# Define "warehouse is active" constraint (WL.3)
def warehouse_active_rule(model, n, m):
return model.x[n,m] <= model.y[n]
model.warehouse_active = pyo.Constraint(model.N, model.M, rule=warehouse_active_rule, doc="Warehouse can only deliver if it is built")
# Define warehouse capacity constraint (WL.4)
def num_warehouses_rule(model):
return sum(model.y[n] for n in model.N) <= P
model.num_warehouses = pyo.Constraint(rule=num_warehouses_rule, doc="Can open at most P warehouses")
Now that we've finished writing the model in Pyomo, we need to solve it.
The code below sets the solver to "glpk"
, which corresponds to GLPK, the open source GNU Linear Programming Kit solver.
Then the code solves model
and saves the results object in results
.
Note that in this example, we set tee=False
, which means we will not see all the output from the solver.
# Set the solver to be GLPK
opt = pyo.SolverFactory("glpk")
# Call the solver
# Note that we set tee = False,
# so we don't see all the output from the solver
results = opt.solve(model, tee=False)
# Check to see if the solver termination condition was optimal
# If it's optimal, print the optimal objective value and solution
if str(results.solver.termination_condition) == "optimal":
print("Results\n")
# Print the objective function value
print(f"Minimum total cost is {pyo.value(model.obj):0.2f}.\n")
print("Open Warehouses")
print("---------------")
# Below we loop over all the elements of the set N
for n in model.N:
# Here we check to see if the decision variable value is
# not equal to 0 for this value of n
if model.y[n].value != 0:
# If the decision variable is not equal to 0, we
# print the value of n. This will tell us which
# warehouses the model decided to open.
print(f"{n:15s}")
print("\n")
print("Customer Warehouse Fraction of Demand Satistifed")
print("---------- ---------- -----------------------------")
# Below we loop over all the elements of the sets N and M
for n in model.N:
for m in model.M:
# Here we check to see if the decision variable value is
# not equal to 0 for this combination of n and m.
if model.x[n,m].value != 0:
# If the decision variable is not equal to 0, we
# print the values m, n, and the decision variable.
# This will tell us how much demand was satisfied for
# each customer from each servicing warehouse.
print(f"{m:10s} {n:10s} {model.x[n,m].value:<.2f}")
# Otherwise, print a warning message
else:
print("Solver termination condition was not optimal.")
Results Minimum total cost is 2745.00. Open Warehouses --------------- Harlingen Ashland Customer Warehouse Fraction of Demand Satistifed ---------- ---------- ----------------------------- LA Harlingen 1.00 Houston Harlingen 1.00 NYC Ashland 1.00 Chicago Ashland 1.00
Read pages 29-38 from POMP for an even more detailed treatment of this example.
There are a few minor differences in the code above over what is presented in the book.
pyo
as an alias when importing pyomo.environ
. This choice makes the model slightly more complicated, but it prevents pollution of the name space and is therefore recommended as good Python practice.Based on the example above, you should be able to solve Problems 1 and 2 below.
In Problem 2, you will grab the necessary data from an Excel workbook. You'll need to use of xlwings
, so you may want to refer back to Lesson 10 on working with spreadsheets.
Before proceeding with Problems 1 and 2, do Problem 0 first.
Problem 0. Read the sections on "Line Continuation", "Implicit Line Continuation", and "Explicit Line Continuation" in this article on Real Python.
"Line continuation" might sound scary, but it's just about syntax. In particular, what if you have a line of Python code that is really, really long - so long that it can't all be displayed at once? This article outlines some techniques you can write your code to avoid this problem and make your code more readable.
You'll find that the solutions to this lesson in particular, uses some of these techniques. You'll also find that we've already been using some of these techniques stealthily in previous lessons.
Problem 1. (Winston, Operations Research: Applications and Algorithms, Fourth Edition, Exercise 8 on page 93) Highland's TV-Radio Store must determine how many TVs and radios to keep in stock. A TV requires 10 sq. ft. of floorspace, whereas a radio requires 4 sq. ft. 200 sq. ft. of floorspace is available. A TV will earn Highland \$60 in profits and a radio will earn \$20. The store stocks only TVs and radios. Marketing requirements dictate that at least 60% of all appliances in stock be radios. TVs tie up \$200 in capital and radios tie up \$50. Highland wants to have at most \$3000 worth of capital tied up at any time.
Below is a linear program without set-based notation that can be used to maximize Highland's profit.
Linear program without set-based notation.
Decision variables. Let $t$ equal the number of TVs to keep in stock and let $r$ equal the number of radios to keep in stock.
Formulation. \begin{alignat}{3} \max \quad & 60t + 20r \\ \mbox{s.t.} \quad & 10t + 4r \le 200 &\qquad& (\rm{floorspace}) &\qquad& \\ & r \ge 0.60(t+r) & & (\rm{marketing}) \\ & 200t + 50r \le 3000 && (\rm{capital})\\ & t,r \ge 0 && \end{alignat}
Note that we can rewrite the marketing constraint as
\begin{alignat}{3} & 0.60 t - 0.40 r \le 0 &\qquad & (\rm{ marketing}) \end{alignat}Here is the same linear program, but with set-based notation. Make sure you understand how these two formulations are related.
Linear program with set-based notation.
Indices and sets. \begin{array}{ll} \mbox{$p \in P$} & \mbox{products}, P=\{\mbox{tv, radio} \} \\ \end{array}
Data. \begin{array}{ll} \mbox{profit}_p & \mbox{profit earned for product $p$} \\ \mbox{floorspace}_p & \mbox{sq. ft. of floorspace required for product $p$} \\ \mbox{marketing}_p & \mbox{marketing percentage for product $p$} \\ \mbox{capital}_p & \mbox{capital dollars tied up for product $p$} \\ \mbox{floorspace_avail} & \mbox{total sq of floorspace available} \\ \mbox{capital_avail} & \mbox{total dollars of capital available} \\ \end{array}
Decision variables. \begin{array}{ll} \mbox{PRODUCTS}_p & \mbox{number of items of product $p$ to put on the floor} \\ \end{array}
Formulation. \begin{alignat}{3} \max \quad & \sum_{p\in P} \mbox{profit}_p * \mbox{PRODUCTS}_p \\ \mbox{s.t.} \quad & \sum_{p\in P} \mbox{floorspace}_p * \mbox{PRODUCTS}_p \le \mbox{floorspace_avail} &\qquad& &\qquad& \\ & \sum_{p \in P} \mbox{marketing}_p * \mbox{PRODUCTS}_p \le 0 && \\ & \sum_{p \in P} \mbox{capital}_p * \mbox{PRODUCTS}_p \le \mbox{capital_avail} && \\ & \mbox{PRODUCTS}_p \ge 0 && \forall p\in P \end{alignat}
Implement the linear program with set-based notation in Pyomo, and solve it.
Verify that both $(t,r) = (6,35)$ and $(t,r) = (7,32)$ are feasible solutions. Which of these gives the higher profit?
Find how close the feasible integer solutions above are to the upper bound on the maximum possible profit given by the linear program's optimal value. You should report your answer as a percentage, e.g.
The answer with t = ... and r = ... is within 92.6% of the upper bound on the maximum profit.
# Import Pyomo
import pyomo.environ as pyo
# Define data
productslist = ['tv', 'radio']
profitdict = {'tv': 60, 'radio': 20}
floorspacedict = {'tv': 10, 'radio': 4}
marketingdict = {'tv': 0.6, 'radio': -0.4}
capitaldict = {'tv': 200, 'radio': 50}
floorspace_avail = 200
capital_avail = 3000
# Define the model as a Concrete Model in Pyomo
model = pyo.ConcreteModel()
# Define sets
model.P = pyo.Set(initialize=productslist, doc='Set of Products')
# Define parameters
model.profit = pyo.Param(model.P, initialize=profitdict,
doc='profit for item p')
model.floorspace = pyo.Param(model.P, initialize=floorspacedict,
doc='floorspace for item p')
model.marketing = pyo.Param(model.P, initialize=marketingdict,
doc='marketing for item p')
model.capital = pyo.Param(model.P, initialize=capitaldict,
doc='capital for item p')
# Define decision variables
model.PRODUCTS = pyo.Var(model.P, within=pyo.NonNegativeReals,
doc='number of items to put on floor')
# Define objective
# Note that the default objective sense is minimize in Pyomo
def profit_rule(model):
return sum(model.profit[p] * model.PRODUCTS[p] for p in model.P)
model.cost = pyo.Objective(rule=profit_rule, sense=pyo.maximize)
# Define constraints
def floorspace_rule(model):
return (
sum(model.PRODUCTS[p] * model.floorspace[p] for p in model.P)
<= floorspace_avail
)
model.avail_floorspace_rule = pyo.Constraint(rule=floorspace_rule,
doc='Floor space constraint')
def marketing_rule(model):
return (
sum(model.PRODUCTS[p] * model.marketing[p] for p in model.P)
<= 0
)
model.avail_marketing_rule = pyo.Constraint(rule=marketing_rule,
doc='Marketing constraint')
def capital_rule(model):
return (
sum(model.PRODUCTS[p] * model.capital[p] for p in model.P)
<= capital_avail
)
model.avail_capital_rule = pyo.Constraint(rule=capital_rule,
doc='Capital constraint')
# Set the solver to be GLPK
opt = pyo.SolverFactory("glpk")
# Call the solver
results = opt.solve(model, tee=False)
# If it's optimal, print the optimal objective value and solution
if str(results.solver.termination_condition) == "optimal":
print("Results\n")
print(f"Maximum Total Profit is {pyo.value(model.cost):0.2f} dollars.\n")
print("Product Number Proft Made ")
print("------------- ----------- -----------")
for p in model.P:
if model.PRODUCTS[p].value != 0:
# Note that multiple strings (plain or f-string) in a print statement
# (with nothing in between) will be concatenated
print(f"{p:13s} {model.PRODUCTS[p].value:11.2f} "
f"{model.PRODUCTS[p].value*model.profit[p]:11.2f}")
Results Maximum Total Profit is 1066.67 dollars. Product Number Proft Made ------------- ----------- ----------- tv 6.67 400.00 radio 33.33 666.67
# Note that multiple strings (plain or f-string) in a print statement
# (with nothing in between) will be concatenated
print(f'Stocking 7 televisions and 32 radios is feasible and gives '
f'profit equal to ${str(7 * 60 + 32 * 20)}.')
print(f'Stocking 6 televisions and 35 radios is feasible and gives '
f'profit equal to ${str(6 * 60 + 35 * 20)}.')
print(f'Both stocking plans are within just $6.67 of the LP relaxation.')
print(f'Either integer solution is at least {100 * 1060 / 1066.67:.2f}% optimal.')
Stocking 7 televisions and 32 radios is feasible and gives profit equal to $1060. Stocking 6 televisions and 35 radios is feasible and gives profit equal to $1060. Both stocking plans are within just $6.67 of the LP relaxation. Either integer solution is at least 99.37% optimal.
If we wanted to ensure that the number of each product is integer, we can change the type of decision variable on line 30 below.
While the change in code is simple, how the resulting model is solved is very different from the one above. You'll learn more about this in SA405.
# Import Pyomo
import pyomo.environ as pyo
# Define data
productslist = ['tv', 'radio']
profitdict = {'tv': 60, 'radio': 20}
floorspacedict = {'tv': 10, 'radio': 4}
marketingdict = {'tv': 0.6, 'radio': -0.4}
capitaldict = {'tv': 200, 'radio': 50}
floorspace_avail = 200
capital_avail = 3000
# Define the model as a Concrete Model in Pyomo
model = pyo.ConcreteModel()
# Define sets
model.P = pyo.Set(initialize=productslist, doc='Set of Products')
# Define parameters
model.profit = pyo.Param(model.P, initialize=profitdict,
doc='profit for item p')
model.floorspace = pyo.Param(model.P, initialize=floorspacedict,
doc='floorspace for item p')
model.marketing = pyo.Param(model.P, initialize=marketingdict,
doc='marketing for item p')
model.capital = pyo.Param(model.P, initialize=capitaldict,
doc='capital for item p')
# Define decision variables
model.PRODUCTS = pyo.Var(model.P, within=pyo.NonNegativeIntegers,
doc='number of items to put on floor')
# Define objective
# Note that the default objective sense is minimize in Pyomo
def profit_rule(model):
return sum(model.profit[p] * model.PRODUCTS[p] for p in model.P)
model.cost = pyo.Objective(rule=profit_rule, sense=pyo.maximize)
# Define constraints
def floorspace_rule(model):
return (
sum(model.PRODUCTS[p] * model.floorspace[p] for p in model.P)
<= floorspace_avail
)
model.avail_floorspace_rule = pyo.Constraint(rule=floorspace_rule,
doc='Floor space constraint')
def marketing_rule(model):
return (
sum(model.PRODUCTS[p] * model.marketing[p] for p in model.P)
<= 0
)
model.avail_marketing_rule = pyo.Constraint(rule=marketing_rule,
doc='Marketing constraint')
def capital_rule(model):
return (
sum(model.PRODUCTS[p] * model.capital[p] for p in model.P)
<= capital_avail
)
model.avail_capital_rule = pyo.Constraint(rule=capital_rule,
doc='Capital constraint')
# Set the solver to be GLPK
opt = pyo.SolverFactory("glpk")
# Call the solver
results = opt.solve(model, tee=False)
# If it's optimal, print the optimal objective value and solution
if str(results.solver.termination_condition) == "optimal":
print("Results\n")
print(f"Maximum Total Profit is {pyo.value(model.cost):0.2f} dollars.\n")
print("Product Number Proft Made ")
print("------------- ----------- -----------")
for p in model.P:
if model.PRODUCTS[p].value != 0:
# Note that multiple strings (plain or f-string) in a print statement
# (with nothing in between) will be concatenated
print(f"{p:13s} {model.PRODUCTS[p].value:11.2f} "
f"{model.PRODUCTS[p].value*model.profit[p]:11.2f}")
Results Maximum Total Profit is 1060.00 dollars. Product Number Proft Made ------------- ----------- ----------- tv 6.00 360.00 radio 35.00 700.00
Problem 2. A small engineering consulting firm has 3 senior designers available to work on the firm's 4 current projects over the next 2 weeks. Each designer has 80 hours to split among the projects, and the following table shows the manager's scoring (0 = nil to 100 = perfect) of the capability of each designer to contribute to each project, along with his estimate of the hours that each project will require.
\begin{array}{l|rrrr} & \mbox{Project} \\ \hline \mbox{Designer} & \mbox{1} & \mbox{2} & \mbox{3} & \mbox{4} \\ \hline \mbox{1} & 90 & 80 & 10 & 50 \\ \mbox{2} & 60 & 70 & 50 & 65 \\ \mbox{3} & 70 & 40 & 80 & 85 \\ \hline \mbox{Required} & 70 & 50 & 85 & 35 \end{array}Below is a linear program using set-based notation that maximizes the number of hours spent by the most qualified designers for a given project, while ensuring that all projects are completed according to the given requirements.
engineering_consulting.xlsx
, located in the same folder as this notebook. Use xlwings to read the data from the Excel Workbook on the Input
sheet.xlwings
to write the results back to the Excel spreadsheet on the Output
worksheet:A2
.B7:E9
.Indices and sets. \begin{array}{ll} \mbox{$d \in D$} & \mbox{designer}, D=\{\mbox{1, 2, 3} \} \\ \mbox{$p \in P$} & \mbox{project}, P=\{\mbox{1, 2, 3, 4} \} \\ \end{array}
Data. \begin{array}{ll} \mbox{score}_{d,p} & \mbox{manager's score for designer $d$ working on project $p$} \\ \mbox{required}_p & \mbox{manager's estimate of hours required to complete project $p$} \\ \mbox{hours}_d & \mbox{number of hours designer $d$ has to split among projects} \\ \end{array}
Decision variables. \begin{array}{ll} \mbox{$H_{d,p}$} & \mbox{number of hours worked by designer $d$ on project $p$} \\ \end{array}
Formulation. \begin{alignat}{3} \max \quad & \sum_{d \in D}\sum_{p\in P} \mbox{score}_{d,p} * H_{d,p} \\ \mbox{s.t.} \quad & \sum_{p\in P} H_{d,p} \le \mbox{hours}_d && \forall d \in D \\ & \sum_{d \in D} H_{d,p} \ge req_p && \forall p \in P\\ & H_{d,p} \ge 0 && \forall d \in D, p\in P \end{alignat}
# Import xlwings
import xlwings as xw
# Define the model as a Concrete Model in Pyomo
model = pyo.ConcreteModel()
# Define sets
model.D = pyo.Set(initialize=[1,2,3], doc='Set of designers')
model.P = pyo.Set(initialize=[1,2,3,4], doc='Set of projects')
# Define data and parameters
#
# Load data from Excel
#
# Load data workbook
wb = xw.Book(r'engineering_consulting.xlsx')
sinput = wb.sheets['Input']
soutput = wb.sheets['Output']
# Read ALL the entries from Excel into memory for score,
# then assign them to the scoredict dictionary in a nested loop
scoredict = {}
rng1 = sinput.range('C3').expand('table')
for d in model.D:
for p in model.P:
scoredict[d, p] = rng1[d - 1, p - 1].value
model.score = pyo.Param(
model.D, model.P, initialize=scoredict,
doc='score if designer d is working on project p'
)
# Read ALL the entries from Excel into memory for required,
# then assign them to the requireddict dictionary in a loop
requireddict = {}
rng2 = sinput.range('C7').expand('right')
for p in model.P:
requireddict[p] = rng2[p - 1].value
model.required = pyo.Param(
model.P, initialize=requireddict,
doc='manager estimate for time required to complete project p'
)
# Below we assign a value of 80 for EACH designer.
# This could change based on input from the client,
# so having it be "easy" to change as a parameter
# is potentially helpful for later.
model.hours = pyo.Param(
model.D, initialize=80,
doc='number of hours designer d has to split among projects'
)
# Define decision variables
model.H = pyo.Var(
model.D, model.P, within=pyo.NonNegativeReals,
doc='number of hours worked by designer d on project p'
)
# Define objective
# Recall that default sense is minimize in Pyomo
def qualified_rule(model):
return sum(model.score[d, p] * model.H[d, p]
for d in model.D for p in model.P)
model.cost = pyo.Objective(rule=qualified_rule, sense=pyo.maximize)
# Define constraints
def hours_rule(model, d):
return sum(model.H[d, p] for p in model.P) <= model.hours[d]
model.avail_hours_rule = pyo.Constraint(model.D, rule=hours_rule,
doc='Designer hours bound')
def project_rule(model, p):
return sum(model.H[d, p] for d in model.D) >= model.required[p]
model.req_project_rule = pyo.Constraint(model.P, rule=project_rule,
doc='Project hours bound')
# Set the solver to be GLPK
opt = pyo.SolverFactory("glpk")
# Call the solver
results = opt.solve(model, tee=False)
# If it's optimal, print the optimal objective value and solution
if str(results.solver.termination_condition)=="optimal":
print("Results\n")
print(f"Maximum Objective is {pyo.value(model.cost):0.2f}.\n")
print(" Project")
print("Designer 1 2 3 4")
print("-------- ----- ----- ----- -----")
for d in model.D:
print(f"{d:8d} {model.H[d, 1].value:5.1f} "
f"{model.H[d, 2].value:5.1f} {model.H[d, 3].value:5.1f} "
f"{model.H[d, 4].value:5.1f}")
# Now use xlwings to write output back to spreadsheet
# Clear current output in worksheet
wb.sheets['Output'].range('A2').clear()
wb.sheets['Output'].range('B7:E9').clear()
# Output objective value to Excel
soutput.range('A2').value = pyo.value(model.cost)
# Output solution to Excel
counter = 0
for d in model.D:
soutput.range(f'B{counter + 7}').value = model.H[d,1].value
soutput.range(f'C{counter + 7}').value = model.H[d,2].value
soutput.range(f'D{counter + 7}').value = model.H[d,3].value
soutput.range(f'E{counter + 7}').value = model.H[d,4].value
counter = counter + 1
Results Maximum Objective is 18825.00. Project Designer 1 2 3 4 -------- ----- ----- ----- ----- 1 70.0 10.0 0.0 0.0 2 0.0 40.0 5.0 35.0 3 0.0 0.0 80.0 0.0
The problem below is supplemental material. You should only work on it after you have completed the two probelms above.
Problem 3. (Winston, Operations Research: Applications and Algorithms, Fourth Edition, Example 7 on page 72)
A post office requires different numbers of full-time employees on different days of the week. The number of full-time employees required on each day is given in the Excel workbook post_office.xlsx
, located in the same folder as this notebook. Union rules state that each full-time employee must work five consecutive days and then receive two days off. For example, an employee who works Monday to Friday must be off on Saturday and Sunday. The post office wants to meet its daily requirements using only full-time employees.
Sheet1
.# Let X[0] = number of workers starting on Monday, ...,
# X[6] = number of workers starting on Sunday.
# Then want to minimize X[0] + X[1] + X[2] + X[3] + X[4] + X[5] + X[6]
# subject to:
# X[0] + X[3] + X[4] + X[5] + X[6] >= MON workers required
# X[0] + X[1] + X[4] + X[5] + X[6] >= TUE workers required
# X[0] + X[1] + X[2] + X[5] + X[6] >= WED workers required
# X[0] + X[1] + X[2] + X[3] + X[6] >= THU workers required
# X[0] + X[1] + X[2] + X[3] + X[4] >= FRI workers required
# X[1] + X[2] + X[3] + X[4] + X[5] >= SAT workers required
# X[2] + X[3] + X[4] + X[5] + X[6] >= SUN workers required
# X[0], ..., X[6] >= 0
# Define the model as a Concrete Model in Pyomo
model = pyo.ConcreteModel()
# Define sets
model.D = pyo.Set(initialize=[0, 1, 2, 3, 4, 5, 6],
doc='Set of days of the week -> 0 is Monday')
# Define data and parameters
# Read data from Excel
wb = xw.Book('post_office.xlsx')
sht = wb.sheets('Sheet1')
reqs = sht.range('B3').expand('down').value
# Read ALL the entries from Excel for number of
# required workers on each day of the week and
# then assign them to the reqsdict dictionary
reqsdict = {}
rng1 = sht.range('B3').expand('down')
for d in model.D:
reqsdict[d] = rng1[d].value
model.reqs = pyo.Param(model.D, initialize=reqsdict,
doc='number of full time employees required on day d')
# Define decision variables
model.X = pyo.Var(model.D, within=pyo.NonNegativeReals,
doc='number of workers starting on day d')
# Define objective
# Recall that default sense is minimize in Pyomo
def total_workers_rule(model):
return sum(model.X[d] for d in model.D)
model.cost = pyo.Objective(rule=total_workers_rule, sense=pyo.minimize)
# Define constraints
def monday_rule(model):
return model.X[0] + model.X[3] + model.X[4] + model.X[5] + model.X[6] \
>= model.reqs[0]
model.avail_monday_rule = pyo.Constraint(rule=monday_rule,
doc='Obey Monday workers constraint')
def tuesday_rule(model):
return model.X[0] + model.X[1] + model.X[4] + model.X[5] + model.X[6] \
>= model.reqs[1]
model.avail_tuesday_rule = pyo.Constraint(rule=tuesday_rule,
doc='Obey Tuesday workers constraint')
def wednesday_rule(model):
return model.X[0] + model.X[1] + model.X[2] + model.X[5] + model.X[6] \
>= model.reqs[2]
model.avail_wednesday_rule = pyo.Constraint(rule=monday_rule,
doc='Obey Wednesday workers constraint')
def thursday_rule(model):
return model.X[0] + model.X[1] + model.X[2] + model.X[3] + model.X[6] \
>= model.reqs[3]
model.avail_thursday_rule = pyo.Constraint(rule=thursday_rule,
doc='Obey Thursday workers constraint')
def friday_rule(model):
return model.X[0] + model.X[1] + model.X[2] + model.X[3] + model.X[4] \
>= model.reqs[4]
model.avail_friday_rule = pyo.Constraint(rule=friday_rule,
doc='Obey Friday workers constraint')
def saturday_rule(model):
return model.X[1] + model.X[2] + model.X[3] + model.X[4] + model.X[5] \
>= model.reqs[5]
model.avail_saturday_rule = pyo.Constraint(rule=saturday_rule,
doc='Obey Saturday workers constraint')
def sunday_rule(model):
return model.X[2] + model.X[3] + model.X[4] + model.X[5] + model.X[6] \
>= model.reqs[6]
model.avail_sunday_rule = pyo.Constraint(rule=sunday_rule,
doc='Obey Sunday workers constraint')
# Set the solver to be GLPK
opt = pyo.SolverFactory("glpk")
# Call the solver
results = opt.solve(model, tee=False)
# If it's optimal, print the optimal objective value and solution
if str(results.solver.termination_condition) == "optimal":
# display output
model.X.display()
# Now use xlwings to write output back to spreadsheet
# Clear current output in worksheet
sht.range('C3:C9').clear()
# Output objective value to Excel
sht.range('B11').value = pyo.value(model.cost)
# Output solution to Excel
counter = 0
for d in model.D:
sht.range(f'C{counter + 3}').value = model.X[d].value
counter = counter + 1
X : number of workers starting on day d Size=7, Index=D Key : Lower : Value : Upper : Fixed : Stale : Domain 0 : 0 : 2.0 : None : False : False : NonNegativeReals 1 : 0 : 7.0 : None : False : False : NonNegativeReals 2 : 0 : 0.0 : None : False : False : NonNegativeReals 3 : 0 : 9.0 : None : False : False : NonNegativeReals 4 : 0 : 6.0 : None : False : False : NonNegativeReals 5 : 0 : 0.0 : None : False : False : NonNegativeReals 6 : 0 : 0.0 : None : False : False : NonNegativeReals
Now solve the same problem but use the data on Sheet2
. What difficulties do you anticipate in using the optimal solution to staff the post office? How can you modify the Pyomo model to provide a solution that the post office can actually use for this second set of data?
# Let X[0] = number of workers starting on Monday, ...,
# X[6] = number of workers starting on Sunday.
# Then want to minimize X[0] + X[1] + X[2] + X[3] + X[4] + X[5] + X[6]
# subject to:
# X[0] + X[3] + X[4] + X[5] + X[6] >= MON workers required
# X[0] + X[1] + X[4] + X[5] + X[6] >= TUE workers required
# X[0] + X[1] + X[2] + X[5] + X[6] >= WED workers required
# X[0] + X[1] + X[2] + X[3] + X[6] >= THU workers required
# X[0] + X[1] + X[2] + X[3] + X[4] >= FRI workers required
# X[1] + X[2] + X[3] + X[4] + X[5] >= SAT workers required
# X[2] + X[3] + X[4] + X[5] + X[6] >= SUN workers required
# X[0], ..., X[6] >= 0
# Define the model as a Concrete Model in Pyomo
model = pyo.ConcreteModel()
# Define sets
model.D = pyo.Set(initialize=[0, 1, 2, 3, 4, 5, 6],
doc='Set of days of the week -> 0 is Monday')
# Define data and parameters
# Read data from Excel
wb = xw.Book('post_office.xlsx')
sht = wb.sheets('Sheet2')
reqs = sht.range('B3').expand('down').value
# Read ALL the entries from Excel for number of
# required workers on each day of the week and
# then assign them to the reqsdict dictionary
reqsdict = {}
rng1 = sht.range('B3').expand('down')
for d in model.D:
reqsdict[d] = rng1[d].value
model.reqs = pyo.Param(model.D, initialize=reqsdict,
doc='number of full time employees required on day d')
# Define decision variables
model.X = pyo.Var(model.D, within=pyo.NonNegativeReals,
doc='number of workers starting on day d')
# Define objective
# Recall that default sense is minimize in Pyomo
def total_workers_rule(model):
return sum(model.X[d] for d in model.D)
model.cost = pyo.Objective(rule=total_workers_rule, sense=pyo.minimize)
# Define constraints
def monday_rule(model):
return model.X[0] + model.X[3] + model.X[4] + model.X[5] + model.X[6] \
>= model.reqs[0]
model.avail_monday_rule = pyo.Constraint(rule=monday_rule,
doc='Obey Monday workers constraint')
def tuesday_rule(model):
return model.X[0] + model.X[1] + model.X[4] + model.X[5] + model.X[6] \
>= model.reqs[1]
model.avail_tuesday_rule = pyo.Constraint(rule=tuesday_rule,
doc='Obey Tuesday workers constraint')
def wednesday_rule(model):
return model.X[0] + model.X[1] + model.X[2] + model.X[5] + model.X[6] \
>= model.reqs[2]
model.avail_wednesday_rule = pyo.Constraint(rule=monday_rule,
doc='Obey Wednesday workers constraint')
def thursday_rule(model):
return model.X[0] + model.X[1] + model.X[2] + model.X[3] + model.X[6] \
>= model.reqs[3]
model.avail_thursday_rule = pyo.Constraint(rule=thursday_rule,
doc='Obey Thursday workers constraint')
def friday_rule(model):
return model.X[0] + model.X[1] + model.X[2] + model.X[3] + model.X[4] \
>= model.reqs[4]
model.avail_friday_rule = pyo.Constraint(rule=friday_rule,
doc='Obey Friday workers constraint')
def saturday_rule(model):
return model.X[1] + model.X[2] + model.X[3] + model.X[4] + model.X[5] \
>= model.reqs[5]
model.avail_saturday_rule = pyo.Constraint(rule=saturday_rule,
doc='Obey Saturday workers constraint')
def sunday_rule(model):
return model.X[2] + model.X[3] + model.X[4] + model.X[5] + model.X[6] \
>= model.reqs[6]
model.avail_sunday_rule = pyo.Constraint(rule=sunday_rule,
doc='Obey Sunday workers constraint')
# Set the solver to be GLPK
opt = pyo.SolverFactory("glpk")
# Call the solver
results = opt.solve(model, tee=False)
# If it's optimal, print the optimal objective value and solution
if str(results.solver.termination_condition) == "optimal":
# display output
model.X.display()
# Now use xlwings to write output back to spreadsheet
# Clear current output in worksheet
sht.range('C3:C9').clear()
# Output objective value to Excel
sht.range('B11').value = pyo.value(model.cost)
# Output solution to Excel
counter = 0
for d in model.D:
sht.range(f'C{counter + 3}').value = model.X[d].value
counter = counter + 1
X : number of workers starting on day d Size=7, Index=D Key : Lower : Value : Upper : Fixed : Stale : Domain 0 : 0 : 5.66666666666667 : None : False : False : NonNegativeReals 1 : 0 : 4.66666666666667 : None : False : False : NonNegativeReals 2 : 0 : 0.0 : None : False : False : NonNegativeReals 3 : 0 : 8.66666666666667 : None : False : False : NonNegativeReals 4 : 0 : 2.66666666666667 : None : False : False : NonNegativeReals 5 : 0 : 0.0 : None : False : False : NonNegativeReals 6 : 0 : 0.0 : None : False : False : NonNegativeReals
Now fix the domain of the decision variables to be NonNegativeIntegers
so we can get a solution the post office can actually use.
# Let X[0] = number of workers starting on Monday, ...,
# X[6] = number of workers starting on Sunday.
# Then want to minimize X[0] + X[1] + X[2] + X[3] + X[4] + X[5] + X[6]
# subject to:
# X[0] + X[3] + X[4] + X[5] + X[6] >= MON workers required
# X[0] + X[1] + X[4] + X[5] + X[6] >= TUE workers required
# X[0] + X[1] + X[2] + X[5] + X[6] >= WED workers required
# X[0] + X[1] + X[2] + X[3] + X[6] >= THU workers required
# X[0] + X[1] + X[2] + X[3] + X[4] >= FRI workers required
# X[1] + X[2] + X[3] + X[4] + X[5] >= SAT workers required
# X[2] + X[3] + X[4] + X[5] + X[6] >= SUN workers required
# X[0], ..., X[6] >= 0
# Define the model as a Concrete Model in Pyomo
model = pyo.ConcreteModel()
# Define sets
model.D = pyo.Set(initialize=[0, 1, 2, 3, 4, 5, 6],
doc='Set of days of the week -> 0 is Monday')
# Define data and parameters
# Read data from Excel
wb = xw.Book('post_office.xlsx')
sht = wb.sheets('Sheet2')
reqs = sht.range('B3').expand('down').value
# Read ALL the entries from Excel for number of
# required workers on each day of the week and
# then assign them to the reqsdict dictionary
reqsdict = {}
rng1 = sht.range('B3').expand('down')
for d in model.D:
reqsdict[d] = rng1[d].value
model.reqs = pyo.Param(model.D, initialize=reqsdict,
doc='number of full time employees required on day d')
# Define decision variables
model.X = pyo.Var(model.D, within=pyo.NonNegativeIntegers,
doc='number of workers starting on day d')
# Define objective
# Recall that default sense is minimize in Pyomo
def total_workers_rule(model):
return sum(model.X[d] for d in model.D)
model.cost = pyo.Objective(rule=total_workers_rule, sense=pyo.minimize)
# Define constraints
def monday_rule(model):
return model.X[0] + model.X[3] + model.X[4] + model.X[5] + model.X[6] \
>= model.reqs[0]
model.avail_monday_rule = pyo.Constraint(rule=monday_rule,
doc='Obey Monday workers constraint')
def tuesday_rule(model):
return model.X[0] + model.X[1] + model.X[4] + model.X[5] + model.X[6] \
>= model.reqs[1]
model.avail_tuesday_rule = pyo.Constraint(rule=tuesday_rule,
doc='Obey Tuesday workers constraint')
def wednesday_rule(model):
return model.X[0] + model.X[1] + model.X[2] + model.X[5] + model.X[6] \
>= model.reqs[2]
model.avail_wednesday_rule = pyo.Constraint(rule=monday_rule,
doc='Obey Wednesday workers constraint')
def thursday_rule(model):
return model.X[0] + model.X[1] + model.X[2] + model.X[3] + model.X[6] \
>= model.reqs[3]
model.avail_thursday_rule = pyo.Constraint(rule=thursday_rule,
doc='Obey Thursday workers constraint')
def friday_rule(model):
return model.X[0] + model.X[1] + model.X[2] + model.X[3] + model.X[4] \
>= model.reqs[4]
model.avail_friday_rule = pyo.Constraint(rule=friday_rule,
doc='Obey Friday workers constraint')
def saturday_rule(model):
return model.X[1] + model.X[2] + model.X[3] + model.X[4] + model.X[5] \
>= model.reqs[5]
model.avail_saturday_rule = pyo.Constraint(rule=saturday_rule,
doc='Obey Saturday workers constraint')
def sunday_rule(model):
return model.X[2] + model.X[3] + model.X[4] + model.X[5] + model.X[6] \
>= model.reqs[6]
model.avail_sunday_rule = pyo.Constraint(rule=sunday_rule,
doc='Obey Sunday workers constraint')
# Set the solver to be GLPK
opt = pyo.SolverFactory("glpk")
# Call the solver
results = opt.solve(model, tee=False)
# If it's optimal, print the optimal objective value and solution
if str(results.solver.termination_condition) == "optimal":
# display output
model.X.display()
# Now use xlwings to write output back to spreadsheet
# Clear current output in worksheet
sht.range('C3:C9').clear()
# Output objective value to Excel
sht.range('B11').value = pyo.value(model.cost)
# Output solution to Excel
counter = 0
for d in model.D:
sht.range(f'C{counter + 3}').value = model.X[d].value
counter = counter + 1
X : number of workers starting on day d Size=7, Index=D Key : Lower : Value : Upper : Fixed : Stale : Domain 0 : 0 : 6.0 : None : False : False : NonNegativeIntegers 1 : 0 : 4.0 : None : False : False : NonNegativeIntegers 2 : 0 : 0.0 : None : False : False : NonNegativeIntegers 3 : 0 : 9.0 : None : False : False : NonNegativeIntegers 4 : 0 : 3.0 : None : False : False : NonNegativeIntegers 5 : 0 : 0.0 : None : False : False : NonNegativeIntegers 6 : 0 : 0.0 : None : False : False : NonNegativeIntegers