Pyomo - Getting started

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

pip install pyomo



In [ ]:

    
from pyomo.environ import *

Definition of the problem

$$ \begin{align} \min_{x_0,x_1} & \quad -x_0 + 4 x_1 \\ \text{s.t.} & \quad -3 x_0 + x_1 \leq 6 \\ & \quad -x_0 - 2 x_1 \geq -4 \\ & \quad x_1 \geq -3 \end{align} $$

Solution

Version 1: naive implementation



In [ ]:

    
model = ConcreteModel(name="Getting started")

c = [-1, 4]
b = [ 6, 4, 3]
A = [[-3,  1],
     [ 1,  2],
     [ 0, -1]]

N = list(range(len(c)))    # num decision variables
M = list(range(len(A)))    # num constraints

model.x = Var(N)

model.obj = Objective(expr = sum(c[n] * model.x[n] for n in N) )

m = 0
model.const_0 = Constraint(expr = sum(A[m][n] * model.x[n] for n in N) <= b[m] )
m = 1
model.const_1 = Constraint(expr = sum(A[m][n] * model.x[n] for n in N) <= b[m] )
m = 2
model.const_2 = Constraint(expr = sum(A[m][n] * model.x[n] for n in N) <= b[m] )

model.pprint()

# @tail:
print()
print("-" * 60)
print()

opt = SolverFactory('glpk')

results = opt.solve(model)

model.display()

print()
print("Optimal solution: ", [value(model.x[n]) for n in N])
print("Gain of the optimal solution: ", value(model.obj))
# @:tail

Version 2: without construction rules



In [ ]:

    
model = ConcreteModel(name="Getting started")

c = [-1, 4]
b = [ 6, 4, 3]
A = [[-3,  1],
     [ 1,  2],
     [ 0, -1]]

N = list(range(len(c)))    # num decision variables
M = list(range(len(A)))    # num constraints

model.x = Var(N)

model.obj = Objective(expr = sum(c[n] * model.x[n] for n in N) )

for m in M:
    name = "const_{}".format(m)
    val = Constraint(expr = sum(A[m][n] * model.x[n] for n in N) <= b[m])
    model.add_component(name=name, val=val)

model.pprint()

# @tail:
print()
print("-" * 60)
print()

opt = SolverFactory('glpk')

results = opt.solve(model)

model.display()

print()
print("Optimal solution: ", [value(model.x[n]) for n in N])
print("Gain of the optimal solution: ", value(model.obj))
# @:tail

Version 3: with construction rules



In [ ]:

    
model = ConcreteModel(name="Getting started")

c = [-1, 4]
b = [ 6, 4, 3]
A = [[-3,  1],
     [ 1,  2],
     [ 0, -1]]

N = list(range(len(c)))    # num decision variables
M = list(range(len(A)))    # num constraints

model.x = Var(N)

model.obj = Objective(expr = sum(c[n] * model.x[n] for n in N) )

def constraint_fn(model, m):
    return sum(A[m][n] * model.x[n] for n in N) <= b[m]

model.constraint = Constraint(M, rule=constraint_fn)

model.pprint()

# @tail:
print()
print("-" * 60)
print()

opt = SolverFactory('glpk')

results = opt.solve(model)

model.display()

print()
print("Optimal solution: ", [value(model.x[n]) for n in N])
print("Gain of the optimal solution: ", value(model.obj))
# @:tail

Version 4: construction rules on Objective too



In [ ]:

    
model = ConcreteModel(name="Getting started")

c = [-1, 4]
b = [ 6, 4, 3]
A = [[-3,  1],
     [ 1,  2],
     [ 0, -1]]

N = list(range(len(c)))    # num decision variables
M = list(range(len(A)))    # num constraints

model.x = Var(N)


def objective_fn(model):
    return sum(c[n] * model.x[n] for n in N)

model.obj = Objective(rule=objective_fn)


def constraint_fn(model, m):
    return sum(A[m][n] * model.x[n] for n in N) <= b[m]

model.constraint = Constraint(M, rule=constraint_fn)


model.pprint()

# @tail:
print()
print("-" * 60)
print()

opt = SolverFactory('glpk')

results = opt.solve(model)

model.display()

print()
print("Optimal solution: ", [value(model.x[n]) for n in N])
print("Gain of the optimal solution: ", value(model.obj))
# @:tail