Search code examples
pythonpyomoamplglpk

Converting from AMPL to Pyomo


I'm trying to convert an AMPL model to Pyomo (something I have no experience with using). I'm finding the syntax hard to adapt to, especially the constraint and objective sections. I've already linked my computer together with python, anaconda, Pyomo, and GLPK, and just need to get the actual code down. I'm a beginner coder so forgive me if my code is poorly written. Still trying to get the hang of this!

Here is the data from the AMPL code:

set PROD := 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30;

set PROD1:= 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30;

ProdCost    414 3   46  519 876 146 827 996 922 308 568 176 58  13  20  974 121 751 130 844 280 123 275 843 717 694 72  413 65  631

HoldingCost 606 308 757 851 148 867 336 44  364 960 69  428 778 485 285 938 980 932 199 175 625 513 536 965 366 950 632 88  698 744

Demand  105 70  135 67  102 25  147 69  23  84  32  41  81  133 180 22  174 80  24  156 28  125 23  137 180 151 39  138 196 69

And here is the model:

set PROD;  # set of production amounts
set PROD1; # set of holding amounts

param ProdCost {PROD} >= 0;     # parameter set of production costs

param Demand {PROD} >= 0;     # parameter set of demand at each time

param HoldingCost {PROD} >= 0;     # parameter set of holding costs

var Inventory {PROD1} >= 0;     # variable that sets inventory amount at each time

var Make {p in PROD} >= 0;  # variable of amount produced at each time

minimize Total_Cost: sum {p in PROD} ((ProdCost[p] * Make[p]) + (Inventory[p] * HoldingCost[p]));

               # Objective: minimize total cost from all production and holding cost

subject to InventoryConstraint {p in PROD}: Inventory[p] = Inventory[p-1] + Make[p] - Demand[p];

                # excess production transfers to inventory

subject to MeetDemandConstraint {p in PROD}: Make[p] >= Demand[p] - Inventory[p-1];

               # Constraint: holding and production must exceed demand

subject to InitialInventoryConstraint: Inventory[0] = 0;

                # Constraint: Inventory must start at 0

Here's what I have so far. Not sure if it's right or not:

from pyomo.environ import *

demand=[105,70,135,67,102,25,147,69,23,84,32,41,81,133,180,22,174,80,24,156,28,125,23,137,180,151,39,138,196,69]

holdingcost=[606,308,757,851,148,867,336,44,364,960,69,428,778,485,285,938,980,932,199,175,625,513,536,965,366,950,632,88,698,744]

productioncost=[414,3,46,519,876,146,827,996,922,308,568,176,58,13,20,974,121,751,130,844,280,123,275,843,717,694,72,413,65,631]

model=ConcreteModel()

model.I=RangeSet(1,30)
model.J=RangeSet(0,30)
model.x=Var(model.I, within=NonNegativeIntegers)
model.y=Var(model.J, within=NonNegativeIntegers)

model.obj = Objective(expr = sum(model.x[i]*productioncost[i]+model.y[i]*holdingcost[i] for i in model.I))

def InventoryConstraint(model, i):
    return model.y[i-1] + model.x[i] - demand[i] <= model.y[i]
InvCont = Constraint(model, rule=InventoryConstraint)

def MeetDemandConstraint(model, i):
    return model.x[i] >= demand[i] - model.y[i-1]
DemCont = Constraint(model, rule=MeetDemandConstraint)

def Initial(model):
    return model.y[0] == 0
model.Init = Constraint(rule=Initial)

opt = SolverFactory('glpk')
results = opt.solve(model,load_solutions=True)
model.solutions.store_to(results)
results.write()

Thanks!


Solution

  • The only issues I see are in some of your constraint declarations. You need to attach the constraints to the model and the first argument passed in should be the indexing set (which I'm assuming should be model.I).

    def InventoryConstraint(model, i):
        return model.y[i-1] + model.x[i] - demand[i] <= model.y[i]
    model.InvCont = Constraint(model.I, rule=InventoryConstraint)
    
    def MeetDemandConstraint(model, i):
        return model.x[i] >= demand[i] - model.y[i-1]
    model.DemCont = Constraint(model.I, rule=MeetDemandConstraint)
    

    The syntax that you're using to solve the model is a little out-dated but should work. Another option would be:

    opt = SolverFactory('glpk')
    opt.solve(model,tee=True) # The 'tee' option prints the solver output to the screen
    model.display() # This will print a summary of the model solution
    

    Another command that is useful for debugging is model.pprint(). This will display the entire model including the expressions for Constraints and Objectives.