Defining the value of one variable in a constraint in relation to another variable without making the problem nonlinear in Pyomo
Background references: [1] Generating constraints with conditional summing over several pyomo subsets
[2] Debugging why some requests are unexpectedly not satisfied in pyomo
I am getting a correct solution for the problem detailed above with the solutions provided by @Airsquid, however I am facing long runtimes for the case where there is a small amount of energy supply available (relative to the number of requests).
I have added some early stopping settings, which is fine for the majority of cases as I don't necessarily need the totally optimal solution just one which is good enough. solver.options = { 'sec': 600, 'threads':6, 'ratio': 5} #early stopping settings
I wanted to try to refactor the code per the suggestion given by @Airsquid in [2]: "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."
Since I am still new to linear programming and pyomo, I am having some problems with formulating the syntax for this.
I have only been able to think of two approaches to do this:
Approach 1: Define the value of the m.dispatch[t,r] variable between the time when m.start[t,r] == 1 and m.start[t+request_length_periods,r], a crude attempt at the syntax for this is below! I know this is not correct and why it is not correct but I don't know what the solution for it is.
@m.Constraint(m.windows_flat)
def request_length(m, t, r):
request_length_periods = int(m.duration_periods[r])
request_periods = set(list(range(t,t+request_length_periods+1)))
return (m.dispatch[t,r] == m.request_power_size[r] for t in request_periods if m.start[t,r] == 1)
This unsuprisingly gives me an error:
ERROR: Rule failed when generating expression for Constraint request_length
with index (910192, 88): ValueError: Constraint 'request_length[910192,88]'
does not have a proper value. Found '<generator object
generalOptimisation.<locals>.request_length.<locals>.<genexpr> at
0x168022b20>' Expecting a tuple or relational expression. Examples:
sum(model.costs) == model.income (0, model.price[item], 50)
ERROR: Constructing component 'request_length' from data=None failed:
ValueError: Constraint 'request_length[910192,88]' does not have a proper
value. Found '<generator object
I know how I have defined request_periods is not correct, but I don't know how to define this set without knowing when m.start[t,r] == 1?
Approach 2: m.start[t,r] *(m.dispatch[t,r]...m.dispatch[t+request_length_periods,r]) == 1 (I don't know what the syntax is for this in Pyomo and I know this would make the solution nonlinear)
If anyone has any suggestions on how to formulate this correctly per @Airsquid's suggestion, I would be very grateful, currently I am planning to relax the early stopping further to work around this which is not entirely desirable.
Any input on the early stopping settings is also welcome - perhaps the system time out at 600 seconds is not realistic!
Here's a refactor of your original code that uses 1 binary variable for the dispatch. There's a couple of things in here that could probably be refactored a bit, but it shows one approach to this. It also makes the assumption that the power delivered to each request is constant across periods.
Captured a rather unique result after a few pings (shown)
Code
# energy assignment
import sys
from collections import defaultdict
from random import randint
import pyomo.environ as pyo
from matplotlib import pyplot as plt
### DATA
num_periods = 12
periods = tuple(range(num_periods))
# (earliest, latest)
request_timespan = {
"r1": (0, 7),
"r2": (2, 3),
"r3": (1, 5),
}
request_power_per_period = {"r1": 5, "r2": 4, "r3": 8}
request_duration = {"r1": 4, "r2": 8, "r3": 5}
power_avail = dict(zip(periods, (randint(8, 20) for _ in range(num_periods))))
# prevent mind-numbing errors:
assert request_timespan.keys() == request_power_per_period.keys()
assert request_timespan.keys() == request_duration.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.
eiligible_starts = defaultdict(list)
for r, (start, end) in request_timespan.items():
for t in [
t for t in range(start, end + 1) if t + request_duration[r] <= num_periods
]:
eiligible_starts[t].append(r)
# check it...
print(eiligible_starts)
### MODEL BUILD
m = pyo.ConcreteModel("dispatch")
### SETS
m.T = pyo.Set(initialize=periods)
m.R = pyo.Set(initialize=tuple(request_timespan.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=eiligible_starts,
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 eiligible_starts for r in eiligible_starts[t]},
within=m.T * m.R,
)
### PARAMS
m.request_size = pyo.Param(m.R, initialize=request_power_per_period)
m.power_limit = pyo.Param(m.T, initialize=power_avail)
### VARS
m.dispatch = pyo.Var(
m.windows_flat, domain=pyo.Binary, doc="dispatch power in timeslot t to request r"
)
### OBJ
m.obj = pyo.Objective(
expr=sum(m.dispatch[t, r] for (t, r) in m.windows_flat), sense=pyo.maximize
)
### CONSTRAINTS
@m.Constraint(m.R)
def satisfy_only_once(m, r):
return sum(m.dispatch[t, rr] for (t, rr) in m.windows_flat if r == rr) <= 1
@m.Constraint(m.T)
def supply_limit(m, t):
# we need to sum across all dispatches that could be running in this period
possible_dispatches = {
(tt, r) for (tt, r) in m.windows_flat if 0 <= t - tt < request_duration[r]
}
if not possible_dispatches:
return pyo.Constraint.Skip
return (
sum(m.dispatch[tt, r] * m.request_size[r] for tt, r in possible_dispatches)
<= m.power_limit[t]
)
# check it
m.pprint()
# solve it
solver = pyo.SolverFactory("cbc")
res = solver.solve(m)
print(res)
m.dispatch.display()
plt.step(periods, [power_avail[p] for p in periods], color="g")
assigned_periods = {}
for t, r in m.dispatch:
if pyo.value(m.dispatch[t, r]) > 0.5: # request was assigned
assigned_periods[r] = list(range(t, t + request_duration[r]))
print("hit", t, r)
total_period_power_assigned = []
for t in m.T:
assigned = 0
for r in m.R:
if t in assigned_periods.get(r, set()):
assigned += request_power_per_period[r]
total_period_power_assigned.append(assigned)
print(total_period_power_assigned)
plt.step(periods, total_period_power_assigned)
plt.show()
Output [truncated]
dispatch : dispatch power in timeslot t to request r
Size=15, Index=windows_flat
Key : Lower : Value : Upper : Fixed : Stale : Domain
(0, 'r1') : 0 : 0.0 : 1 : False : False : Binary
(1, 'r1') : 0 : 1.0 : 1 : False : False : Binary
(1, 'r3') : 0 : 0.0 : 1 : False : False : Binary
(2, 'r1') : 0 : 0.0 : 1 : False : False : Binary
(2, 'r2') : 0 : 0.0 : 1 : False : False : Binary
(2, 'r3') : 0 : 0.0 : 1 : False : False : Binary
(3, 'r1') : 0 : 0.0 : 1 : False : False : Binary
(3, 'r2') : 0 : 1.0 : 1 : False : False : Binary
(3, 'r3') : 0 : 0.0 : 1 : False : False : Binary
(4, 'r1') : 0 : 0.0 : 1 : False : False : Binary
(4, 'r3') : 0 : 1.0 : 1 : False : False : Binary
(5, 'r1') : 0 : 0.0 : 1 : False : False : Binary
(5, 'r3') : 0 : 0.0 : 1 : False : False : Binary
(6, 'r1') : 0 : 0.0 : 1 : False : False : Binary
(7, 'r1') : 0 : 0.0 : 1 : False : False : Binary
hit 1 r1
hit 4 r3
hit 3 r2
{'r1': [1, 2, 3, 4], 'r3': [4, 5, 6, 7, 8], 'r2': [3, 4, 5, 6, 7, 8, 9, 10]}
[0, 5, 5, 9, 17, 12, 12, 12, 12, 4, 4, 0]
