I currently try to use this unit commitment example to build my own model with pyomo. After defining switch-on and switch-off variables I struggle to implement the following equation:Equation
The yalmip example is pretty straight forward:
for k = 2:Horizon
for unit = 1:Nunits
% indicator will be 1 only when switched on
indicator = onoff(unit,k)-onoff(unit,k-1);
range = k:min(Horizon,k+minup(unit)-1);
% Constraints will be redundant unless indicator = 1
Constraints = [Constraints, onoff(unit,range) >= indicator];
end
end
Right now I am only looking into one unit, which gives me this model.
model = ConcreteModel()
p = prices
ts = timesteps
ut = min_uptime1
model.x = Var(ts, within = Binary) #onoff
model.v = Var(ts, within = Binary) #switch_on
model.w = Var(ts, within = Binary) #switch_off
def obj_rule(model):
return sum(p[t] * model.x[t] - 0.001 * (model.v[t] + model.w[t]) for t in ts)
model.revenue = Objective(rule = obj_rule, sense = maximize)
#start-up, shut-down costs will be added
def daily_uptime_rule (model):
return sum(model.x[t] for t in ts) == 12
model.daily_uptime_rule = \
Constraint(rule = daily_uptime_rule)
def switch_on(model, t):
if t == ts[0]:
return model.v[t] >= 1 - (1 - model.x[t])
else:
return model.v[t] >= 1 - model.x[t-1] - (1 - model.x[t])
model.switch_on = \
Constraint(ts, rule = switch_on)
def switch_off(model, t):
if t == ts[23]:
return model.w[t] >= model.x[t]
else:
return model.w[t] >= 1 - model.x[t+1] + (model.x[t] - 1)
model.switch_off = \
Constraint(ts, rule = switch_off)
def min_ut(model, t):
a = list(range(t, (min(ts[23], t+ut-1)+1)))
for i in a:
return model.x[i] >= model.v[t]
model.min_ut = \
Constraint(ts, rule = min_ut)
My problem here is, that i can't access the variable x the same way in pyomo. For every timestep t we need constraints for t+1, t+2, .. t+min_up -1, but I can't use ranges with variables (model.x). Can I use the yalmip example in pyomo or do i need a new formulation?
Ok, so it seems the fundamental issue here is that the index of summation that you would like to do is dependent on the RHS of the inequality. You can construct the indices of the summation in a couple ways. You just need to be careful that the values you construct are valid. Here is an idea that might help you. This toy model tries to maximize the sum of x[t], but limits x[t] <= x[t-1] + x[t-2] just for giggles. Note the construction of the summation range "on the fly" from the passed value of t:
from pyomo.environ import *
m = ConcreteModel()
m.t = Set(initialize=range(5))
m.x = Var(m.t)
# constrain x_t to be less than the sum of x_(t-1), x_(t-2)
def x_limiter(m, t):
if t==0:
return m.x[t] <= 1 # limit the first value
# upperlimit on remainder is sum of previous 2
return sum(m.x[i] for i in range(t-2, t) if i in m.t) >= m.x[t]
m.c1 = Constraint(m.t, rule=x_limiter)
# try to maximize x
m.OBJ = Objective(expr=sum(m.x[t] for t in m.t), sense=maximize)
solver = SolverFactory('glpk')
m.pprint()
solver.solve(m)
m.display()
It gets the job done:
Variables:
x : Size=5, Index=t
Key : Lower : Value : Upper : Fixed : Stale : Domain
0 : None : 1.0 : None : False : False : Reals
1 : None : 1.0 : None : False : False : Reals
2 : None : 2.0 : None : False : False : Reals
3 : None : 3.0 : None : False : False : Reals
4 : None : 5.0 : None : False : False : Reals
Objectives:
OBJ : Size=1, Index=None, Active=True
Key : Active : Value
None : True : 12.0
This recent post also has a similar idea: