I have an energy allocation problem which has been implemented in Pyomo. The objective is to allocate the available energy (which is limited and time-varying) across a set of requests in a way which maximises the utilisation of the available energy, the problem is further described here: https://stackoverflow.com/questions/77982332/generating-constraints-with-conditional-summing-over-several-pyomo-subsets/77982847#77982847.
Each request has power (kW) and energy (kWh) requirements and constraints on the timeframe in which it can be satisfied.
I have implemented this successfully (with much thanks to @Airsquid) for 99% of requests, but I have several cases where some requests are seemingly never being satisfied and I cannot figure out why this is.
I have reproduced a minimal example below. For context, the set of time is over 24 hours with a 30 minute resolution (48 periods indexed from 910176 to 910223 inclusive).
The single request in the example below requires 3.25kW of power delivered over 10 periods (5 hours) to receive 16.25kWh, and it needs to receive this energy in a single block within the interval from 910177 to 910211 (which is a 35 period interval or 17.5 hours).
When this code is executed, the request is not satisfied i.e. m.satisfied == 0 after the execution. However, if I change the latest end time from 910211 to 910205 (i.e. reduce the valid interval to deliver the energy within from 35 periods to 29 periods), the request does become satisfied i.e. m.satisfied == 1.
Making this change does not make sense to me as by reducing the latest end time, we are theoretically making the problem more difficult to solve.
As mentioned, when I execute my market, 99% of requests are satisfied with this code as expected, but there are some which are not satisfied even when I execute the market with only the problematic requests included. I have not been able to identify anything specific to do with the requests which may be causing issues, for example there is a mix of large/small power requests, large/small energy requests, large/small windows over which the energy is requested.
At this point I'm very much at a loss why this is happening and I'm wondering if it is something in my implementation, or if it's possible that this is somehow a problem relating to Pyomo. Any suggestions that can be made on how to move past this problem would be sincerely appreciated, and general suggestions on best practise for how to improve the implementation also welcome!
import pyomo.environ as pyo
import sys
import time
from collections import defaultdict
import pandas as pd
import numpy as np
sampling_period_hrs=0.5 #our market uses a time resolution of 30 minutes / 0.5 hours, we use this later to relate request power delivered to request energy required
request_time_constraints_dict = {}
request_max_power_dict = {}
request_max_energy_dict={}
request_duration_num_periods_dict = {}
#this is the available supply (kW) to satisfy requests
supply_index = list(range(910176, 910223+1)) #this is the timeframe over which power is available to be allocated
supply_vals = np.ones(len(supply_index))*2500000 #this is the amount of power (kW) available, this has been made arbitrarily large for the purposes of this illustration
supply_unix_ts = pd.Series(supply_vals, index=supply_index)
request_earliest_start_time = 910177 #earliest start time request can be satisfied
request_latest_end_time = 910211 #latest end time request can be satisfied, TEST: change this to 910205 to see request m.satisfied goes to 1
request_time_constraints_dict[0] = (int(request_earliest_start_time), int(request_latest_end_time))
request_max_power_dict[0] = 3.25 #kw
request_max_energy_dict[0] = 16.25 #kwh
request_duration_num_periods_dict[0] = 10.0 #periods (10 x 0.5h periods = 5 hours)
eligible_requests = defaultdict(list)
for r, (start, end) in request_time_constraints_dict.items():
for t in range(start, end + 1):
eligible_requests[t].append(r)
### MODEL BUILD
m = pyo.ConcreteModel('dispatch')
### SETS
# m.T = pyo.Set(initialize=tuple(supply_unix_ts.index))
m.T = pyo.Set(initialize=tuple(range(910176,910223+1)))
m.R = pyo.Set(initialize=tuple(request_max_power_dict.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.power_limit = pyo.Param(m.T, initialize=supply_unix_ts.to_dict()) #total amount of power (kW) available for a specific period
m.request_power_size = pyo.Param(m.R, initialize=request_max_power_dict) #amount of power (kW) required for a request
m.request_energy_size = pyo.Param(m.R, initialize=request_max_energy_dict) #amount of energy (kWh) required for a request
m.duration_periods = pyo.Param(m.R, initialize=request_duration_num_periods_dict) #number of periods (note this is calculated by relating the energy to the power to ensure there are no conflicts)
### VARS
m.booked_supply = pyo.Var(m.T, domain = pyo.NonNegativeReals, doc='total amount of power used in a time period')
m.dispatch = pyo.Var(m.windows_flat, domain=pyo.NonNegativeReals,
doc='dispatch power in timeslot t to request r')
m.start = pyo.Var(m.windows_flat, domain=pyo.Binary, doc='binary for period in which request starts to be satisfied')
m.running = pyo.Var(m.windows_flat, domain=pyo.Binary, doc='binary for if request being satisfie in interval')
m.satisfied = pyo.Var(m.R, domain=pyo.Binary, doc='binary for if request is fully satisfied')
### OBJ
def obj_rule(m):
return sum(m.booked_supply[t] for t in m.T)
m.obj = pyo.Objective(rule = obj_rule, sense = pyo.maximize)
### CONSTRAINTS
#constraint 1
@m.Constraint(m.T)
def supply_limit(m, t):
return m.booked_supply[t] <= m.power_limit[t]
#constraint 2
@m.Constraint(m.T)
def booked_supply_limit(m, t):
# total power delivered in each interval
return sum(m.running[t,r]*m.request_power_size[r] for r in m.windows[t]) == m.booked_supply[t]
#constraint 3
@m.Constraint(m.R)
def request_satisfied(m, r):
# dictates whether a request is satisfied and relates power to energy
timeslots = {t for t in m.T if r in m.windows[t]}
return sum(m.running[t, r]*m.request_power_size[r] for t in timeslots)*sampling_period_hrs == m.satisfied[r] * m.request_energy_size[r]
# constraint 4
@m.Constraint(m.R)
def run_once(m, r):
#request can only run once
return(sum(m.start[t,r] for t in m.T if r in m.windows[t]) <= 1)
#constraint 5
@m.Constraint(m.windows_flat)
def link_running(m, t, r):
#request can only be satisfied if it was running in prior period or started in this one (we need request to be satisfied in a single block rather than turning off/on)
timeslots = {t for t in m.T if r in m.windows[t]}
if t <= list(timeslots)[0]:
return m.running[t,r] == m.start[t,r]
return m.running[t,r] <= m.running[t-1,r] + m.start[t,r]
# constraint 6
@m.Constraint(m.windows_flat)
def keep_track_power(m, t, r):
#useful for debugging - keep track of actual power values for each request in each time interval
return m.running[t, r]*m.request_power_size[r] == m.dispatch[t, r]
# solve it
solver = pyo.SolverFactory('cbc')
res = solver.solve(m)
temp = sys.stdout
f = open(str(int(time.time()))+'.txt', 'w')
sys.stdout = f
m.pprint()
print(res)
f.close()
sys.stdout = temp
You have an error in your link_running constraint that is causing all the trouble. If you look at some of the output from the failed runs, you'll note in the printout of that constraint, the "first period" equation for running = start seems to show up all over the place....
That's because you are trying to get the first eligible start time in the window with this statement:
if t <= list(timeslots)[0]:
timeslots is a set so there is absolutely no guarantee that the [0] element of the list conversion is the lowest value. Recall that a set is unordered (with some caveats....but this will generally save your bacon.) You want to do this:
if t == min(timeslots):
to properly identify the first eligible period.
On a second, unrelated note, you could probably re-factor this and make it simpler after you get the current form working. If desired. If it runs fine and timely, then don't bother.
However, it appears to be the case in your model that the amount of power per request is fixed in all the periods used to satisfy it. Is that true? Or is your intent to make a more flexible model where the power allocated per request can vary, and thus the number of periods used to fulfill the request can vary? Those are 2 way different models.
In your current build, it appears to be the former. In that case, you "know everything about the request" just from the start variable. So, you could do a big refactor and probably get rid of some of the helper variables. If it starts in period 5 and the request is for 10 periods with a certain fixed amount of power, You don't need a variable to know if it is running in period 8--it is! Similarly for the allocation within the period. You'd still need a variable to keep track of how much power is allocated to each request.
Just something to think about--