Reformat source code

pull/21/head
Alinson S. Xavier 4 years ago
parent 03268dd3df
commit 360308ef4a

@ -23,7 +23,10 @@ Example
import UnitCommitment
instance = UnitCommitment.read_benchmark("matpower/case3375wp/2017-02-01")
"""
function read_benchmark(name::AbstractString; quiet::Bool=false)::UnitCommitmentInstance
function read_benchmark(
name::AbstractString;
quiet::Bool = false,
)::UnitCommitmentInstance
basedir = dirname(@__FILE__)
filename = "$basedir/../../instances/$name.json.gz"
url = "$INSTANCES_URL/$name.json.gz"
@ -225,13 +228,17 @@ function _from_json(json; repair = true)
# Read spinning, up-flexiramp, and down-flexiramp reserve requirements
reserves = Reserves(zeros(T), zeros(T), zeros(T))
if "Reserves" in keys(json)
reserves.spinning =
timeseries(json["Reserves"]["Spinning (MW)"], default = zeros(T))
reserves.upflexiramp =
timeseries(json["Reserves"]["Up-flexiramp (MW)"], default = zeros(T))
reserves.dwflexiramp =
timeseries(json["Reserves"]["Down-flexiramp (MW)"], default = zeros(T))
end
reserves.spinning =
timeseries(json["Reserves"]["Spinning (MW)"], default = zeros(T))
reserves.upflexiramp = timeseries(
json["Reserves"]["Up-flexiramp (MW)"],
default = zeros(T),
)
reserves.dwflexiramp = timeseries(
json["Reserves"]["Down-flexiramp (MW)"],
default = zeros(T),
)
end
# Read transmission lines
if "Transmission lines" in keys(json)

@ -32,12 +32,20 @@ function build_model(;
formulation = Formulation(),
variable_names::Bool = false,
)::JuMP.Model
if formulation.ramping ==WanHob2016.Ramping() && instance.reserves.spinning!=zeros(instance.time)
error("Spinning reserves are not supported by the WanHob2016 ramping formulation")
if formulation.ramping == WanHob2016.Ramping() &&
instance.reserves.spinning != zeros(instance.time)
error(
"Spinning reserves are not supported by the WanHob2016 ramping formulation",
)
end
@show formulation.ramping
if formulation.ramping !== WanHob2016.Ramping() && (instance.reserves.upflexiramp!=zeros(instance.time) || instance.reserves.dwflexiramp!=zeros(instance.time))
error("Flexiramp is supported only by the WanHob2016 ramping formulation")
if formulation.ramping !== WanHob2016.Ramping() && (
instance.reserves.upflexiramp != zeros(instance.time) ||
instance.reserves.dwflexiramp != zeros(instance.time)
)
error(
"Flexiramp is supported only by the WanHob2016 ramping formulation",
)
end
@info "Building model..."

@ -5,12 +5,12 @@
function _add_flexiramp_vars!(model::JuMP.Model, g::Unit)::Nothing
upflexiramp = _init(model, :upflexiramp)
upflexiramp_shortfall = _init(model, :upflexiramp_shortfall)
mfg=_init(model,:mfg)
mfg = _init(model, :mfg)
dwflexiramp = _init(model, :dwflexiramp)
dwflexiramp_shortfall = _init(model, :dwflexiramp_shortfall)
for t in 1:model[:instance].time
# maximum feasible generation, \bar{g_{its}} in Wang & Hobbs (2016)
mfg[g.name,t]=@variable(model, lower_bound = 0)
mfg[g.name, t] = @variable(model, lower_bound = 0)
if g.provides_flexiramp_reserves[t]
upflexiramp[g.name, t] = @variable(model) # up-flexiramp, ur_{it} in Wang & Hobbs (2016)
dwflexiramp[g.name, t] = @variable(model) # down-flexiramp, dr_{it} in Wang & Hobbs (2016)
@ -28,8 +28,6 @@ function _add_flexiramp_vars!(model::JuMP.Model, g::Unit)::Nothing
return
end
function _add_ramp_eqs!(
model::JuMP.Model,
g::Unit,
@ -43,110 +41,145 @@ function _add_ramp_eqs!(
RU = g.ramp_up_limit
RD = g.ramp_down_limit
gn = g.name
minp=g.min_power
maxp=g.max_power
initial_power=g.initial_power
minp = g.min_power
maxp = g.max_power
initial_power = g.initial_power
is_on = model[:is_on]
prod_above = model[:prod_above]
upflexiramp=model[:upflexiramp]
dwflexiramp=model[:dwflexiramp]
mfg=model[:mfg]
upflexiramp = model[:upflexiramp]
dwflexiramp = model[:dwflexiramp]
mfg = model[:mfg]
for t in 1:model[:instance].time
@constraint(model, prod_above[gn, t] + (is_on[gn,t]*minp[t])
<=mfg[gn,t]) # Eq. (19) in Wang & Hobbs (2016)
@constraint(model, mfg[gn,t]<= is_on[gn,t]* maxp[t]) # Eq. (22) in Wang & Hobbs (2016)
if t!=model[:instance].time
@constraint(model, minp[t] * (is_on[gn,t+1]+is_on[gn,t]-1) <=
prod_above[gn, t] - dwflexiramp[gn,t] +(is_on[gn,t]*minp[t])
) # first inequality of Eq. (20) in Wang & Hobbs (2016)
@constraint(model, prod_above[gn, t] - dwflexiramp[gn,t] + (is_on[gn,t]*minp[t]) <=
mfg[gn,t+1]
+ (maxp[t] * (1-is_on[gn,t+1]))
) # second inequality of Eq. (20) in Wang & Hobbs (2016)
@constraint(model, minp[t] * (is_on[gn,t+1]+is_on[gn,t]-1) <=
prod_above[gn, t] + upflexiramp[gn,t] + (is_on[gn,t]*minp[t])
) # first inequality of Eq. (21) in Wang & Hobbs (2016)
@constraint(model, prod_above[gn, t] + upflexiramp[gn,t] +(is_on[gn,t]*minp[t]) <=
mfg[gn,t+1] + (maxp[t] * (1-is_on[gn,t+1]))
) # second inequality of Eq. (21) in Wang & Hobbs (2016)
if t!=1
@constraint(model, mfg[gn,t]<=prod_above[gn,t-1] + (is_on[gn,t-1]*minp[t])
+ (RU * is_on[gn,t-1])
+ (SU*(is_on[gn,t] - is_on[gn,t-1]))
+ maxp[t] * (1-is_on[gn,t])
) # Eq. (23) in Wang & Hobbs (2016)
@constraint(model, (prod_above[gn,t-1] + (is_on[gn,t-1]*minp[t]))
- (prod_above[gn,t] + (is_on[gn,t]*minp[t]))
<= RD * is_on[gn,t]
+ SD * (is_on[gn,t-1] - is_on[gn,t])
+ maxp[t] * (1-is_on[gn,t-1])
) # Eq. (25) in Wang & Hobbs (2016)
@constraint(
model,
prod_above[gn, t] + (is_on[gn, t] * minp[t]) <= mfg[gn, t]
) # Eq. (19) in Wang & Hobbs (2016)
@constraint(model, mfg[gn, t] <= is_on[gn, t] * maxp[t]) # Eq. (22) in Wang & Hobbs (2016)
if t != model[:instance].time
@constraint(
model,
minp[t] * (is_on[gn, t+1] + is_on[gn, t] - 1) <=
prod_above[gn, t] - dwflexiramp[gn, t] +
(is_on[gn, t] * minp[t])
) # first inequality of Eq. (20) in Wang & Hobbs (2016)
@constraint(
model,
prod_above[gn, t] - dwflexiramp[gn, t] +
(is_on[gn, t] * minp[t]) <=
mfg[gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1]))
) # second inequality of Eq. (20) in Wang & Hobbs (2016)
@constraint(
model,
minp[t] * (is_on[gn, t+1] + is_on[gn, t] - 1) <=
prod_above[gn, t] +
upflexiramp[gn, t] +
(is_on[gn, t] * minp[t])
) # first inequality of Eq. (21) in Wang & Hobbs (2016)
@constraint(
model,
prod_above[gn, t] +
upflexiramp[gn, t] +
(is_on[gn, t] * minp[t]) <=
mfg[gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1]))
) # second inequality of Eq. (21) in Wang & Hobbs (2016)
if t != 1
@constraint(
model,
mfg[gn, t] <=
prod_above[gn, t-1] +
(is_on[gn, t-1] * minp[t]) +
(RU * is_on[gn, t-1]) +
(SU * (is_on[gn, t] - is_on[gn, t-1])) +
maxp[t] * (1 - is_on[gn, t])
) # Eq. (23) in Wang & Hobbs (2016)
@constraint(
model,
(prod_above[gn, t-1] + (is_on[gn, t-1] * minp[t])) -
(prod_above[gn, t] + (is_on[gn, t] * minp[t])) <=
RD * is_on[gn, t] +
SD * (is_on[gn, t-1] - is_on[gn, t]) +
maxp[t] * (1 - is_on[gn, t-1])
) # Eq. (25) in Wang & Hobbs (2016)
else
@constraint(model, mfg[gn,t]<=initial_power
+ (RU * is_initially_on)
+ (SU*(is_on[gn,t] - is_initially_on))
+ maxp[t] * (1-is_on[gn,t])
) # Eq. (23) in Wang & Hobbs (2016) for the first time period
@constraint(model, initial_power
- (prod_above[gn,t] + (is_on[gn,t]*minp[t]))
<= RD * is_on[gn,t]
+ SD * (is_initially_on - is_on[gn,t])
+ maxp[t] * (1-is_initially_on)
) # Eq. (25) in Wang & Hobbs (2016) for the first time period
@constraint(
model,
mfg[gn, t] <=
initial_power +
(RU * is_initially_on) +
(SU * (is_on[gn, t] - is_initially_on)) +
maxp[t] * (1 - is_on[gn, t])
) # Eq. (23) in Wang & Hobbs (2016) for the first time period
@constraint(
model,
initial_power -
(prod_above[gn, t] + (is_on[gn, t] * minp[t])) <=
RD * is_on[gn, t] +
SD * (is_initially_on - is_on[gn, t]) +
maxp[t] * (1 - is_initially_on)
) # Eq. (25) in Wang & Hobbs (2016) for the first time period
end
@constraint(model, mfg[gn,t]<=
(SD*(is_on[gn,t] - is_on[gn,t+1]))
+ (maxp[t] * is_on[gn,t+1])
) # Eq. (24) in Wang & Hobbs (2016)
@constraint(model, -RD * is_on[gn,t+1]
-SD * (is_on[gn,t]-is_on[gn,t+1])
-maxp[t] * (1-is_on[gn,t])
<= upflexiramp[gn,t]
) # first inequality of Eq. (26) in Wang & Hobbs (2016)
@constraint(model, upflexiramp[gn,t] <=
RU * is_on[gn,t]
+ SU * (is_on[gn,t+1]-is_on[gn,t])
+ maxp[t] * (1-is_on[gn,t+1])
) # second inequality of Eq. (26) in Wang & Hobbs (2016)
@constraint(model, -RU * is_on[gn,t]
-SU * (is_on[gn,t+1]-is_on[gn,t])
-maxp[t] * (1-is_on[gn,t+1])
<= dwflexiramp[gn,t]
) # first inequality of Eq. (27) in Wang & Hobbs (2016)
@constraint(model, dwflexiramp[gn,t] <=
RD * is_on[gn,t+1]
+ SD * (is_on[gn,t]-is_on[gn,t+1])
+ maxp[t] * (1-is_on[gn,t])
) # second inequality of Eq. (27) in Wang & Hobbs (2016)
@constraint(model, -maxp[t] * is_on[gn,t]
+minp[t] * is_on[gn,t+1]
<= upflexiramp[gn,t]
) # first inequality of Eq. (28) in Wang & Hobbs (2016)
@constraint(model, upflexiramp[gn,t] <=
maxp[t] * is_on[gn,t+1]
) # second inequality of Eq. (28) in Wang & Hobbs (2016)
@constraint(model, -maxp[t] * is_on[gn,t+1]
<= dwflexiramp[gn,t]
) # first inequality of Eq. (29) in Wang & Hobbs (2016)
@constraint(model, dwflexiramp[gn,t] <=
(maxp[t] * is_on[gn,t])
-(minp[t] * is_on[gn,t+1])
) # second inequality of Eq. (29) in Wang & Hobbs (2016)
@constraint(
model,
mfg[gn, t] <=
(SD * (is_on[gn, t] - is_on[gn, t+1])) +
(maxp[t] * is_on[gn, t+1])
) # Eq. (24) in Wang & Hobbs (2016)
@constraint(
model,
-RD * is_on[gn, t+1] - SD * (is_on[gn, t] - is_on[gn, t+1]) -
maxp[t] * (1 - is_on[gn, t]) <= upflexiramp[gn, t]
) # first inequality of Eq. (26) in Wang & Hobbs (2016)
@constraint(
model,
upflexiramp[gn, t] <=
RU * is_on[gn, t] +
SU * (is_on[gn, t+1] - is_on[gn, t]) +
maxp[t] * (1 - is_on[gn, t+1])
) # second inequality of Eq. (26) in Wang & Hobbs (2016)
@constraint(
model,
-RU * is_on[gn, t] - SU * (is_on[gn, t+1] - is_on[gn, t]) -
maxp[t] * (1 - is_on[gn, t+1]) <= dwflexiramp[gn, t]
) # first inequality of Eq. (27) in Wang & Hobbs (2016)
@constraint(
model,
dwflexiramp[gn, t] <=
RD * is_on[gn, t+1] +
SD * (is_on[gn, t] - is_on[gn, t+1]) +
maxp[t] * (1 - is_on[gn, t])
) # second inequality of Eq. (27) in Wang & Hobbs (2016)
@constraint(
model,
-maxp[t] * is_on[gn, t] + minp[t] * is_on[gn, t+1] <=
upflexiramp[gn, t]
) # first inequality of Eq. (28) in Wang & Hobbs (2016)
@constraint(model, upflexiramp[gn, t] <= maxp[t] * is_on[gn, t+1]) # second inequality of Eq. (28) in Wang & Hobbs (2016)
@constraint(model, -maxp[t] * is_on[gn, t+1] <= dwflexiramp[gn, t]) # first inequality of Eq. (29) in Wang & Hobbs (2016)
@constraint(
model,
dwflexiramp[gn, t] <=
(maxp[t] * is_on[gn, t]) - (minp[t] * is_on[gn, t+1])
) # second inequality of Eq. (29) in Wang & Hobbs (2016)
else
@constraint(model, mfg[gn,t]<=prod_above[gn,t-1] + (is_on[gn,t-1]*minp[t])
+ (RU * is_on[gn,t-1])
+ (SU*(is_on[gn,t] - is_on[gn,t-1]))
+ maxp[t] * (1-is_on[gn,t])
) # Eq. (23) in Wang & Hobbs (2016) for the last time period
@constraint(model, (prod_above[gn,t-1] + (is_on[gn,t-1]*minp[t]))
- (prod_above[gn,t] + (is_on[gn,t]*minp[t]))
<= RD * is_on[gn,t]
+ SD * (is_on[gn,t-1] - is_on[gn,t])
+ maxp[t] * (1-is_on[gn,t-1])
) # Eq. (25) in Wang & Hobbs (2016) for the last time period
@constraint(
model,
mfg[gn, t] <=
prod_above[gn, t-1] +
(is_on[gn, t-1] * minp[t]) +
(RU * is_on[gn, t-1]) +
(SU * (is_on[gn, t] - is_on[gn, t-1])) +
maxp[t] * (1 - is_on[gn, t])
) # Eq. (23) in Wang & Hobbs (2016) for the last time period
@constraint(
model,
(prod_above[gn, t-1] + (is_on[gn, t-1] * minp[t])) -
(prod_above[gn, t] + (is_on[gn, t] * minp[t])) <=
RD * is_on[gn, t] +
SD * (is_on[gn, t-1] - is_on[gn, t]) +
maxp[t] * (1 - is_on[gn, t-1])
) # Eq. (25) in Wang & Hobbs (2016) for the last time period
end
end
end
end

@ -12,7 +12,6 @@ module WanHob2016
import ..RampingFormulation
struct Ramping <: RampingFormulation end
end

@ -57,7 +57,7 @@ function _add_reserve_eqs!(model::JuMP.Model)::Nothing
end
function _add_flexiramp_eqs!(model::JuMP.Model)::Nothing
# Note: The flexpramp requirements in Wang & Hobbs (2016) are imposed as hard constraints
# Note: The flexpramp requirements in Wang & Hobbs (2016) are imposed as hard constraints
# through Eq. (17) and Eq. (18). The constraints eq_min_upflexiramp[t] and eq_min_dwflexiramp[t]
# provided below are modified versions of Eq. (17) and Eq. (18), respectively, in that
# they include slack variables for flexiramp shortfall, which are penalized in the
@ -71,15 +71,19 @@ function _add_flexiramp_eqs!(model::JuMP.Model)::Nothing
eq_min_upflexiramp[t] = @constraint(
model,
sum(model[:upflexiramp][g.name, t] for g in instance.units) +
(flexiramp_shortfall_penalty >= 0 ? model[:upflexiramp_shortfall][t] : 0.0) >=
instance.reserves.upflexiramp[t]
(
flexiramp_shortfall_penalty >= 0 ?
model[:upflexiramp_shortfall][t] : 0.0
) >= instance.reserves.upflexiramp[t]
)
# Eq. (18) in Wang & Hobbs (2016)
eq_min_dwflexiramp[t] = @constraint(
model,
sum(model[:dwflexiramp][g.name, t] for g in instance.units) +
(flexiramp_shortfall_penalty >= 0 ? model[:dwflexiramp_shortfall][t] : 0.0) >=
instance.reserves.dwflexiramp[t]
(
flexiramp_shortfall_penalty >= 0 ?
model[:dwflexiramp_shortfall][t] : 0.0
) >= instance.reserves.dwflexiramp[t]
)
# Account for flexiramp shortfall contribution to objective
@ -87,7 +91,10 @@ function _add_flexiramp_eqs!(model::JuMP.Model)::Nothing
add_to_expression!(
model[:obj],
flexiramp_shortfall_penalty,
(model[:upflexiramp_shortfall][t]+model[:dwflexiramp_shortfall][t]),
(
model[:upflexiramp_shortfall][t] +
model[:dwflexiramp_shortfall][t]
),
)
end
end

@ -50,22 +50,25 @@ function solution(model::JuMP.Model)::OrderedDict
sol["Is on"] = timeseries(model[:is_on], instance.units)
sol["Switch on"] = timeseries(model[:switch_on], instance.units)
sol["Switch off"] = timeseries(model[:switch_off], instance.units)
if instance.reserves.upflexiramp != zeros(T) || instance.reserves.dwflexiramp != zeros(T)
if instance.reserves.upflexiramp != zeros(T) ||
instance.reserves.dwflexiramp != zeros(T)
# Report flexiramp solutions only if either of the up-flexiramp and
# down-flexiramp requirements is not a default array of zeros
sol["Up-flexiramp (MW)"] = timeseries(model[:upflexiramp], instance.units)
sol["Up-flexiramp (MW)"] =
timeseries(model[:upflexiramp], instance.units)
sol["Up-flexiramp shortfall (MW)"] = OrderedDict(
t =>
(instance.flexiramp_shortfall_penalty[t] >= 0) ?
round(value(model[:upflexiramp_shortfall][t]), digits = 5) : 0.0 for
t in 1:instance.time
round(value(model[:upflexiramp_shortfall][t]), digits = 5) :
0.0 for t in 1:instance.time
)
sol["Down-flexiramp (MW)"] = timeseries(model[:dwflexiramp], instance.units)
sol["Down-flexiramp (MW)"] =
timeseries(model[:dwflexiramp], instance.units)
sol["Down-flexiramp shortfall (MW)"] = OrderedDict(
t =>
(instance.flexiramp_shortfall_penalty[t] >= 0) ?
round(value(model[:dwflexiramp_shortfall][t]), digits = 5) : 0.0 for
t in 1:instance.time
round(value(model[:dwflexiramp_shortfall][t]), digits = 5) :
0.0 for t in 1:instance.time
)
else
# Report spinning reserve solutions only if both up-flexiramp and
@ -74,10 +77,9 @@ function solution(model::JuMP.Model)::OrderedDict
sol["Reserve shortfall (MW)"] = OrderedDict(
t =>
(instance.shortfall_penalty[t] >= 0) ?
round(value(model[:reserve_shortfall][t]), digits = 5) : 0.0 for
t in 1:instance.time
)
round(value(model[:reserve_shortfall][t]), digits = 5) :
0.0 for t in 1:instance.time
)
end
sol["Net injection (MW)"] =
timeseries(model[:net_injection], instance.buses)

@ -339,13 +339,15 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
err_count += 1
end
upflexiramp =
sum(solution["Up-flexiramp (MW)"][g.name][t] for g in instance.units)
upflexiramp = sum(
solution["Up-flexiramp (MW)"][g.name][t] for g in instance.units
)
upflexiramp_shortfall =
(instance.flexiramp_shortfall_penalty[t] >= 0) ?
solution["Up-flexiramp shortfall (MW)"][t] : 0
if upflexiramp + upflexiramp_shortfall < instance.reserves.upflexiramp[t] - tol
if upflexiramp + upflexiramp_shortfall <
instance.reserves.upflexiramp[t] - tol
@error @sprintf(
"Insufficient up-flexiramp at time %d (%.2f + %.2f should be %.2f)",
t,
@ -356,13 +358,15 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
err_count += 1
end
dwflexiramp =
sum(solution["Down-flexiramp (MW)"][g.name][t] for g in instance.units)
dwflexiramp = sum(
solution["Down-flexiramp (MW)"][g.name][t] for g in instance.units
)
dwflexiramp_shortfall =
(instance.flexiramp_shortfall_penalty[t] >= 0) ?
solution["Down-flexiramp shortfall (MW)"][t] : 0
if dwflexiramp + dwflexiramp_shortfall < instance.reserves.dwflexiramp[t] - tol
if dwflexiramp + dwflexiramp_shortfall <
instance.reserves.dwflexiramp[t] - tol
@error @sprintf(
"Insufficient down-flexiramp at time %d (%.2f + %.2f should be %.2f)",
t,

Loading…
Cancel
Save