Generating constraints with conditional summing over several pyomo subsets

Question

New to pyomo and to linear optimisation generally. I have an energy allocation problem where I have a fixed amount of energy over a timeframe, I have defined the relevant timeframe as a set:

model.t = pyomo.Set(initialize = (supply_unix_ts.index)) 

And the total amount of energy available for each time period is defined as a parameter:

model.available_supply = pyomo.Param(model.t, initialize = supply_unix_ts.to_dict(), domain=pyomo.NonNegativeReals) 

I then have some requests to consume the available energy defined as another set:

model.r = pyomo.Set(initialize = request_ids) 

Each request has some specific parameters associated with it, specifically:

• The total amount of energy needed (kWh)
• The maximum power the request can be provided with (kW)
• The earliest time that the request can start to receive energy and the latest time the request can receive energy by (this interval is different for each request and is a subset of t)

 model.request_earliest_start_time = pyomo.Param(model.r, initialize = request_earliest_start_times_dict, domain=pyomo.NonNegativeReals) #parameter for earliest start times

        model.request_latest_end_time = pyomo.Param(model.r, initialize = request_latest_end_time_dict, domain=pyomo.NonNegativeReals) #parameter for latest end times

        model.request_max_power_dict = pyomo.Param(model.r, initialize = request_max_power_dict, domain=pyomo.NonNegativeReals) #parameter for max power

        model.request_energy_needed_dict = pyomo.Param(model.r, initialize = request_energy_needed_dict, domain=pyomo.NonNegativeReals)

I also have another parameter for the sampling period which I need to convert between power and energy in some places to follow.

model.sampling_period_dict = pyomo.Param(model.t, initialize = sampling_period_dict, domain=pyomo.NonNegativeReals)

My objective is to maximise the utilisation of the available energy, not all requests need to be satisfied.

    def obj_rule(model):
                return sum(model.booked_supply[t] for t in model.t())
            model.obj = pyomo.Objective(rule = obj_rule, sense = pyomo.maximize)

The decision variables I have are as follows, two of which are binary:

        model.booked_supply = pyomo.Var(model.t, domain = pyomo.NonNegativeReals) #total supply that is used for a specific time
        model.request_status = pyomo.Var(model.t, model.r, domain = pyomo.Binary) #if request is allocated energy during time period t
        model.request_satisfied = pyomo.Var(model.r, domain = pyomo.Binary) #if request is fully satisfied

I have implemented a constraint to ensure that the booked supply (model.booked_supply) cannot be greater than the available supply for each time period (model.available_supply), and a constraint to set the booked supply equal to the sum of the requests satisfied for each time period.

The problem that I am having is implementing a constraint to say that the total amount of energy allocated to a request should be equal to the total amount of energy needed by that request (model.request_energy_needed_dict), within the interval specified by the request (i.e. between model.request_earliest_start_time and model.request_latest_end_time) i.e. this constraint: Math form of constraint

Where [E_r,L_r] is the time interval relevant for the request.

I have tried looking into lots of different approaches but without any success, since I am new I am unclear about both what the best approach is, and the required syntax to implement this (and if I have defined everything correctly in pyomo).

My two main approaches have been:

1. Instead of defining [E_r and L_r] (i.e. model.request_earliest_start_time and model.request_latest_end_time ) as parameters, instead try to define these as subsets of the model.t Set and then implement an individual Constraint for each request using the relevant subset. The problem with this is that I have potentially 1000s of requests, I have seen it is possible to dynamically create constraints using ConstraintList but have not seen a way to dynamically create large numbers of individual sets like a SetList (doesn't seem to exist?).

2. I hoped it would be possible to implement this in a single constraint, but so far have struggled to write down how this should work. Some very naive attempts are below which have generated various errors...

def requestMaxEnergyConstraint(model, t, r):
return sum(model.request_confirmed_power[r,t]*model.request_satisfied[r,t]*model.sampling_period_dict[t] for t in model.t if (t > model.request_earliest_start_time[r] and t < model.request_latest_end_time[r]))== model.request_max_energy_dict[r]*model.request_satisfied[r] 

def requestMaxEnergyConstraint(model, t, r):
    maxEnergy=0
    for r in model.r:
        for t in model.t:
            if((t > model.request_earliest_start_time[r,t]) and t < (model.request_earliest_start_time[r,t])):
                maxEnergy += model.request_confirmed_power[r,t]
    return maxEnergy == model.request_max_energy_dict[r]*model.request_satisfied[r] 

Any advice (ideally a worked example) would be really appreciated on how this can be achieved.

Answer 1

Here are a couple of concepts that I think will help you out....

First, I'm quite sure you are going to need to index your power delivery variable by both request and timeslice because you will need to know in your constraints how much power is going to individual requests in each timeslice.

You should also use an indexed set, which is kinda like a python dictionary, where you have a set of sets that is indexed by a set, in your case, you could go for the set of periods by request or vice-versa. Or alternatively, you could just hold the min/max period in 2 different parameters, and create a list of eligible time periods within a constraint by filtering or such.

If you go the way shown below, then you can make a flat set of just the request - time slot pairs that are "legal" for consideration. This has a HUGE benefit if your model grows because you will be omitting all of the ineligible pairs within the power supply variable, which will be "sparse" instead of the full R x T set.

Anyhow, there are probably a couple variants of this, but this gets it done pretty efficiently:

Code:

# energy assignment
from collections import defaultdict

import pyomo.environ as pyo

### DATA

periods = tuple(range(4))
period_limit = 20

# (earliest, latest)
request_timespan = {
    'r1': (0, 3),
    'r2': (2, 3)
}
request_power = {
    'r1': 60,
    'r2': 10
}
# prevent mind-numbing errors:
assert request_timespan.keys() == request_power.keys()

# we know that we are going to need to sum power across each period at some point, so one
# approach is to determine which request are eligible to be satisfied in each period.  So,
# you could probably either make a dictionary of requests: {timeslots} or perhaps better:
# timeslot: {eligible requests}... either way can work.

eligible_requests = defaultdict(list)
for r, (start, end) in request_timespan.items():
    for t in range(start, end + 1):
        eligible_requests[t].append(r)

# check it...
print(eligible_requests)

### MODEL BUILD

m = pyo.ConcreteModel('dispatch')

### SETS
m.T = pyo.Set(initialize=periods)
m.R = pyo.Set(initialize=tuple(request_power.keys()))
# Note:  This will be an "indexed" set, where we have sets indexed by some index, in this case, m.T
m.windows = pyo.Set(m.T, initialize=eligible_requests, within=m.R,
                    doc='requests eligible in each timeslot')
# We can make a flat-set here which will be a good basis for the dispatch variable
m.windows_flat = pyo.Set(initialize={(t, r) for t in eligible_requests for r in
                                     eligible_requests[t]},
                         within=m.T * m.R)

### PARAMS
m.request_size = pyo.Param(m.R, initialize=request_power)
m.power_limit = pyo.Param(m.T, initialize=period_limit)

### VARS
m.satisfied = pyo.Var(m.R, domain=pyo.Binary)
m.dispatch = pyo.Var(m.windows_flat, domain=pyo.NonNegativeReals,
                     doc='dispatch power in timeslot t to request r')

### OBJ
m.obj = pyo.Objective(expr=sum(m.satisfied[r] for r in m.R), sense=pyo.maximize)


### CONSTRAINTS
@m.Constraint(m.R)
def request_satisfied(m, r):
    # this is kind of a reverse-lookup in the dict, but it works...  note, we could have
    # made a 2nd dictionary that works in reverse, but it wouldn't make things much faster
    # and this demo's an alternate approach.  Won't make much difference for overall time...
    timeslots = {t for t in m.T if r in m.windows[t]}
    return sum(m.dispatch[t, r] for t in timeslots) >= m.satisfied[r] * m.request_size[r]


@m.Constraint(m.T)
def supply_limit(m, t):
    return sum(m.dispatch[t, r] for r in m.windows[t]) <= m.power_limit[t]


# check it
m.pprint()

# solve it
solver = pyo.SolverFactory('cbc')
res = solver.solve(m)
print(res)

m.dispatch.display()

Output:

defaultdict(<class 'list'>, {0: ['r1'], 1: ['r1'], 2: ['r1', 'r2'], 3: ['r1', 'r2']})
5 Set Declarations
    R : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    2 : {'r1', 'r2'}
    T : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    4 : {0, 1, 2, 3}
    windows : requests eligible in each timeslot
        Size=4, Index=T, Ordered=Insertion
        Key : Dimen : Domain : Size : Members
          0 :     1 :      R :    1 : {'r1',}
          1 :     1 :      R :    1 : {'r1',}
          2 :     1 :      R :    2 : {'r1', 'r2'}
          3 :     1 :      R :    2 : {'r1', 'r2'}
    windows_flat : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain              : Size : Members
        None :     2 : windows_flat_domain :    6 : {(0, 'r1'), (2, 'r2'), (2, 'r1'), (3, 'r2'), (3, 'r1'), (1, 'r1')}
    windows_flat_domain : Size=1, Index=None, Ordered=True
        Key  : Dimen : Domain : Size : Members
        None :     2 :    T*R :    8 : {(0, 'r1'), (0, 'r2'), (1, 'r1'), (1, 'r2'), (2, 'r1'), (2, 'r2'), (3, 'r1'), (3, 'r2')}

2 Param Declarations
    power_limit : Size=4, Index=T, Domain=Any, Default=None, Mutable=False
        Key : Value
          0 :    20
          1 :    20
          2 :    20
          3 :    20
    request_size : Size=2, Index=R, Domain=Any, Default=None, Mutable=False
        Key : Value
         r1 :    60
         r2 :    10

2 Var Declarations
    dispatch : dispatch power in timeslot t to request r
        Size=6, Index=windows_flat
        Key       : Lower : Value : Upper : Fixed : Stale : Domain
        (0, 'r1') :     0 :  None :  None : False :  True : NonNegativeReals
        (1, 'r1') :     0 :  None :  None : False :  True : NonNegativeReals
        (2, 'r1') :     0 :  None :  None : False :  True : NonNegativeReals
        (2, 'r2') :     0 :  None :  None : False :  True : NonNegativeReals
        (3, 'r1') :     0 :  None :  None : False :  True : NonNegativeReals
        (3, 'r2') :     0 :  None :  None : False :  True : NonNegativeReals
    satisfied : Size=2, Index=R
        Key : Lower : Value : Upper : Fixed : Stale : Domain
         r1 :     0 :  None :     1 : False :  True : Binary
         r2 :     0 :  None :     1 : False :  True : Binary

1 Objective Declarations
    obj : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : maximize : satisfied[r1] + satisfied[r2]

2 Constraint Declarations
    request_satisfied : Size=2, Index=R, Active=True
        Key : Lower : Body                                                                                   : Upper : Active
         r1 :  -Inf : 60*satisfied[r1] - (dispatch[0,r1] + dispatch[1,r1] + dispatch[2,r1] + dispatch[3,r1]) :   0.0 :   True
         r2 :  -Inf :                                   10*satisfied[r2] - (dispatch[2,r2] + dispatch[3,r2]) :   0.0 :   True
    supply_limit : Size=4, Index=T, Active=True
        Key : Lower : Body                            : Upper : Active
          0 :  -Inf :                  dispatch[0,r1] :  20.0 :   True
          1 :  -Inf :                  dispatch[1,r1] :  20.0 :   True
          2 :  -Inf : dispatch[2,r1] + dispatch[2,r2] :  20.0 :   True
          3 :  -Inf : dispatch[3,r1] + dispatch[3,r2] :  20.0 :   True

12 Declarations: T R windows windows_flat_domain windows_flat request_size power_limit satisfied dispatch obj request_satisfied supply_limit

Problem: 
- Name: unknown
  Lower bound: 2.0
  Upper bound: 2.0
  Number of objectives: 1
  Number of constraints: 4
  Number of variables: 6
  Number of binary variables: 2
  Number of integer variables: 2
  Number of nonzeros: 2
  Sense: maximize
Solver: 
- Status: ok
  User time: -1.0
  System time: 0.0
  Wallclock time: 0.0
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
    Black box: 
      Number of iterations: 0
  Error rc: 0
  Time: 0.00905919075012207
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

dispatch : dispatch power in timeslot t to request r
    Size=6, Index=windows_flat
    Key       : Lower : Value : Upper : Fixed : Stale : Domain
    (0, 'r1') :     0 :  20.0 :  None : False : False : NonNegativeReals
    (1, 'r1') :     0 :  20.0 :  None : False : False : NonNegativeReals
    (2, 'r1') :     0 :  10.0 :  None : False : False : NonNegativeReals
    (2, 'r2') :     0 :  10.0 :  None : False : False : NonNegativeReals
    (3, 'r1') :     0 :  10.0 :  None : False : False : NonNegativeReals
    (3, 'r2') :     0 :   0.0 :  None : False : False : NonNegativeReals