Pyomo - Getting started

Pyomo installation: see http://www.pyomo.org/installation

pip install pyomo



In [ ]:

    
from pyomo.environ import *

Ice cream example

This example is taken from the following book: Pyomo - Optimization Modeling in Python by W. E. Hart & al. , Second Edition, Springer (p.19)

$$ \begin{align} \max_{x} & \quad \sum_{i \in \mathcal{A}} h_i (1 - u/d_i^2) x_i \\ \text{s.t.} & \quad \sum_{i \in \mathcal{A}} c_i x_i \leq b \\ & \quad 0 \leq x_i \leq u_i, \quad i \in \mathcal{A} \end{align} $$

Concrete Model

ConcHLinScript.py in https://github.com/Pyomo/pyomo/tree/master/examples/doc/pyomobook/optimization-ch



In [ ]:

    
instance = ConcreteModel(name="Linear (H)")

A = ['I_C_Scoops', 'Peanuts']
h = {'I_C_Scoops': 1,    'Peanuts': 0.1}
d = {'I_C_Scoops': 5,    'Peanuts': 27}
c = {'I_C_Scoops': 3.14, 'Peanuts': 0.2718}
b = 12
u = {'I_C_Scoops': 100,  'Peanuts': 40.6}

def x_bounds(m, i):
    return (0,u[i])

instance.x = Var(A, bounds=x_bounds)

def obj_rule(instance):
    return sum(h[i]*(1 - u[i]/d[i]**2) * instance.x[i] for i in A)

instance.z = Objective(rule=obj_rule, sense=maximize)

instance.budgetconstr = Constraint(expr = sum(c[i] * instance.x[i] for i in A) <= b)

# @tail:
opt = SolverFactory('glpk')

results = opt.solve(instance)    # solves and updates instance

instance.display()
# @:tail

Abstract Model

AbstHLinScript.py in https://github.com/Pyomo/pyomo/tree/master/examples/doc/pyomobook/optimization-ch



In [ ]:

    
DATA_STR = """# Pyomo data file for AbstractH.py
set A := I_C_Scoops Peanuts ;
param h := I_C_Scoops 1 Peanuts 0.1 ;
param d := 
  I_C_Scoops 5
  Peanuts 27 ;
param c := I_C_Scoops 3.14 Peanuts 0.2718 ;
param b := 12 ;
param u := I_C_Scoops 100 Peanuts 40.6 ;
"""

with open("AbstractH.dat", "w") as fd:
    print(DATA_STR, file=fd)



In [ ]:

    
!cat AbstractH.dat



In [ ]:

    
model = AbstractModel(name="Simple Linear (H)")

model.A = Set()

model.h = Param(model.A)
model.d = Param(model.A)
model.c = Param(model.A)
model.b = Param()
model.u = Param(model.A)

def xbounds_rule(model, i):
    return (0, model.u[i])
model.x = Var(model.A, bounds=xbounds_rule)

def obj_rule(model):
    return sum(model.h[i] * (1 - model.u[i]/model.d[i]**2) * model.x[i] for i in model.A)

model.z = Objective(rule=obj_rule, sense=maximize)

def budget_rule(model):
    return summation(model.c, model.x) <= model.b

model.budgetconstr = Constraint(rule=budget_rule)

# @tail:
opt = SolverFactory('glpk')

instance = model.create_instance("AbstractH.dat")
results = opt.solve(instance)     # solves and updates instance

instance.display()
# @:tail