is it possible to get a subset of a 2-dimensional Set in pyomo

Question

I want to determin optimal charging for BEV along a grid line. Therefore I use pyomo to build a LP model describing the optimization. I have a Set of times (one day in quarter hour steps):

model.times = pe.Set(initialize=list(range(96)))

and a Set of nodes:

model.buses = pe.Set(initialize=list(range(6)))

I use these both sets to index the Variables and Constraints. For example, I have one Constraint tracking the SOCs of the BEVs:

def track_socs_rule(model, t, b):
    if t >= model.bevs[b].t_start and t <= model.bevs[b].t_target:
        return (model.SOC[t, b] + model.I[t, b] * model.voltages[b] * model.resolution/60 / 1000
                / model.bevs[b].e_bat*100 - model.SOC[t+1, b]) == 0
    else:
        return pe.Constraint.Skip

So i have to take care inside the rule defining the Constraint, that only valid timesteps are included (namely when the BEV at the according node is charging). It works, and I get proper results - but there has to be a more elegant way. I thought of defining a Set that only contains valid timesteps, something like:

model.valid_subset = pe.Set(within=model.times*model.charger_buses,
                            initialize=[[t for t in model.times if t >= model.bevs[b].t_start and
                                         t <= model.bevs[b].t_target] for b in model.charger_buses]

But this gives me an Error:

ValueError: Cannot add value([bla, bla, ...], [bla, bla...], ...) to Set valid_subset.
The value is not in the domain valid_subset_domain

Here is a minimal working example that should reproduce the error:

import pyomo.environ as pe

class BEV:
    def __init__(self, home_bus, e_bat, soc_start, soc_target, t_start, t_target, resolution):
        self.resolution = resolution
        self.home_bus = home_bus
        self.e_bat = e_bat
        self.soc_start = soc_start
        self.soc_target = soc_target
        self.t_start = int(t_start * 60/self.resolution)
        self.t_target = int(t_target * 60/self.resolution)


RESOLUTION = 15 # minutes

home_buses = [0, 2, 3, 4]
e_bats = [50, 50, 50, 50] # kWh
soc_starts = [20, 25, 20, 30] # %
soc_targets = [80, 75, 90, 100] # %
t_starts = [10, 12, 9, 15] # hours
t_targets = [19, 18, 15, 23] # hours

bevs = []
for pos in range(4):
    bev = BEV(home_bus=home_buses[pos], e_bat=e_bats[pos], soc_start=soc_starts[pos],
              soc_target=soc_targets[pos], t_start=t_starts[pos],
              t_target=t_targets[pos], resolution=RESOLUTION)
    bevs.append(bev)

solver = pe.SolverFactory('glpk')
model = pe.ConcreteModel()

model.bevs = {bev.home_bus: bev for bev in bevs}
model.resolution = RESOLUTION
model.times = pe.Set(initialize=list(range(int(24 * 60/RESOLUTION))))
model.buses = pe.Set(initialize=list(range(6)))
model.charger_buses = pe.Set(within=model.buses, initialize=[bus for bus in model.bevs])
#model.valid_subset = pe.Set(within=model.times*model.charger_buses,
                            #initialize=[[t for t in model.times if t >= model.bevs[b].t_start
                                         #and t <= model.bevs[b].t_target] for b in model.charger_buses])

model.voltages = pe.Param(model.buses, initialize={i: 400-i/2 for i in model.buses})

model.SOC = pe.Var(model.times*model.charger_buses, domain=pe.PositiveReals)
model.I = pe.Var(model.times*model.charger_buses, domain=pe.PositiveReals)

def track_socs_rule(model, t, b):
    # instead leave if (once indexed in model.valid_subset)
    if t >= model.bevs[b].t_start and t <= model.bevs[b].t_target:
        return (model.SOC[t, b] + model.I[t, b] * model.voltages[b] * model.resolution/60 /1000
                /model.bevs[b].e_bat*100 - model.SOC[t+1, b]) == 0

    else:
        return pe.Constraint.Skip

# instead index in model.valid_subset (once it works...)
model.track_socs = pe.Constraint(model.times*model.charger_buses, rule=track_socs_rule)
model.track_socs.pprint()

Once one uncomments the line with the model.valid_subset the error occurs. My hopings are, that once I use this "clever" index, I can leave away the if in track_socs_rule, because the components are allready properly indexed. This will also help me with the upper and lower bounds for I and SOC (they aren't even included in the code - i left them for the sake of simplicity).

So after a lot of explanation, my question:

is there a way to construct a subset of a 2-d Set in pyomo?

Many thanks for any help in advance!

Answer 1

It surely is a way to construct such a desired subset. In fact, your formulation is kinda fine, you need to refine some mistakes.

1. In the line model.valid_subset = pe.Set(within=model.times*model.charger_buses,... you're telling the model to constraint the subset to be within the cross product of model.times and model.charger_buses, but you're trying to initialize the subset just with t, therefore, you cannot add such a subset, since domain is (t_i,b_i) and you're just adding t_i.

2. For a Set component, initialize arg should be a list of actual values (I think that any np.array-like iterable should work, but not fully sure), but you're passing a list of lists, then, and error will raise. You need to flatten your list of lists into a single list.

To flatten your list of list, I will be using itertools as in this answer, but you can do it the way you want it.

import itertools
#All your model
...
model.valid_subset = pe.Set(within=model.times*model.charger_buses,
                            initialize=list(itertools.chain(*[[(t,b) for t in model.times if t >= model.bevs[b].t_start
                                         and t <= model.bevs[b].t_target] for b in model.charger_buses])))

You may want to construct any sub-set outside the actual initialize arg. It may be the same list comprenhension, but not in the actual arg of Set (or any) component. This would help to make it more readable.