How to implement conditional summing within Pyomo constraint

Question

I am writing a model in Pyomo that optimizes for biomass production across ~100 different biomass types (model.Biomass = corn, wheat, wood scraps, etc.) on a county basis (model.SourceCounty). One of the constraints I am trying to write requires that the biomass output from my model equals production values that I've already obtained from another model. This other model does not have the granularity that my Pyomo model will have, however. It has biomass production values only on a regional (not county) basis (model.Zone) across more general biomass groupings (model.SimpBiomass = herbaceous biomass, woody biomass).

What I am trying to do in my constraint is sum up the biomass production decision variable (model.x) over the regions and more general biomass groupings from the other model before requiring that this sum equals the output from the other model so that my model produces a consistent result. However, what I'm learning is that the current way I've written the code (below) doesn't work because Pyomo calls constraints only once, when the value of the decision variables is yet to be solved for. Thus, my for loops with if statements just return a value of 0.

from pyomo.environ import *

# initialize model -- can rewrite as Concrete Model if that's better for what I'm trying to do
model = AbstractModel()

# initialize indices, including some manually
model.SourceCounty = Set() # county
model.Biomass = Set() # biomass type in my model
model.Year = Set(initialize=[2022, 2024, 2026, 2028, 2030, 2032, 2035, 2040])  # year
model.SimpBiomass = Set(initialize=['herbaceous biomass', 'waste biomass', 'woody biomass'])  # type of feedstock resource - simplified (from separate model)
model.Zone = Set(initialize=['midwest','southeast'])  # zones from separate model

# Create data import structure
data = DataPortal()

# Load indices that require data frame
data.load(filename='fips.csv', set=model.SourceCounty)
data.load(filename='Resources.csv', set=model.Biomass)

# initialize parameters
model.EERF = Param(model.SimpBiomass, model.Zone, model.Year)  # biomass production from the other model that I'm trying to match in my model
model.QIJ = Param(model.SourceCounty)  # mapping of county FIPS code from my model to zones from other model
model.AC = Param(model.Biomass)  # mapping of specific resource type from my model into less specific from other model (values are those in SimpBiomass)

# load in parameters
data.load(filename="county_to_zone.csv", param=model.QIJ)
data.load(filename="BT16_to_zone_Resource.csv", param=model.AC)

# create decision variables (known as Var in Pyomo)
model.x = Var(model.Biomass, model.SourceCounty, model.Year, domain=PositiveReals)  # feedstock production indexed by feedstock, source county, year

# I leave out the objective function for brevity

# Constraint in question
def feedstock_prod_rule(model, c, q, t):
    expr2 = 0 # initialize summing variable
    
    # for each biomass type (a) in my model, check if it belongs to a biomass category (c) from the other model
    
    for a in model.Biomass:
        if model.AC[a] == c:

            # for each county (i) in my model, check if it belongs to a zone (q) from the other model
            for i in model.SourceCounty:
                if model.QIJ[i] == q:

                    # if it belongs to q and c from other model, add to expr2
                    expr2 += model.x[a, i, t]
    
    # Sum of all biomass production from my model within zone q and biomass type c (expr2 at end of looping) should equal the output of the other model (EERF).

    return expr2 == model.EERF[c, q, t]


# Add as constraint
model.feedstock_prod = Constraint(model.SimpBiomass, model.Zone, model.Year, rule=feedstock_prod_rule)

I need help figuring out a different way to write this constraint such that it doesn't rely on building up an expression that depends on the value of my decision variable model.x that has yet to be solved for. Is there a way to have one line of code in the return line that accomplishes the same thing?

user_10789207 · Answer

I'm not totally sure you diagnosed the problem correctly, and without some snippets of your data files, it is too difficult to recreate. (It is also unnecessary, because even if it was re-creatable, there is a better way. :) )

In general, you are correct that you cannot embed conditional statements into constraints that depend on the value of a variable, which is unknown when the constraint is encoded into the model. However, you are using conditionals based on parameters which are fixed, so that should be OK. However, you are comparing parameter values to items in a set, which is.... a bad plan even if it worked... never tried it. The structural problem with that is that you are working with 1:1 pairings instead of labeling groups of stuff, which is leading to the conditional statements. You have a structure like:

beans : food
lettuce : food
paper : trash

Where you'd really be happier to work with groups like below and avoid the "if" statements:

food: { beans, lettuce }
trash: { paper }

So this can be done, and you can load it into an abstract model. I don't think you can do it from a .csv file, however, as I don't think there is a way to express indexed sets in .csv. You can easily do it in .yaml or .json. See the pyomo dox for more examples. You can even commingle your data sources, so you can retain your other csv's, as long as they are consistent, so you need to be aware of that. Here is a working example, that I think will clean up your model a bunch. Specifically note the indexing and groupings in the constraint at the end):

import pyomo.environ as pe

m = pe.AbstractModel()

# SETS
m.Group_names = pe.Set()
m.Items = pe.Set()
m.Groupings = pe.Set(m.Group_names, within=m.Items)

# PARAMS
m.Cost = pe.Param(m.Group_names)

# VARS
m.X = pe.Var(m.Items)

# Constraint "for each Group"
def C1(m, g_name):
    return sum(m.X[i] for i in m.Groupings[g_name]) <= 10

m.C1 = pe.Constraint(m.Group_names, rule=C1)

# load data from sources
data = pe.DataPortal()
data.load(filename='cost.csv', param=m.Cost)
data.load(filename='data.yaml')
data.load(filename='items.csv', set=m.Items)

instance = m.create_instance(data)
instance.pprint()

The `yaml` file (others not shown, they are straightforward):

Group_names: ['Waste', 'Food', 'Critters']

Groupings:
  'Waste': ['Trash', 'Coal']
  'Food': ['Salad', 'Compost', 'Fish']
  'Critters': ['Snails', 'Worms']

The resultant model:

3 Set Declarations
    Group_names : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    3 : {'Waste', 'Food', 'Critters'}
    Groupings : Size=3, Index=Group_names, Ordered=Insertion
        Key      : Dimen : Domain : Size : Members
        Critters :     1 :  Items :    2 : {'Snails', 'Worms'}
            Food :     1 :  Items :    3 : {'Salad', 'Compost', 'Fish'}
           Waste :     1 :  Items :    2 : {'Trash', 'Coal'}
    Items : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    7 : {'Trash', 'Snails', 'Worms', 'Coal', 'Salad', 'Compost', 'Fish'}

1 Param Declarations
    Cost : Size=2, Index=Group_names, Domain=Any, Default=None, Mutable=False
        Key   : Value
         Food :   8.9
        Waste :   4.2

1 Var Declarations
    X : Size=7, Index=Items
        Key     : Lower : Value : Upper : Fixed : Stale : Domain
           Coal :  None :  None :  None : False :  True :  Reals
        Compost :  None :  None :  None : False :  True :  Reals
           Fish :  None :  None :  None : False :  True :  Reals
          Salad :  None :  None :  None : False :  True :  Reals
         Snails :  None :  None :  None : False :  True :  Reals
          Trash :  None :  None :  None : False :  True :  Reals
          Worms :  None :  None :  None : False :  True :  Reals

1 Constraint Declarations
    C1 : Size=3, Index=Group_names, Active=True
        Key      : Lower : Body                            : Upper : Active
        Critters :  -Inf :            X[Snails] + X[Worms] :  10.0 :   True
            Food :  -Inf : X[Salad] + X[Compost] + X[Fish] :  10.0 :   True
           Waste :  -Inf :              X[Trash] + X[Coal] :  10.0 :   True

How to implement conditional summing within Pyomo constraint

1 Answers

The yaml file (others not shown, they are straightforward):

The resultant model:

The `yaml` file (others not shown, they are straightforward):