Challenge 1 - LP - Powerco (3 Plants, 4 Cities)

from Winston, Goldberg 2004, Pg 127. Nice one!

There are three power-plants and four cities. Each city's demand MUST be met. There is a limit to how much each power-plant CAN supply. There are costs associated with supply to a certain city from a certain power-plant. Find out what the supply breakup should be.

There is an elegant way and an inelegant way.

Inelegant way - use Set only - which means you need to spell out C1, C2, etc.. in a dict. And then you 

model.S1 = pyo.Var(model.i, within=pyo.NonNegativeReals)

S1 = model.S1

model.S2 = pyo.Var(model.i, within=pyo.NonNegativeReals)

S2 = model.S2

model.S3 = pyo.Var(model.i, within=pyo.NonNegativeReals)

S3 = model.S3

# To view final solution:

indices = list(model.i)

for var in [model.S1, model.S2, model.S3]:

    row = [f"{pyo.value(var[i]):.2f}" for i in indices]

    print(" ".join(row))

Elegant way: use RangeSet and then you specify the row and column variable names and maintain the weights (used to compute cost function) in a 2D matrix.

model.i = pyo.RangeSet(1,4) # cities

model.j = pyo.RangeSet(1,3) # plants

model.Price = pyo.Param( model.j, model.i, initialize = {(1,1):8, (1,2):6, (1,3):10, (1,4):9,

                                                         (2,1):9, (2,2):12, (2,3):13, (2,4):7,

                                                         (3,1):14, (3,2):9, (3,3):16, (3,4):5})

model.x = pyo.Var(model.j, model.i, within = pyo.NonNegativeReals)

x = model.x

I thought there might be some elegance involved in not having to spell out 7 constraints (capacity limit and demand requirement) but NS did the same thing as me.

To view final variable values:

for j in model.j:

    row = []

    for i in model.i:

        row.append(pyo.value(model.x[j, i]))

    print(row)

asdf

This is both elegant in that it's nicer than the first one, which uses only a single Set and multiple variable entry statements. But, don't you think Navid could do better than enter seven constraint functions manually? Come on..

chatGPT agrees:

Yes! You're absolutely right — Pyomo supports indexed constraints, which let you generalize patterns like this without needing to write separate functions for each one.

You can replace your 7 separately defined constraints with just two indexed constraints, one for supply and one for demand, using their indices as set elements.

model.S = pyo.Param(model.j, initialize={1: 35, 2: 50, 3: 40})  # supply at each plant

model.D = pyo.Param(model.i, initialize={1: 45, 2: 20, 3: 30, 4: 30})  # demand at each city

# Supply constraint: total shipped from each plant cannot exceed its supply

def SupplyRule(model, j):

    return sum(model.x[j, i] for i in model.i) <= model.S[j]

model.SupplyConstraint = pyo.Constraint(model.j, rule=SupplyRule)

# Demand constraint: total received by each city must meet demand

def DemandRule(model, i):

    return sum(model.x[j, i] for j in model.j) >= model.D[i]

model.DemandConstraint = pyo.Constraint(model.i, rule=DemandRule)