Solving problems using Python¶

Colab notebook

We will solve LP problems using Python. In case you are less familiar with Python, I recommend you to go over the Python Review section.

In this course we will mostly use Gurobi to solve the optimization problems. The Gurobi optimizer is a mathematical optimization software library for solving mixed-integer linear and quadratic optimization problems.

For more information, see: https://www.gurobi.com/academia/academic-program-and-licenses/

Other useful resources can be found here

# install the Gurobi package
!pip install gurobipy

Requirement already satisfied: gurobipy in /usr/local/lib/python3.7/dist-packages (9.5.0)

# Import gurobi library
from gurobipy import * # This command imports the Gurobi functions and classes.

# Create new model
m = Model('') # The Model() constructor creates a model object m. The name of this new model is 'Factory'.
                     # This new model m initially contains no decision variables, constraints, or objective function.

Restricted license - for non-production use only - expires 2023-10-25

# Create decision variables
# This method adds a decision variable to the model object m, one by one; i.e. x1 and then x2. 
# The argument of the method gives the name of added decision variable. 
# The default values are applied here; i.e. the decision variables are of type continuous and non-negative, with no upper bound.
x1 = m.addVar(lb=0, vtype = GRB.CONTINUOUS, name='chairs') 
x2 = m.addVar(lb=0, vtype = GRB.CONTINUOUS, name='tables')

# Define objective function
# This method adds the objective function to the model object m. The first
# argument is a linear expression and the second argument defines
# the sense of the optimization.
m.setObjective(40*x1+50*x2, GRB.MAXIMIZE)

# Add constraints
# This method adds a constraint to the model object m
m.addConstr(1*x1+2*x2<=40, name='wood')
m.addConstr(4*x1+3*x2<=120, name= 'labor')

<gurobi.Constr *Awaiting Model Update*>

# Run optimization engine
# This method runs the optimization engine to solve the LP problem in the model object m
m.optimize()

Gurobi Optimizer version 9.5.0 build v9.5.0rc5 (linux64)
Thread count: 1 physical cores, 2 logical processors, using up to 2 threads
Optimize a model with 2 rows, 2 columns and 4 nonzeros
Model fingerprint: 0x3a526911
Coefficient statistics:
  Matrix range     [1e+00, 4e+00]
  Objective range  [4e+01, 5e+01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [4e+01, 1e+02]
Presolve time: 0.01s
Presolved: 2 rows, 2 columns, 4 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    9.0000000e+31   3.250000e+30   9.000000e+01      0s
       2    1.3600000e+03   0.000000e+00   0.000000e+00      0s

Solved in 2 iterations and 0.02 seconds (0.00 work units)
Optimal objective  1.360000000e+03

#display optimal production plan
for v in m.getVars():
  print(v.varName, v.x)
print('optimal total revenue:', m.objVal)

chairs 24.0
tables 8.0
optimal total revenue: 1360.0

Given the constraints, the maximum revenue is $1360 structured such that we produce 24 chairs and 8 tables.

What if our LP problem has hundreds of thousands varibales and constraints?

The Gurobi python code just presented is too manual and would take too long too build a large scale LP problem. We should use appropriate data structures, python functions and objects to abstract the problem to any size.

Python list comprehension:

#List comprehension is compact way to create lists
sqrd = [i*i for i in range(10)]
print(sqrd) 

#Can be used to create subsequences that satisfy certain conditions (ex: filtering a list)
bigsqrd = [i*i for i in range(10) if i*i >= 5]
print(bigsqrd) 

#Can be used with multiple for loops (ex: all combinations)
prod = [i+j for i in range(3) for j in range(3)]
print(prod) 

#Generator expression is similar, but no brackets (ex: argument to aggregate sum)
sum_sq = sum(i for i in range(11))
print(sum_sq)

[0, 1, 4, 9, 16, 25, 36, 49, 64, 81]
[9, 16, 25, 36, 49, 64, 81]
[0, 1, 2, 1, 2, 3, 2, 3, 4]
55

# Resource data
# The multidict function returns a list which maps each resource (key) to its capacity value.
resources, capacity = multidict({ 
    'wood':  40,
    'labor': 120 })
print(resources, capacity)

['wood', 'labor'] {'wood': 40, 'labor': 120}

# Products data
# This multidict function returns a list which maps each product (key) to its price value.
products, price = multidict({
    'chair': 40,
    'table': 50 })
print(products, price)

['chair', 'table'] {'chair': 40, 'table': 50}

# Bill of materials: resources required by each product
# This dictionary has a 2-tuple as a key, mapping the resource required by a product with its quantity per.
bom={
('wood','chair'):1,
('wood','table'):2,
('labor','chair'):4,
('labor','table'):3
}
print(bom)

{('wood', 'chair'): 1, ('wood', 'table'): 2, ('labor', 'chair'): 4, ('labor', 'table'): 3}

m = Model('Factory')

# This method adds decision variables to the model object m
make = m.addVars(products, name='make')

# This method adds constraints to the model object m
res = m.addConstrs(((sum(bom[r,p]*make[p] for p in products) <= capacity[r]) for r in resources), name='R')

# This method adds the objective function to the model object m.
# The first argument is a linear expression which is generated by the 'prod' method. 
# The 'prod' method is the product of the object (revenue) with the object (make) 
# for each product p in the set (products). The second argument defines the sense of the optimization.
m.setObjective(make.prod(price), GRB.MAXIMIZE)

#save model for inspection
m.write('factory.lp')

cat factory.lp

\ Model Factory
\ LP format - for model browsing. Use MPS format to capture full model detail.
Maximize
  40 make[chair] + 50 make[table]
Subject To
 R[wood]: make[chair] + 2 make[table] <= 40
 R[labor]: 4 make[chair] + 3 make[table] <= 120
Bounds
End

# run optimization engine
m.optimize()

Gurobi Optimizer version 9.5.0 build v9.5.0rc5 (linux64)
Thread count: 1 physical cores, 2 logical processors, using up to 2 threads
Optimize a model with 2 rows, 2 columns and 4 nonzeros
Model fingerprint: 0x3a526911
Coefficient statistics:
  Matrix range     [1e+00, 4e+00]
  Objective range  [4e+01, 5e+01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [4e+01, 1e+02]
Presolve time: 0.01s
Presolved: 2 rows, 2 columns, 4 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    9.0000000e+31   3.250000e+30   9.000000e+01      0s
       2    1.3600000e+03   0.000000e+00   0.000000e+00      0s

Solved in 2 iterations and 0.03 seconds (0.00 work units)
Optimal objective  1.360000000e+03

#display optimal production plan
for v in m.getVars():
  print(v.varName, v.x)
print('optimal total revenue:', m.objVal)

make[chair] 24.0
make[table] 8.0
optimal total revenue: 1360.0

OR/MS in Python

Solving problems using Python¶

Sensitivity analysis of LP problems¶

Objective value coefficients¶

Constraint quantity values¶

Shadow prices¶

Sensitivity of the shadow price¶

Adding new variable and/or new constraint¶

Problem example: C2Q1¶

Problem example: C2Q2¶

Problem example: C4Q8¶

Problem example: C4Q32¶

Problem example: C4Q33¶