Pyomo Homework 1

Pyomo Homework 1#

# This code cell installs packages on Colab

import sys
if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
    import helper
    helper.easy_install()
    helper.install_glpk()
else:
    sys.path.insert(0, '../')
    import helper
helper.set_plotting_style()

## IMPORT LIBRARIES
import pyomo.environ as pyo
import pandas as pd

Special thanks to the Pyomo team for create these excercises as part of their excellent PyomoFest workshop.

1 Pyomo Fundamentals#

1.1 Knapsack example#

You want to fill a knapsack (a.k.a. bag). You can choose from a hammer, wrench, screwdriver, and towel. Each item has a different weight and value. You want to maximize the value (benefit) of the collection of items constrained by a total weight limit. Let’s formulate this as an optimization problem.

Sets

\[\mathcal{A} = \{\text{hammer},~\text{wrench},~\text{screwdriver},~\text{towel} \}\]

Parameters (Data)

Let \(b_i\) and \(w_i\) represent the benefit and weight of item \(i\), respectfully.

Item (\(i\))

Benefit (\(b_i\))

Weigth (\(w_i\))

hammer

8

5

wrench

3

7

screwdriver

6

4

towel

11

3

Let \(W_{max} = 14\) be the maximum weight.

Variables

Let \(x_i \in \{0,1\}\) (binary) represent whether or not we include item \(i\) in the knapsack. For now, we will consider only being able to choose either none or one of each item.

Objective and Constraints

\[\begin{split} \begin{equation} \begin{split} \max_{x} \quad & \sum_{i\in{\mathcal{A}}}b_i x_i \\ \text{s.t.} \quad & \sum_{i\in{\mathcal{A}}}w_ix_i \leq W_{max} \\ & x_i \in \{0,1\}, \quad \forall i \in \mathcal{A} \end{split} \end{equation} \end{split}\]

Pyomo

Solve the knapsack problem given below using GLPK and answer the following questions:

  1. Which items are acquired in the optimal solution?

  2. Why does this solution make sense? (Write ~2 sentences.)

A = ['hammer', 'wrench', 'screwdriver', 'towel']
b = {'hammer':8, 'wrench':3, 'screwdriver':6, 'towel':11}
w = {'hammer':5, 'wrench':7, 'screwdriver':4, 'towel':3}
W_max = 14

model = pyo.ConcreteModel()
model.x = pyo.Var( A, within=pyo.Binary )

model.obj = pyo.Objective(
    expr = sum( b[i]*model.x[i] for i in A ), 
    sense = pyo.maximize )

model.weight_con = pyo.Constraint(
    expr = sum( w[i]*model.x[i] for i in A ) <= W_max )

# Add your solution here

model.display()

Question Answers

  1. Fill in here

  2. Fill in here

1.2 Knapsack example with improve printing#

Complete the missing lines in the code below to produce formatted output: print the total weight, the value of the items selected (the objective), and the items acquired in the optimal solution. Note, the Pyomo value function should be used to get the floating point value of Pyomo modeling components (e.g., print(value(model.x[i])).

A = ['hammer', 'wrench', 'screwdriver', 'towel']
b = {'hammer':8, 'wrench':3, 'screwdriver':6, 'towel':11}
w = {'hammer':5, 'wrench':7, 'screwdriver':4, 'towel':3}
W_max = 14

model = pyo.ConcreteModel()
model.x = pyo.Var( A, within=pyo.Binary )

model.obj = pyo.Objective(
    expr = sum( b[i]*model.x[i] for i in A ), 
    sense = pyo.maximize )

model.weight_con = pyo.Constraint(
    expr = sum( w[i]*model.x[i] for i in A ) <= W_max )

opt = pyo.SolverFactory('glpk')
opt_success = opt.solve(model)

total_weight = sum( w[i]*pyo.value(model.x[i]) for i in A )
# Add your solution here

print('%12s %12s' % ('Item', 'Selected'))
print('=========================')
for i in A:
    # Add your solution here
print('-------------------------')

1.3 Changing data#

Using your code from Question 1.2, if we were to increase the value of the wrench, at what point would it become selected as part of the optimal solution?

# Add your solution here

Question Answer

Fill in here

1.4 Loading data from Excel#

In the code above, the data is hardcoded at the top of the file. Instead of hardcoding the data, use Python to load the data from a difference source. You may use Pandas to load data from ‘knapsack_data.xlsx’ into a dataframe. You will then need to write code to obtain a dictionary from the dataframe.

df_items = pd.read_excel(https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/data/knapsack_data.xlsx', sheet_name='data', header=0, index_col=0)
W_max = 14

A = df_items.index.tolist()
# Add your solution here

model = pyo.ConcreteModel()
model.x = pyo.Var( A, within=pyo.Binary )

model.obj = pyo.Objective(
    expr = sum( b[i]*model.x[i] for i in A ), 
    sense = pyo.maximize )

model.weight_con = pyo.Constraint(
    expr = sum( w[i]*model.x[i] for i in A ) <= W_max )

opt = pyo.SolverFactory('glpk', tee=True)
opt_success = opt.solve(model)

total_weight = sum( w[i]*pyo.value(model.x[i]) for i in A )
print('Total Weight:', total_weight)
print('Total Benefit:', pyo.value(model.obj))

print('%12s %12s' % ('Item', 'Selected'))
print('=========================')
for i in A:
    acquired = 'No'
    if pyo.value(model.x[i]) >= 0.5:
        acquired = 'Yes'
    print('%12s %12s' % (i, acquired))
print('-------------------------')

1.5 NLP vs. MIP#

Solve the knapsack problem with IPOPT instead of glpk. Print the solution values for model.x. What happened? Why?

Hint: Switch glpk to ipopt in the call to SolverFactory.

A = ['hammer', 'wrench', 'screwdriver', 'towel']
b = {'hammer':8, 'wrench':3, 'screwdriver':6, 'towel':11}
w = {'hammer':5, 'wrench':7, 'screwdriver':4, 'towel':3}
W_max = 14

model = pyo.ConcreteModel()
model.x = pyo.Var( A, within=pyo.Binary )

model.obj = pyo.Objective(
    expr = sum( b[i]*model.x[i] for i in A ), 
    sense = pyo.maximize )

model.weight_con = pyo.Constraint(
    expr = sum( w[i]*model.x[i] for i in A ) <= W_max )

# Add your solution here
opt_success = opt.solve(model)

model.pprint()

Question Answers

Fill in here

More Pyomo Fundamentals#

2.1 Knapsack problem with rules#

Rules are important for defining indexed constraints, however, they can also be used for single (i.e. scalar) constraints. Reimplement the knapsack model from Question 1.1 using rules for the objective and the constraints.

A = ['hammer', 'wrench', 'screwdriver', 'towel']
b = {'hammer':8, 'wrench':3, 'screwdriver':6, 'towel':11}
w = {'hammer':5, 'wrench':7, 'screwdriver':4, 'towel':3}
W_max = 14

model = pyo.ConcreteModel()
model.x = pyo.Var( A, within=pyo.Binary )

# Add your solution here

2.2 Integer formulation of the knapsack problem#

Consider again the knapsack problem. Assume now that we can acquire multiple items of the same type. In this new formulation, \(x_i\) is now an integer variable instead of a binary variable. One way to formulate this problem is as follows:

\[\begin{split} \begin{equation} \begin{split} \max_{x} \quad & \sum_{i\in{\mathcal{A}}}b_i x_i \\ \text{s.t.} \quad & \sum_{i\in{\mathcal{A}}}w_i x_i \leq W_{max} \\ & x_i=\sum_{j=0}^Njq_{i,j}, \quad \forall i \in \mathcal{A} \\ & 0 \leq x_i \leq N, \quad \forall i \in \mathcal{A} \\ & q_{i,j} \in \{0,1\}, \quad \forall i \in \mathcal{A}, j \in \{0,...,N\} \end{split} \end{equation} \end{split}\]

One could optionally add the following constraint select only one \(q_{i,j}\) for each \(i\), although it is not stricly neccessary to yield an integer solution. $$

(1)#\[\begin{equation} \sum_{j=0}^N q_{i,j} = 1, \quad \forall i \in \mathcal{A} \end{equation}\]

$$

Starting with your code from Question 2.1, implement this new formulation and solve. Is the solution surprising?

A = ['hammer', 'wrench', 'screwdriver', 'towel']
b = {'hammer':8, 'wrench':3, 'screwdriver':6, 'towel':11}
w = {'hammer':5, 'wrench':7, 'screwdriver':4, 'towel':3}
W_max = 14
N = range(6) # create a list from 0-5

model = pyo.ConcreteModel()

model.x = pyo.Var( A )
model.q = pyo.Var( A, N, within=pyo.Binary )

def obj_rule(m):
    return sum( b[i]*m.x[i] for i in A )
model.obj = pyo.Objective(rule=obj_rule, sense = pyo.maximize )

def weight_con_rule(m):
    return sum( w[i]*m.x[i] for i in A ) <= W_max
model.weight_con = pyo.Constraint(rule=weight_con_rule)

# Add your solution here

Question Answer

Fill in here