Sören Petrat, Operations Research, Fall 2021

This file is a generalization of the basic "WYNDOR GLASS CO." example file which shows how to use more abstract models. (The code is based on previous examples by Marcel Oliver.)

from pyomo.environ import *
from pyomo.opt import *
opt = solvers.SolverFactory("glpk")

Put model data directly in the code (as in previous example)

model = ConcreteModel()

model.x = Var([1,2], within=NonNegativeReals)

model.c1 = Constraint(expr = model.x[1] <= 4)
model.c2 = Constraint(expr = 2*model.x[2] <= 12)
model.c3 = Constraint(expr = 3*model.x[1] + 2*model.x[2] <= 18)
model.z = Objective(expr = 3*model.x[1] + 5*model.x[2], sense=maximize)

results = opt.solve(model)

model.x.get_values()

{1: 2.0, 2: 6.0}

model.z.expr()

36.0

Better: Define data separately beforehand

model2 = ConcreteModel()

A = [1,2] #set of activities
R = ['Plant1','Plant2','Plant3']  #set of resources/production facilities
p = {1:3, 2:5}  #profit
w = {(1,'Plant1'):1, (2,'Plant1'):0,
     (1,'Plant2'):0, (2,'Plant2'):2,
     (1,'Plant3'):3, (2,'Plant3'):2}  #workloads
c = {'Plant1':4, 'Plant2':12, 'Plant3':18}  #capactities

model2.x = Var(A, within=NonNegativeReals)

#Maximize Profit
def profit_rule(model2):
    return sum(p[i]*model2.x[i] for i in A)
model2.profit = Objective(rule=profit_rule, sense=maximize)

#Constraints
def capacity_rule(model2,j):
    return sum(w[i,j]*model2.x[i] for i in A) <= c[j]
model2.con = Constraint(R,rule=capacity_rule)

results2 = opt.solve(model2)

model2.x.get_values()

{1: 2.0, 2: 6.0}

model2.profit.expr()

36.0

Even better: Define an abstract model first and then fill it with data from a file

model3 = AbstractModel()

#Products
model3.P = Set()
#Resources, Machines, etc
model3.M = Set()

#Available time, capacities etc
model3.t3 = Param(model3.M, within=PositiveReals)
# Profit for each product
model3.p3 = Param(model3.P, within=PositiveReals)
# Work capacity requirement for each machine and each product
model3.w3 = Param(model3.M, model3.P, within=NonNegativeReals)

# Number of units to be produced for each product
model3.x3 = Var(model3.P, within=NonNegativeReals)

#Profit
def profit_rule3(model3):
    return sum(model3.p3[i]*model3.x3[i] for i in model3.P)
model3.profit3 = Objective(rule=profit_rule3, sense=maximize)

#Constraints
def capacity_rule3(model3, j):
    return sum(model3.w3[j,i]*model3.x3[i] for i in model3.P) <= model3.t3[j]
model3.con3 = Constraint(model3.M, rule=capacity_rule3)
#model3.pprint()

instance = model3.create_instance('production_data.dat') #creates a specific instance of the abstract model
#instance.pprint()

results3 = opt.solve(instance) #solve the instance (not the model)

instance.x3.get_values()

{'Glass_Door': 2.0, 'Wood_Window': 6.0}

instance.profit3.expr()

36.0

instance.t3.extract_values()

{'Aluminum': 4, 'Wood': 12, 'Assembly': 18}

instance.p3.extract_values()

{'Glass_Door': 3, 'Wood_Window': 5}

instance.w3.extract_values()

{('Aluminum', 'Glass_Door'): 1,
 ('Wood', 'Glass_Door'): 0,
 ('Assembly', 'Glass_Door'): 3,
 ('Aluminum', 'Wood_Window'): 0,
 ('Wood', 'Wood_Window'): 2,
 ('Assembly', 'Wood_Window'): 2}