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:
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:
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?).
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.
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:
# 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()
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