An example of using Coopr (Pymo) and PySP with ConcreteModel

farmer_concrete.py

#  _________________________________________________________________________
#
#  Pyomo: Python Optimization Modeling Objects
#  Copyright (c) 2014 Sandia Corporation.
#  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
#  the U.S. Government retains certain rights in this software.
#  This software is distributed under the BSD License.
#  _________________________________________________________________________
#
# Farmer: rent out version has a scalar root node var
# note: this will minimize
#
# Imports
#
""" Modify the example to the coopr version with 3.5.8787 """
from coopr.core import *
from coopr.pyomo import *
from coopr.opt import SolverStatus, TerminationCondition, SolutionStatus
from coopr.opt.base import SolverFactory
import coopr.environ

#
# Model
#

model = ConcreteModel()

#
# Parameters
#

model.CROPS = Set(initialize=['WHEAT', 'CORN', 'SUGAR_BEETS'])

model.TOTAL_ACREAGE = 500.0

model.PriceQuota = {'WHEAT': 1000.0, 'CORN': 1000.0, 'SUGAR_BEETS': 6000.0}
model.SubQuotaSellingPrice = {'WHEAT': 170.0, 'CORN': 150.0, 'SUGAR_BEETS': 36.0}
model.SuperQuotaSellingPrice = {'WHEAT': 0.0, 'CORN': 0.0, 'SUGAR_BEETS': 10.0}
model.CattleFeedRequirement = {'WHEAT': 200.0, 'CORN': 240.0, 'SUGAR_BEETS': 0.0}
model.PurchasePrice = {'WHEAT': 238.0, 'CORN': 210.0, 'SUGAR_BEETS': 100000.0}
model.PlantingCostPerAcre = {'WHEAT': 150.0, 'CORN': 230.0, 'SUGAR_BEETS': 260.0}
model.Yield = Param(model.CROPS, within=NonNegativeReals, initialize=0.0, mutable=True)

#
# Variables
#

model.DevotedAcreage = Var(model.CROPS, bounds=(0.0, model.TOTAL_ACREAGE))
model.QuantitySubQuotaSold = Var(model.CROPS, bounds=(0.0, None))
model.QuantitySuperQuotaSold = Var(model.CROPS, bounds=(0.0, None))
model.QuantityPurchased = Var(model.CROPS, bounds=(0.0, None))

#
# Constraints
#


def ConstrainTotalAcreage_rule(model):
    return summation(model.DevotedAcreage) <= model.TOTAL_ACREAGE
model.ConstrainTotalAcreage = Constraint(rule=ConstrainTotalAcreage_rule)


def EnforceCattleFeedRequirement_rule(model, i):
    return model.CattleFeedRequirement[i] <= (model.Yield[i] * model.DevotedAcreage[i]) + model.QuantityPurchased[i] - model.QuantitySubQuotaSold[i] - model.QuantitySuperQuotaSold[i]
model.EnforceCattleFeedRequirement = Constraint(model.CROPS, rule=EnforceCattleFeedRequirement_rule)


def LimitAmountSold_rule(model, i):
    return model.QuantitySubQuotaSold[i] + model.QuantitySuperQuotaSold[i] - (model.Yield[i] * model.DevotedAcreage[i]) <= 0.0
model.LimitAmountSold = Constraint(model.CROPS, rule=LimitAmountSold_rule)


def EnforceQuotas_rule(model, i):
    return (0.0, model.QuantitySubQuotaSold[i], model.PriceQuota[i])
model.EnforceQuotas = Constraint(model.CROPS, rule=EnforceQuotas_rule)

#
# Stage-specific cost computations
#


def ComputeFirstStageCost_rule(model):
    return summation(model.PlantingCostPerAcre, model.DevotedAcreage)
model.FirstStageCost = Expression(expr=ComputeFirstStageCost_rule)


def ComputeSecondStageCost_rule(model):
    expr = summation(model.PurchasePrice, model.QuantityPurchased)
    expr -= summation(model.SubQuotaSellingPrice, model.QuantitySubQuotaSold)
    expr -= summation(model.SuperQuotaSellingPrice, model.QuantitySuperQuotaSold)
    return expr
model.SecondStageCost = Expression(expr=ComputeSecondStageCost_rule)

#
# PySP Auto-generated Objective
#
# minimize: sum of StageCostVariables
#
# An active scenario objective equivalent to that generated by PySP is
# included here for informational purposes.


def total_cost_rule(model):
    return model.FirstStageCost + model.SecondStageCost
model.Total_Cost_Objective = Objective(rule=total_cost_rule, sense=minimize)


#
# Stochastic Data
#
Yield = {}
Yield['BelowAverageScenario'] = \
    {'WHEAT': 2.0, 'CORN': 2.4, 'SUGAR_BEETS': 16.0}
Yield['AverageScenario'] = \
    {'WHEAT': 2.5, 'CORN': 3.0, 'SUGAR_BEETS': 20.0}
Yield['AboveAverageScenario'] = \
    {'WHEAT': 3.0, 'CORN': 3.6, 'SUGAR_BEETS': 24.0}


def pysp_instance_creation_callback(scenario_name, node_names):

    instance = model.clone()
    instance.Yield.store_values(Yield[scenario_name])

    return instance


# model.pprint()

instance = model.create()
opt = SolverFactory("cbc")
results = opt.solve(instance, tee=True)
# model.pprint()
instance.load(results)
print [k for k in dir(instance)]
print instance.Total_Cost_Objective