I am working on an optimisation problem, where there is a device that is shiftable and can be scheduled for other hours rather than its initial schedule. The problem is that this device has a specific duty cycle, I mean if it is ON, it has to be ON for at least T_ON hours and cannot be on for more than T_ON hours as well.
There is a binary variable, I(t) showing ON (1) or OFF (0) status of the device. In Pyomo, as far as I know, we cannot define an if statement on variables, so I am quite stuck in how to formulise this problem in Pyomo. I also tried to make a graphic version of what I am looking for in figure below.
I really appreciate it if anyone can help me with an idea on how to formulise this problem.
Cheers! :)
If the "if" statement was working in Pyomo, I would write it as below:
model.Gen = Set(initialize=RangeSet(0,End_uses -1), doc='set of enduses') #index of end uses
model.T = Set(initialize=RangeSet(0,50-1), doc='Time') #index of data samples
model.taw = SetOf(model.T) # references whatever is in T
#Yes/No of each item (binary variable) --> include or exclude an item
model.I_test=Var(model.Gen,model.T,within=Binary)
#Min Up time
def MinUp(model,g,t):
if t > 0:
if model.I_test[g,t]-model.I_test[g,t-1]>0:
return sum(model.I_test[g,t2] for t2 in model.taw if t2>=t and t2<=t+T_on[g]) >= T_on[g]*(model.I_test[g,t]-model.I_test[g,t-1])
else:
return Constraint.Skip
model.cons5 = Constraint(model.Gen,model.T,rule=MinUp)
I am looking for a help on formulising the problem or even a similar piece of code that can help me to solve this issue.
The piece you are missing is an additional variable to keep track of which periods the machine was started in. This will allow several other types of constraints to be drawn using some integer arithmetic. Below is an example of several.
It isn't super clear from your problem which of these would be most useful to you. You seem to imply that the min_runtime == max_runtime and then you mention a duty cycle, but don't define it.
If your specification simplifies down to exactly n hours on followed by minimum of n hours off, you could re-write one of these types of constraints below such that the n periods after a start are "on" and then n+1 to 2n periods are off and delete the other constraints.
# cycling machine example
import pyomo.environ as pyo
# machine characteristics
min_runtime = 2 # minimum number of periods to run after startup
min_reset_time = 3 # off periods after being turned off
duty_cycle = (4, 8) # 50% duty cycle evaluated over 8 periods "4 of 8"
time_periods = 12
demand = 7 # periods of 'runtime' from the machine
m = pyo.ConcreteModel()
# SETS
m.T = pyo.Set(initialize=range(time_periods))
# VARS
m.start = pyo.Var(m.T, domain=pyo.Binary, doc='start machine')
m.running = pyo.Var(m.T, domain=pyo.Binary, doc='machine running')
# OBJ
m.obj = pyo.Objective(expr=sum(m.running[t] for t in m.T)) # minimze the runtime
# CONSTRAINTS
# satisfy the demand
m.satisfy_demand = pyo.Constraint(expr=sum(m.running[t] for t in m.T) >= demand)
@m.Constraint(m.T)
def link_running(m, t):
"""machine can only be running if it was running in prior period or started in this one"""
if t == m.T.first():
return pyo.Constraint.Skip
return m.running[t] <= m.running[t-1] + m.start[t]
@m.Constraint(m.T)
def min_runtime(m, t):
"""force the minimum runtime after a start event"""
# make a subset of the time periods of interest...
next_time_periods = {t + offset for offset in range(min_runtime) if t + offset < time_periods}
return sum(m.running[tt] for tt in next_time_periods) >= len(next_time_periods) * m.start[t]
@m.Constraint(m.T)
def honor_reset_time(m, t):
if t == m.T.first():
return pyo.Constraint.Skip
previous_time_periods = {t - offset for offset in range(1, min_reset_time + 1) if t - offset >=0}
return len(previous_time_periods) * m.start[t] <= sum(1-m.running[tt] for tt in previous_time_periods)
@m.Constraint(m.T)
def duty_cycle_limit(m, t):
""" evaluate all possible duty cycles and impose limit """
if t + duty_cycle[1] > time_periods:
return pyo.Constraint.Skip
duty_cycle_periods = list(range(t, t + duty_cycle[1]))
return sum(m.running[tt] for tt in duty_cycle_periods) <= duty_cycle[0]
m.pprint()
solver = pyo.SolverFactory('glpk')
soln = solver.solve(m)
print(soln)
m.running.display()
Solver:
- Status: ok
Termination condition: optimal
Statistics:
Branch and bound:
Number of bounded subproblems: 1
Number of created subproblems: 1
Error rc: 0
Time: 0.0057790279388427734
Solution:
- number of solutions: 0
number of solutions displayed: 0
running : machine running
Size=12, Index=T
Key : Lower : Value : Upper : Fixed : Stale : Domain
0 : 0 : 1.0 : 1 : False : False : Binary
1 : 0 : 1.0 : 1 : False : False : Binary
2 : 0 : 1.0 : 1 : False : False : Binary
3 : 0 : 1.0 : 1 : False : False : Binary
4 : 0 : 0.0 : 1 : False : False : Binary
5 : 0 : 0.0 : 1 : False : False : Binary
6 : 0 : 0.0 : 1 : False : False : Binary
7 : 0 : 0.0 : 1 : False : False : Binary
8 : 0 : 1.0 : 1 : False : False : Binary
9 : 0 : 1.0 : 1 : False : False : Binary
10 : 0 : 1.0 : 1 : False : False : Binary
11 : 0 : 0.0 : 1 : False : False : Binary