Pyomo solving Farmer's problem step by step improvement

farmer_lp.py

# -*- coding: utf-8
"""
A deterministic model of farmer's problem:

"""
try:
    from pyomo.core import *
    from pyomo import *
    from pyomo.opt import *
    from pyomo.core.base import *
    from pyomo.environ import *
except:
    # from __future__ import division
    from coopr.pyomo import *
    from coopr.opt import SolverStatus, TerminationCondition, SolutionStatus
    from coopr.opt.base import SolverFactory
    import coopr.environ

model = ConcreteModel()

model.CROPS = Set(initialize=["W", "C", "S"])
model.TOTAL_ACREAGE = 500

""" first stage variables """
model.plant = Var(model.CROPS, within=NonNegativeReals)
model.sell = Var(model.CROPS, within=NonNegativeReals)
model.buy = Var(model.CROPS, within=NonNegativeReals)
model.sell_extra = Var(within=NonNegativeReals)
""" first stage constraints """
model.con_total_arc = Constraint(expr=summation(model.plant) <= model.TOTAL_ACREAGE)
model.con_feed_cattle_W = Constraint(expr=model.plant['W'] * 2.5 - model.sell['W'] + model.buy['W'] >= 200)
model.con_feed_cattle_C = Constraint(expr=model.plant['C'] * 3 - model.sell['C'] + model.buy['C'] >= 240)
model.con_sell_W = Constraint(expr=model.plant['W'] * 2.5 + model.buy['W'] >= model.sell["W"])
model.con_sell_C = Constraint(expr=model.plant['C'] * 3 + model.buy['C'] >= model.sell["C"])
model.con_sell_S_extra = Constraint(expr=model.plant['S'] * 20 >= model.sell["S"] + model.sell_extra)
model.con_sell_S = Constraint(expr=model.sell['S'] <= 6000)


def total_cost_rule(model):
    total_sell_S = 36 * model.sell['S'] + 10 * model.sell_extra
    total_sell_W = model.sell["W"] * 170
    total_sell_C = model.sell["C"] * 150
    total_buy_W = model.buy["W"] * 238
    total_buy_C = model.buy["C"] * 210
    total_plant = model.plant["W"] * 150 + model.plant["C"] * 230 + model.plant["S"] * 260
    return total_plant + total_buy_W + total_buy_C - total_sell_S - total_sell_W - total_sell_C
model.total_cost = Objective(rule=total_cost_rule, sense=minimize)

instance = model.create()
opt = SolverFactory("cbc")
results = opt.solve(instance)
instance.load(results)
print results

farmer_sp.py

# -*- coding: utf-8
"""
A explicit implementation of two-stage SP of farmer's problem
Uncertainty:
* Yield rate for each plant:
    H: 1.2X
    M: 1X
    L: 0.8X

The fisrt stage is the arceage of each plant
The second stage is the amount to buy and sell of each plant
"""
try:
    from pyomo.core import *
    from pyomo import *
    from pyomo.opt import *
    from pyomo.core.base import *
    from pyomo.environ import *
except:
    # from __future__ import division
    from coopr.pyomo import *
    from coopr.opt import SolverStatus, TerminationCondition, SolutionStatus
    from coopr.opt.base import SolverFactory
    import coopr.environ

model = ConcreteModel()

model.CROPS = Set(initialize=["W", "C", "S"])
model.TOTAL_ACREAGE = 500

""" first stage variables """
model.plant = Var(model.CROPS, within=NonNegativeReals)
model.sell = Var(model.CROPS, within=NonNegativeReals)
model.buy = Var(model.CROPS, within=NonNegativeReals)
model.sell_extra = Var(within=NonNegativeReals)

""" first stage constraints """
model.con_total_arc = Constraint(expr=summation(model.plant) <= model.TOTAL_ACREAGE)
model.con_feed_cattle_W = Constraint(expr=model.plant['W'] * 2.5 - model.sell['W'] + model.buy['W'] >= 200)
model.con_feed_cattle_C = Constraint(expr=model.plant['C'] * 3 - model.sell['C'] + model.buy['C'] >= 240)
model.con_sell_W = Constraint(expr=model.plant['W'] * 2.5 + model.buy['W'] >= model.sell["W"])
model.con_sell_C = Constraint(expr=model.plant['C'] * 3 + model.buy['C'] >= model.sell["C"])
model.con_sell_S_extra = Constraint(expr=model.plant['S'] * 20 >= model.sell["S"] + model.sell_extra)
model.con_sell_S = Constraint(expr=model.sell['S'] <= 6000)
# model.con_sell_S = Constraint(expr=model.plant['S']*20 >= model.sell["S"])

model.scenarios = Set(initialize=["H", "M", "L"])  # high and low
""" second stage variables """
model.sell_H = Var(model.CROPS, within=NonNegativeReals)
model.buy_H = Var(model.CROPS, within=NonNegativeReals)
model.sell_extra_H = Var(within=NonNegativeReals)

model.sell_M = Var(model.CROPS, within=NonNegativeReals)
model.buy_M = Var(model.CROPS, within=NonNegativeReals)
model.sell_extra_M = Var(within=NonNegativeReals)

model.sell_L = Var(model.CROPS, within=NonNegativeReals)
model.buy_L = Var(model.CROPS, within=NonNegativeReals)
model.sell_extra_L = Var(within=NonNegativeReals)

""" second stage constraints """
model.con_total_arc_H = Constraint(expr=summation(model.plant) <= model.TOTAL_ACREAGE)
model.con_feed_cattle_W_H = Constraint(expr=model.plant['W'] * 2.5 * 1.2 - model.sell_H['W'] + model.buy_H['W'] >= 200)
model.con_feed_cattle_C_H = Constraint(expr=model.plant['C'] * 3 * 1.2 - model.sell_H['C'] + model.buy_H['C'] >= 240)
model.con_sell_W_H = Constraint(expr=model.plant['W'] * 2.5 * 1.2 + model.buy_H['W'] >= model.sell_H["W"])
model.con_sell_C_H = Constraint(expr=model.plant['C'] * 3 * 1.2 + model.buy_H['C'] >= model.sell_H["C"])
model.con_sell_S_extra_H = Constraint(expr=model.plant['S'] * 20 * 1.2 >= model.sell_H["S"] + model.sell_extra_H)
model.con_sell_S_H = Constraint(expr=model.sell_H['S'] <= 6000)

model.con_total_arc_M = Constraint(expr=summation(model.plant) <= model.TOTAL_ACREAGE)
model.con_feed_cattle_W_M = Constraint(expr=model.plant['W'] * 2.5 - model.sell_M['W'] + model.buy_M['W'] >= 200)
model.con_feed_cattle_C_M = Constraint(expr=model.plant['C'] * 3 - model.sell_M['C'] + model.buy_M['C'] >= 240)
model.con_sell_W_M = Constraint(expr=model.plant['W'] * 2.5 + model.buy_M['W'] >= model.sell_M["W"])
model.con_sell_C_M = Constraint(expr=model.plant['C'] * 3 + model.buy_M['C'] >= model.sell_M["C"])
model.con_sell_S_extra_M = Constraint(expr=model.plant['S'] * 20 >= model.sell_M["S"] + model.sell_extra_M)
model.con_sell_S_M = Constraint(expr=model.sell_M['S'] <= 6000)

model.con_total_arc_L = Constraint(expr=summation(model.plant) <= model.TOTAL_ACREAGE)
model.con_feed_cattle_W_L = Constraint(expr=model.plant['W'] * 2.5 * 0.8 - model.sell_L['W'] + model.buy_L['W'] >= 200)
model.con_feed_cattle_C_L = Constraint(expr=model.plant['C'] * 3 * 0.8 - model.sell_L['C'] + model.buy_L['C'] >= 240)
model.con_sell_W_L = Constraint(expr=model.plant['W'] * 2.5 * 0.8 + model.buy_L['W'] >= model.sell_L["W"])
model.con_sell_C_L = Constraint(expr=model.plant['C'] * 3 * 0.8 + model.buy_L['C'] >= model.sell_L["C"])
model.con_sell_S_extra_L = Constraint(expr=model.plant['S'] * 20 * 0.8 >= model.sell_L["S"] + model.sell_extra_L)
model.con_sell_S_L = Constraint(expr=model.sell_L['S'] <= 6000)


def first_stage_cost(model):
    total_plant = model.plant["W"] * 150 + model.plant["C"] * 230 + model.plant["S"] * 260
    return total_plant
model.first_stage_cost = Expression(expr=first_stage_cost)


def second_stage_cost(model):
    total_sell_S_H = 36 * model.sell_H['S'] + 10 * model.sell_extra_H
    total_sell_W_H = model.sell_H["W"] * 170
    total_sell_C_H = model.sell_H["C"] * 150
    total_buy_W_H = model.buy_H["W"] * 238
    total_buy_C_H = model.buy_H["C"] * 210
    total_H = total_buy_W_H + total_buy_C_H - total_sell_S_H - total_sell_W_H - total_sell_C_H

    total_sell_S_M = 36 * model.sell_M['S'] + 10 * model.sell_extra_M
    total_sell_W_M = model.sell_M["W"] * 170
    total_sell_C_M = model.sell_M["C"] * 150
    total_buy_W_M = model.buy_M["W"] * 238
    total_buy_C_M = model.buy_M["C"] * 210
    total_M = total_buy_W_M + total_buy_C_M - total_sell_S_M - total_sell_W_M - total_sell_C_M

    total_sell_S_L = 36 * model.sell_L['S'] + 10 * model.sell_extra_L
    total_sell_W_L = model.sell_L["W"] * 170
    total_sell_C_L = model.sell_L["C"] * 150
    total_buy_W_L = model.buy_L["W"] * 238
    total_buy_C_L = model.buy_L["C"] * 210
    total_L = total_buy_W_L + total_buy_C_L - total_sell_S_L - total_sell_W_L - total_sell_C_L

    return 1. / 3 * total_H + 1. / 3 * total_M + 1. / 3 * total_L
model.second_stage_cost = Expression(expr=second_stage_cost)


def total_cost_rule(model):
    return model.first_stage_cost + model.second_stage_cost
model.total_expected_cost = Objective(rule=total_cost_rule, sense=minimize)

instance = model.create()
opt = SolverFactory("cbc", tee=True)
results = opt.solve(instance)
instance.load(results)
print results

farmer_sp_2stage.py

# -*- coding: utf-8
"""
A implicit implementation of two-stage SP of farmer's problem
Uncertainty:
* Yield rate for each plant:
    H: 1.2X
    M: 1X
    L: 0.8X

The fisrt stage is the arceage of each plant
The second stage is the amount to buy and sell of each plant
"""
try:
    from pyomo.core import *
    from pyomo import *
    from pyomo.opt import *
    from pyomo.core.base import *
    from pyomo.environ import *
except:
    # from __future__ import division
    from coopr.pyomo import *
    from coopr.opt import SolverStatus, TerminationCondition, SolutionStatus
    from coopr.opt.base import SolverFactory
    import coopr.environ
    # from pyomo.environ import *

model = ConcreteModel()

model.CROPS = Set(initialize=["W", "C", "S"])
model.SCENARIOS = Set(initialize=["H", "M", "L"])  # high and low
model.TOTAL_ACREAGE = 500
model.PROB = Param(model.SCENARIOS, initialize={"H": 1./3, "M": 1./3, "L": 1./3})
""" first stage variables """
model.plant = Var(model.CROPS, within=NonNegativeReals)
model.yield_rate = Var(model.CROPS, initialize={"W": 2.5, "C": 3, "S": 20})

""" first stage constraints """


def con_total_arc(model):
    return summation(model.plant) <= model.TOTAL_ACREAGE
model.con_total_arc = Constraint(rule=con_total_arc)

""" second stage variables """
model.buy_2 = Var(model.CROPS, model.SCENARIOS, within=NonNegativeReals)
model.sell_2 = Var(model.CROPS, model.SCENARIOS, within=NonNegativeReals)
model.sell_extra_2 = Var(model.SCENARIOS, within=NonNegativeReals)
model.total_2 = Var(model.SCENARIOS)

def yield_rate_init(model, crops, scenarios):
    if scenarios == "H":
        return model.yield_rate[crops] * 1.2
    elif scenarios == "M":
        return model.yield_rate[crops] * 1
    elif scenarios == "L":
        return model.yield_rate[crops] * 0.8
model.yield_rate_2 = Param(model.CROPS, model.SCENARIOS, initialize=yield_rate_init)

""" second stage constraints """


def con_feed_cattle_W(model, s):
    return model.plant["W"] * model.yield_rate_2["W", s] - model.sell_2["W", s] + model.buy_2["W", s] >= 200
model.con_feed_cattle_W = Constraint(model.SCENARIOS, rule=con_feed_cattle_W)


def con_feed_cattle_C(model, s):
    return model.plant["C"] * model.yield_rate_2["C", s] - model.sell_2["C", s] + model.buy_2["C", s] >= 240
model.con_feed_cattle_C = Constraint(model.SCENARIOS, rule=con_feed_cattle_C)


def con_sell_W(model, s):
    return model.plant['W'] * model.yield_rate_2["W", s] + model.buy_2['W', s] >= model.sell_2["W", s]
model.con_sell_W = Constraint(model.SCENARIOS, rule=con_sell_W)


def con_sell_C(model, s):
    return model.plant['C'] * model.yield_rate_2["C", s] + model.buy_2['C', s] >= model.sell_2["C", s]
model.con_sell_C = Constraint(model.SCENARIOS, rule=con_sell_C)


def con_sell_S_extra(model, s):
    return model.plant['S'] * model.yield_rate_2["S", s] >= model.sell_2["S", s] + model.sell_extra_2[s]
model.con_sell_S_extra = Constraint(model.SCENARIOS, rule=con_sell_S_extra)


def con_sell_S(model, s):
    return model.sell_2["S", s] <= 6000
model.con_sell_S = Constraint(model.SCENARIOS, rule=con_sell_S)


def con_second_stage_cost(model, s):
    return model.total_2[s] == model.PROB[s] * ((model.buy_2["W", s] * 238) + (model.buy_2["C", s] * 210) - (model.sell_2["C", s] * 150) - (model.sell_2["W", s] * 170) - (36 * model.sell_2['S', s] + 10 * model.sell_extra_2[s]))
model.con_second_stage_cost = Constraint(model.SCENARIOS, rule=con_second_stage_cost)

""" Cost calculations """


def first_stage_cost(model):
    total_plant = model.plant["W"] * 150 + model.plant["C"] * 230 + model.plant["S"] * 260
    return total_plant
model.first_stage_cost = Expression(expr=first_stage_cost)


def second_stage_cost(model):
    return summation(model.total_2)
model.second_stage_cost = Expression(expr=second_stage_cost)


def total_cost_rule(model):
    return model.first_stage_cost + model.second_stage_cost
model.total_expected_cost = Objective(rule=total_cost_rule, sense=minimize)

instance = model.create()
opt = SolverFactory("cbc")
results = opt.solve(instance)
instance.load(results)
# results.write()
# print results
# print instance
# print dict(instance.plant)
plant_r = {}
for key in instance.plant.keys():
    plant_r[key] = instance.plant[key].value
print plant_r
print instance.__dict__.keys()
print value(instance.first_stage_cost)
print value(instance.second_stage_cost)
print value(instance.total_expected_cost)
print [x() for x in instance.total_2.values()]