Pyomo Fundamentals

This notebook contains material from the Pyomo Summer Workshop 2018 by the Pyomo Developement Team and the notebook is developed by David Bernal (dbernaln@purdue.edu), Zedong Peng (peng372@purdue.edu), and Albert Lee (lee4382@purdue.edu); the content is available on Github. The original sources of this material are available on Pyomo.

import pyomo.environ as pyo
import pandas as pd

Exercise 1: Pyomo Fundamentals

1.1. Knapsack example

Solve the knapsack problem shown in tutorial using your IDE (e.g., Spyder) or the command lined: > python knapsack.py. 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=x_index
        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 )
# TODO: INSERT CODE HERE TO PRINT TOTAL WEIGHT AND BENEFIT

print('%12s %12s' % ('Item', 'Selected'))
print('=========================')
for i in A:
    # TODO: INSERT CODE HERE TO PRINT EACH ITEM AND WHETHER OR NOT IT WAS SELECTED
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 = # TODO: WRITE CODE TO GET A DICTIONARY FROM THE DATAFRAME
w = # TODO: WRITE CODE TO GET A DICTIONARY FROM THE DATAFRAME

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 Set Declarations
    x_index : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    4 : {'hammer', 'wrench', 'screwdriver', 'towel'}

1 Var Declarations
    x : Size=4, Index=x_index
        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

4 Declarations: x_index x obj weight_con

Exercise 2: 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. 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 Set Declarations
    x_index : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    4 : {'hammer', 'wrench', 'screwdriver', 'towel'}

1 Var Declarations
    x : Size=4, Index=x_index
        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

4 Declarations: x_index 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{split} \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} \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} \end{split}\]

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

Exercies 3: Pyomo Fundamentals

3.1. Using the 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

# TODO: ADD DECORATOR HERE
def warehouse_active(m, w, c):
    return m.x[w,c] <= m.y[w]

# TODO: ADD DECORATOR HERE
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=y_index
    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=x_index
    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 = # TODO: DEFINE THE PYOMO PARAM HERE

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)

    # TODO: PRINT THE BENEFIT OF THE WRENCH FOR EACH ITERATION IN THE LOOP

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. \]

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. Putting it all together with the 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 Hagen et al. (2001) shown below.

\[\begin{split} \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} \end{split}\]

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}


# TODO: WRITE CODE FOR THE MODEL HERE

# 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

1. Hart, W. E., Laird, C. D., Watson, J. P., Woodrff, D. L., Hackebeil, G. A., Nicholson, B. L., and Siirola, J. D. Pyomo: Optimization Modeling in Python (Second Edition), Vol (67), Springer Verlag, 2017.