Pyomo Fundamentals - Pyomo Tutorial

Open this notebook in Colab.

This notebook is part of the SECQUOIA Research Group Pyomo Tutorial. Install GLPK before running the mixed-integer examples locally; see Setup and Solvers for solver notes.

Learning objectives¶

By the end of this chapter, you should be able to:

build Pyomo models with variables, objectives, and constraints;
load data from Python structures and spreadsheets;
compare continuous and mixed-integer formulations;
use mutable parameters and integer cuts to explore alternative solutions;
formulate a lot-sizing model from algebraic notation.

import pyomo.environ as pyo
import pandas as pd

1. Knapsack Modeling Basics¶

1.1. Knapsack example¶

Solve the knapsack problem shown in the tutorial using a Python environment or by running python knapsack.py from the command line. Which items are acquired in the optimal solution? What is the value of the selected items?

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)

model.display()

Model unknown

  Variables:
    x : Size=4, Index={hammer, wrench, screwdriver, towel}
        Key         : Lower : Value : Upper : Fixed : Stale : Domain
             hammer :     0 :   1.0 :     1 : False : False : Binary
        screwdriver :     0 :   1.0 :     1 : False : False : Binary
              towel :     0 :   1.0 :     1 : False : False : Binary
             wrench :     0 :   0.0 :     1 : False : False : Binary

  Objectives:
    obj : Size=1, Index=None, Active=True
        Key  : Active : Value
        None :   True :  25.0

  Constraints:
    weight_con : Size=1
        Key  : Lower : Body : Upper
        None :  None : 12.0 :  14.0

1.2. Knapsack with improved printing¶

The knapsack.py in the tutorial uses model.pprint() to see the value of the solution variables. Starting with the code in knapsack_print_incomplete.py, complete the missing lines to produce formatted output. Note that the Pyomo value function should be used to get the floating point value of Pyomo modeling components (e.g., print(value(model.x[i]))). Also print the value of the items selected (the objective), and the total weight. (A solution can be found in knapsack_print_soln.py).

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 )
print('Total Weight:', total_weight)
print('Total Benefit:', pyo.value(model.obj))

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

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 )
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('-------------------------')

Total Weight: 12.0
Total Benefit: 25.0
        Item     Selected
=========================
      hammer          Yes
      wrench           No
 screwdriver          Yes
       towel          Yes
-------------------------

1.3. Changing data¶

If we were to increase the value of the wrench, at what point would it become selected as part of the optimal solution? (A solution can be found in knapsack_wrench_soln.py.)

A = ['hammer', 'wrench', 'screwdriver', 'towel']
b = {'hammer':8, 'wrench':8, '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 )
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('-------------------------')

Total Weight: 14.0
Total Benefit: 25.0
        Item     Selected
=========================
      hammer           No
      wrench          Yes
 screwdriver          Yes
       towel          Yes
-------------------------

1.4. Loading data from Excel¶

In the knapsack example shown in the tutorial slides, the data is hardcoded at the top of the file. Instead of hardcoding the data, use Python to load the data from a different source. You can start from the file knapsack_pandas_excel_incomplete.py. (A solution that uses pandas to load the data from Excel is shown in knapsack_pandas_excel_soln.py.)

df_items = pd.read_excel('knapsack_data.xlsx', sheet_name='data', header=0, index_col=0)
W_max = 14

A = df_items.index.tolist()
b = df_items['Benefit'].to_dict()
w = df_items['Weight'].to_dict()

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 )
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('-------------------------')

df_items = pd.read_excel('knapsack_data.xlsx', sheet_name='data', header=0, index_col=0)
W_max = 14

A = df_items.index.tolist()
b = df_items['Benefit'].to_dict()
w = df_items['Weight'].to_dict()

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 )
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('-------------------------')

Total Weight: 12.0
Total Benefit: 25.0
        Item     Selected
=========================
      hammer          Yes
      wrench           No
 screwdriver          Yes
       towel          Yes
-------------------------

1.5. NLP vs MIP¶

Solve the knapsack problem with IPOPT instead of glpk. (Hint: switch glpk to ipopt in the call SolverFactory.) Print the solution values for model.x. What happened? Why?

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('ipopt')
opt_success = opt.solve(model)

model.pprint()

1 Var Declarations
    x : Size=4, Index={hammer, wrench, screwdriver, towel}
        Key         : Lower : Value              : Upper : Fixed : Stale : Domain
             hammer :     0 : 1.0000000095722053 :     1 : False : False : Binary
        screwdriver :     0 : 1.0000000094153156 :     1 : False : False : Binary
              towel :     0 :  1.000000009742034 :     1 : False : False : Binary
             wrench :     0 : 0.2857142884855868 :     1 : False : False : Binary

1 Objective Declarations
    obj : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : maximize : 8*x[hammer] + 3*x[wrench] + 6*x[screwdriver] + 11*x[towel]

1 Constraint Declarations
    weight_con : Size=1, Index=None, Active=True
        Key  : Lower : Body                                                      : Upper : Active
        None :  -Inf : 5*x[hammer] + 7*x[wrench] + 4*x[screwdriver] + 3*x[towel] :  14.0 :   True

3 Declarations: x obj weight_con

2. Rules and Integer Decisions¶

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. Starting with knapsack.py, reimplement the model using rules for the objective and the constraints. (A solution can be found in knapsack_rules_soln.py.)

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 )

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)

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

model.pprint()

1 Var Declarations
    x : Size=4, Index={hammer, wrench, screwdriver, towel}
        Key         : Lower : Value : Upper : Fixed : Stale : Domain
             hammer :     0 :   1.0 :     1 : False : False : Binary
        screwdriver :     0 :   1.0 :     1 : False : False : Binary
              towel :     0 :   1.0 :     1 : False : False : Binary
             wrench :     0 :   0.0 :     1 : False : False : Binary

1 Objective Declarations
    obj : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : maximize : 8*x[hammer] + 3*x[wrench] + 6*x[screwdriver] + 11*x[towel]

1 Constraint Declarations
    weight_con : Size=1, Index=None, Active=True
        Key  : Lower : Body                                                      : Upper : Active
        None :  -Inf : 5*x[hammer] + 7*x[wrench] + 4*x[screwdriver] + 3*x[towel] :  14.0 :   True

3 Declarations: x obj weight_con

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{aligned} \max_{q,x} &~ \sum_{i \in A} v_i x_i \\ \text{s.t.} &~ \sum_{i \in A} w_i x_i \leq W_{\max} \\ &~ x_i - \sum_{j=0}^N j q_{i,j} = 0 \quad \forall i \in A \\ &~ 0 \leq x \leq N \\ &~ q_{i,j} \in \{0,1\} \quad~~ \forall i \in A, j \in \{0,1,\dots,N\} \end{aligned}

(1)

Starting with knapsack_rules.py, implement this new formulation and solve. Is the solution surprising? (A solution can be found in knapsack_integer_soln.py.)

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)

def x_integer_rule(m, i):
    return m.x[i] == sum( j*m.q[i,j] for j in N )
model.x_integer = pyo.Constraint(A, rule=x_integer_rule)

opt = pyo.SolverFactory('glpk')
result_obj = 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:
    print('%12s %12s' % (i, pyo.value(model.x[i])))
print('-------------------------')

Total Weight: 12.0
Total Benefit: 44.0
        Item   # Selected
=========================
      hammer          0.0
      wrench          0.0
 screwdriver          0.0
       towel          4.0
-------------------------

3. Model Extensions¶

3.1. Using decorator notation for rules¶

In the slides, we saw an alternative notation for declaring and defining Pyomo components using decorators. Starting with the warehouse location problem in warehouse_location_decorator_incomplete.py change the model to use the decorator notation. (A solution for this problem can be found in warehouse_location_decorator_soln.py.)

# warehouse_location.py: Warehouse location determination problem
model = pyo.ConcreteModel(name="(WL)")

W = ['Harlingen', 'Memphis', 'Ashland']
C = ['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

model.x = pyo.Var(W, C, bounds=(0,1))
model.y = pyo.Var(W, within=pyo.Binary)

@model.Objective()
def obj(m):
    return sum(d[w,c]*m.x[w,c] for w in W for c in C)

@model.Constraint(C)
def one_per_cust(m, c):
    return sum(m.x[w,c] for w in W) == 1

@model.Constraint(W, C)
def warehouse_active(m, w, c):
    return m.x[w,c] <= m.y[w]

@model.Constraint()
def num_warehouses(m):
    return sum(m.y[w] for w in W) <= P

pyo.SolverFactory('glpk').solve(model)

model.y.pprint()
model.x.pprint()

model = pyo.ConcreteModel(name="(WL)")

W = ['Harlingen', 'Memphis', 'Ashland']
C = ['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

model.x = pyo.Var(W, C, bounds=(0,1))
model.y = pyo.Var(W, within=pyo.Binary)

@model.Objective()
def obj(m):
    return sum(d[w,c]*m.x[w,c] for w in W for c in C)

@model.Constraint(C)
def one_per_cust(m, c):
    return sum(m.x[w,c] for w in W) == 1

@model.Constraint(W,C)
def warehouse_active(m, w, c):
    return m.x[w,c] <= m.y[w]

@model.Constraint()
def num_warehouses(m):
    return sum(m.y[w] for w in W) <= P

pyo.SolverFactory('glpk').solve(model)

model.y.pprint()
model.x.pprint()

y : Size=3, Index={Harlingen, Memphis, Ashland}
    Key       : Lower : Value : Upper : Fixed : Stale : Domain
      Ashland :     0 :   1.0 :     1 : False : False : Binary
    Harlingen :     0 :   1.0 :     1 : False : False : Binary
      Memphis :     0 :   0.0 :     1 : False : False : Binary
x : Size=12, Index={Harlingen, Memphis, Ashland}*{NYC, LA, Chicago, Houston}
    Key                      : Lower : Value : Upper : Fixed : Stale : Domain
      ('Ashland', 'Chicago') :     0 :   1.0 :     1 : False : False :  Reals
      ('Ashland', 'Houston') :     0 :   0.0 :     1 : False : False :  Reals
           ('Ashland', 'LA') :     0 :   0.0 :     1 : False : False :  Reals
          ('Ashland', 'NYC') :     0 :   1.0 :     1 : False : False :  Reals
    ('Harlingen', 'Chicago') :     0 :   0.0 :     1 : False : False :  Reals
    ('Harlingen', 'Houston') :     0 :   1.0 :     1 : False : False :  Reals
         ('Harlingen', 'LA') :     0 :   1.0 :     1 : False : False :  Reals
        ('Harlingen', 'NYC') :     0 :   0.0 :     1 : False : False :  Reals
      ('Memphis', 'Chicago') :     0 :   0.0 :     1 : False : False :  Reals
      ('Memphis', 'Houston') :     0 :   0.0 :     1 : False : False :  Reals
           ('Memphis', 'LA') :     0 :   0.0 :     1 : False : False :  Reals
          ('Memphis', 'NYC') :     0 :   0.0 :     1 : False : False :  Reals

3.2. Changing Parameter values¶

In the tutorial slides, we saw that a parameter could be specified to be mutable. This tells Pyomo that the value of the parameter may change in the future, and allows the user to change the parameter value and resolve the problem without the need to rebuild the entire model each time. We will use this functionality to find a better solution to an earlier exercise. Considering again the knapsack problem, we would like to find when the wrench becomes valuable enough to be a part of the optimal solution. Create a Pyomo Parameter for the value of the items, make it mutable, and then write a loop that prints the solution for different wrench values. Start with the file knapsack_mutable_parameter_incomplete.py. (A solution for this problem can be found in knapsack_mutable_parameter_soln.py.)

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.item_benefit = pyo.Param(A, initialize=b, mutable=True)

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

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

opt = pyo.SolverFactory('glpk')

for wrench_benefit in range(1,11):
    model.item_benefit['wrench'] = wrench_benefit
    result_obj = opt.solve(model)

    print('Wrench benefit:', wrench_benefit, 'Objective:', pyo.value(model.obj))

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.item_benefit = pyo.Param( A, within=pyo.NonNegativeReals, initialize=b, mutable=True)

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

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

opt = pyo.SolverFactory('glpk')

for wrench_benefit in range(1,11):
    model.item_benefit['wrench'] = wrench_benefit
    result_obj = opt.solve(model)

    print('Wrench benefit:', wrench_benefit, "x['wrench']:", pyo.value(model.x['wrench']))

Wrench benefit: 1 x['wrench']: 0.0

Wrench benefit: 2 x['wrench']: 0.0

Wrench benefit: 3 x['wrench']: 0.0

Wrench benefit: 4 x['wrench']: 0.0

Wrench benefit: 5 x['wrench']: 0.0

Wrench benefit: 6 x['wrench']: 0.0

Wrench benefit: 7 x['wrench']: 0.0

Wrench benefit: 8 x['wrench']: 1.0

Wrench benefit: 9 x['wrench']: 1.0

Wrench benefit: 10 x['wrench']: 1.0

3.3. Integer cuts¶

Often, it can be important to find not only the “best” solution, but a number of solutions that are equally optimal, or close to optimal. For discrete optimization problems, this can be done using something known as an integer cut. Consider again the knapsack problem where the choice of which items to select is a discrete variable $x_i \forall i \in A$ . Let $x_i^k$ be a particular set of $x$ values we want to remove from the feasible solution space. We define an integer cut using two sets. The first set $S_0$ contains the indices for those variables whose current solution is 0, and the second set $S_1$ consists of indices for those variables whose current solution is 1. Given these two sets, an integer cut constraint that would prevent such a solution from appearing again is defined by,

\sum_{i \in S_0} x[i] + \sum_{i \in S_1} (1 - x[i]) \geq 1.

(2)

Starting with knapsack_rules.py, write a loop that solves the problem 5 times, adding an integer cut to remove the previous solution, and printing the value of the objective function and the solution at each iteration of the loop. (A solution for this problem can be found in knapsack_integer_cut_soln.py.)

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=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)

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

model.pprint()

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 )

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)

opt = pyo.SolverFactory('glpk')


# create the ConstraintList to hold the integer cuts
model.int_cuts = pyo.ConstraintList()

# loop 5 times
for l in range(5):
    # solve the problem
    result_obj = opt.solve(model)

    # print the solution
    output_str = 'Obj: ' + str(pyo.value(model.obj))
    for i in A:
        output_str += "  x[%s]: %f" % (str(i), pyo.value(model.x[i]))
    print(output_str)           
    
    # add the integer cut based on the current solution
    cut_expr = 0
    for i in A:
        if pyo.value(model.x[i]) < 0.5:
            cut_expr += model.x[i]
        else:
            cut_expr += (1.0 - model.x[i])
    model.int_cuts.add(cut_expr >= 1)

Obj: 25.0  x[hammer]: 1.000000  x[wrench]: 0.000000  x[screwdriver]: 1.000000  x[towel]: 1.000000

Obj: 20.0  x[hammer]: 0.000000  x[wrench]: 1.000000  x[screwdriver]: 1.000000  x[towel]: 1.000000

Obj: 19.0  x[hammer]: 1.000000  x[wrench]: 0.000000  x[screwdriver]: 0.000000  x[towel]: 1.000000

Obj: 17.0  x[hammer]: 0.000000  x[wrench]: 0.000000  x[screwdriver]: 1.000000  x[towel]: 1.000000

Obj: 14.0  x[hammer]: 0.000000  x[wrench]: 1.000000  x[screwdriver]: 0.000000  x[towel]: 1.000000

3.4. Lot sizing example (Hart et al., 2017)¶

We will now write a complete model from scratch using a well-known multi-period optimization problem for optimal lot-sizing adapted from Hart et al. (2017) shown below.

\begin{align} \min &~ \sum_{t \in T} c_t y_t + h_t^+ I_t^+ + h_t^- I_t^- \\ \text{s.t.} &~ I_t = I_{t-1} + X_t -d_t \quad \forall t \in T \\ &~ I_t = I_t^+ - I_t^- \quad \quad \quad ~~ \forall t \in T \\ &~ X_t \leq P y_t \quad \quad \quad \quad \quad \forall t \in T \\ &~ X_t, I_t^+, I_t^- \geq 0 \quad \quad \quad \forall t \in T \\ &~ y_t \in \{0, 1\} \quad \quad \quad \quad ~~ \forall t \in T \\ \end{align}

(3)

Our goal is to find the optimal production $X_t$ given known demands $d_t$ , fixed cost $c_t$ associated with active production in a particular time period, an inventory holding cost $h^+_t$ and a shortage cost $h_t^-$ (cost of keeping a backlog) of orders. The variable $y_t$ (binary) determines if we produce in time $t$ or not, and $I^+_t$ represents inventory that we are storing across time period $t$ , while $I^-_t$ represents the magnitude of the backlog. Note that equation $(4)$ is a constraint that only allows production in time period $t$ if the indicator variable $y_t=1$ .

Write a Pyomo model for this problem and solve it using glpk using the data provided below. You can start with the file lot_sizing_incomplete.py. (A solution is provided in lot_sizing_soln.py.)

Parameter	Description	Value
$c$	fixed cost of production	$4.6$
$I_0^+$	initial value of positive inventory	$5.0$
$I_0^-$	initial value of backlogged orders	$0.0$
$h^+$	cost (per unit) of holding inventory	$0.7$
$h^-$	shortage cost (per unit)	$1.2$
$P$	maximum production amount (big-M value)	$5$
$d$	demand	$[5, 7, 6.2, 3.1, 1.7]$

model = pyo.ConcreteModel()
model.T = pyo.RangeSet(5)    # time periods

i0 = 5.0           # initial inventory
c = 4.6            # setup cost
h_pos = 0.7        # inventory holding cost
h_neg = 1.2        # shortage cost
P = 5.0            # maximum production amount

# demand during period t
d = {1: 5.0, 2:7.0, 3:6.2, 4:3.1, 5:1.7}


model.x = pyo.Var(model.T, within=pyo.NonNegativeReals)
model.y = pyo.Var(model.T, within=pyo.Binary)
model.i_pos = pyo.Var(model.T, within=pyo.NonNegativeReals)
model.i_neg = pyo.Var(model.T, within=pyo.NonNegativeReals)

def inventory_balance(m, t):
    previous = i0 if t == m.T.first() else m.i_pos[t - 1] - m.i_neg[t - 1]
    return m.i_pos[t] - m.i_neg[t] == previous + m.x[t] - d[t]
model.inventory_balance = pyo.Constraint(model.T, rule=inventory_balance)

def production_indicator(m, t):
    return m.x[t] <= P * m.y[t]
model.production_indicator = pyo.Constraint(model.T, rule=production_indicator)

model.obj = pyo.Objective(
    expr=sum(c * model.y[t] + h_pos * model.i_pos[t] + h_neg * model.i_neg[t] for t in model.T),
    sense=pyo.minimize,
)

# solve the problem
solver = pyo.SolverFactory('glpk')
solver.solve(model)

# print the results
for t in model.T:
    print('Period: {0}, Prod. Amount: {1}'.format(t, pyo.value(model.x[t])))

model = pyo.ConcreteModel()
model.T = pyo.RangeSet(5)    # time periods

i0 = 5.0           # initial inventory
c = 4.6            # setup cost
h_pos = 0.7        # inventory holding cost
h_neg = 1.2        # shortage cost
P = 5.0            # maximum production amount

# demand during period t
d = {1: 5.0, 2:7.0, 3:6.2, 4:3.1, 5:1.7}

# define the variables
model.y = pyo.Var(model.T, domain=pyo.Binary)
model.x = pyo.Var(model.T, domain=pyo.NonNegativeReals)
model.i = pyo.Var(model.T)
model.i_pos = pyo.Var(model.T, domain=pyo.NonNegativeReals)
model.i_neg = pyo.Var(model.T, domain=pyo.NonNegativeReals)

# define the inventory relationships
def inventory_rule(m, t):
    if t == m.T.first():
        return m.i[t] == i0 + m.x[t] - d[t]
    return m.i[t] == m.i[t-1] + m.x[t] - d[t]
model.inventory = pyo.Constraint(model.T, rule=inventory_rule)

def pos_neg_rule(m, t):
    return m.i[t] == m.i_pos[t] - m.i_neg[t]
model.pos_neg = pyo.Constraint(model.T, rule=pos_neg_rule)

# create the big-M constraint for the production indicator variable
def prod_indicator_rule(m,t):
    return m.x[t] <= P*m.y[t]
model.prod_indicator = pyo.Constraint(model.T, rule=prod_indicator_rule)

# define the cost function
def obj_rule(m):
    return sum(c*m.y[t] + h_pos*m.i_pos[t] + h_neg*m.i_neg[t] for t in m.T)
model.obj = pyo.Objective(rule=obj_rule)

# solve the problem
solver = pyo.SolverFactory('glpk')
solver.solve(model)

# print the results
for t in model.T:
    print('Period: {0}, Prod. Amount: {1}'.format(t, pyo.value(model.x[t])))

Period: 1, Prod. Amount: 3.0
Period: 2, Prod. Amount: 5.0
Period: 3, Prod. Amount: 5.0
Period: 4, Prod. Amount: 5.0
Period: 5, Prod. Amount: 0.0

References¶

This notebook adapts material from the public Pyomo Summer Workshop 2018. The Pyomo modeling examples and lot-sizing discussion draw on Hart et al., 2017Bynum et al., 2021.

References¶

Hart, W. E., Laird, C. D., Watson, J.-P., Woodruff, D. L., Hackebeil, G. A., Nicholson, B. L., & Siirola, J. D. (2017). Pyomo–optimization modeling in python (Second, Vol. 67). Springer Science & Business Media.
Bynum, M. L., Hackebeil, G. A., Hart, W. E., Laird, C. D., Nicholson, B. L., Siirola, J. D., Watson, J.-P., & Woodruff, D. L. (2021). Pyomo–optimization modeling in python (Third, Vol. 67). Springer Science & Business Media.