Defining the value of one variable in a constraint in relation to another variable without making the problem nonlinear in Pyomo

Question

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!

Answer 1

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]