new formatting

pull/21/head
oyurdakul 4 years ago
parent febb4f1aad
commit b4bc50c865

@ -129,15 +129,19 @@ formulations = Dict(
const gap_limit = parse(Float64, args["--gap"])
const time_limit = parse(Float64, args["--time-limit"])
methods = Dict(
"default" =>
XavQiuWanThi2019.Method(time_limit = time_limit, gap_limit = gap_limit),
"default" => XavQiuWanThi2019.Method(
time_limit = time_limit,
gap_limit = gap_limit,
),
)
# MIP solvers
# -----------------------------------------------------------------------------
optimizers = Dict(
"gurobi" =>
optimizer_with_attributes(Gurobi.Optimizer, "Threads" => Threads.nthreads()),
"gurobi" => optimizer_with_attributes(
Gurobi.Optimizer,
"Threads" => Threads.nthreads(),
),
)
# Parse command line arguments

@ -41,7 +41,7 @@ function read_egret_solution(path::String)::OrderedDict
startup_cost[gen_name] = zeros(T)
production_cost[gen_name] = zeros(T)
if "commitment_cost" in keys(gen_dict)
for t = 1:T
for t in 1:T
x = gen_dict["commitment"]["values"][t]
commitment_cost = gen_dict["commitment_cost"]["values"][t]
prod_above_cost = gen_dict["production_cost"]["values"][t]

@ -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"
@ -62,7 +65,9 @@ function read(path::AbstractString)::UnitCommitmentInstance
end
function _read(file::IO)::UnitCommitmentInstance
return _from_json(JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)))
return _from_json(
JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)),
)
end
function _read_json(path::String)::OrderedDict
@ -92,7 +97,8 @@ function _from_json(json; repair = true)
end
time_horizon !== nothing || error("Missing parameter: Time horizon (h)")
time_step = scalar(json["Parameters"]["Time step (min)"], default = 60)
(60 % time_step == 0) || error("Time step $time_step is not a divisor of 60")
(60 % time_step == 0) ||
error("Time step $time_step is not a divisor of 60")
time_multiplier = 60 ÷ time_step
T = time_horizon * time_multiplier
@ -102,23 +108,23 @@ function _from_json(json; repair = true)
function timeseries(x; default = nothing)
x !== nothing || return default
x isa Array || return [x for t = 1:T]
x isa Array || return [x for t in 1:T]
return x
end
# Read parameters
power_balance_penalty = timeseries(
json["Parameters"]["Power balance penalty (\$/MW)"],
default = [1000.0 for t = 1:T],
default = [1000.0 for t in 1:T],
)
# Penalty price for shortage in meeting system-wide flexiramp requirements
flexiramp_shortfall_penalty = timeseries(
json["Parameters"]["Flexiramp penalty (\$/MW)"],
default = [500.0 for t = 1:T],
default = [500.0 for t in 1:T],
)
shortfall_penalty = timeseries(
json["Parameters"]["Reserve shortfall penalty (\$/MW)"],
default = [-1.0 for t = 1:T],
default = [-1.0 for t in 1:T],
)
# Read buses
@ -140,14 +146,17 @@ function _from_json(json; repair = true)
# Read production cost curve
K = length(dict["Production cost curve (MW)"])
curve_mw = hcat([timeseries(dict["Production cost curve (MW)"][k]) for k = 1:K]...)
curve_cost =
hcat([timeseries(dict["Production cost curve (\$)"][k]) for k = 1:K]...)
curve_mw = hcat(
[timeseries(dict["Production cost curve (MW)"][k]) for k in 1:K]...,
)
curve_cost = hcat(
[timeseries(dict["Production cost curve (\$)"][k]) for k in 1:K]...,
)
min_power = curve_mw[:, 1]
max_power = curve_mw[:, K]
min_power_cost = curve_cost[:, 1]
segments = CostSegment[]
for k = 2:K
for k in 2:K
amount = curve_mw[:, k] - curve_mw[:, k-1]
cost = (curve_cost[:, k] - curve_cost[:, k-1]) ./ amount
replace!(cost, NaN => 0.0)
@ -158,10 +167,13 @@ function _from_json(json; repair = true)
startup_delays = scalar(dict["Startup delays (h)"], default = [1])
startup_costs = scalar(dict["Startup costs (\$)"], default = [0.0])
startup_categories = StartupCategory[]
for k = 1:length(startup_delays)
for k in 1:length(startup_delays)
push!(
startup_categories,
StartupCategory(startup_delays[k] .* time_multiplier, startup_costs[k]),
StartupCategory(
startup_delays[k] .* time_multiplier,
startup_costs[k],
),
)
end
@ -174,7 +186,8 @@ function _from_json(json; repair = true)
else
initial_status !== nothing ||
error("unit $unit_name has initial power but no initial status")
initial_status != 0 || error("unit $unit_name has invalid initial status")
initial_status != 0 ||
error("unit $unit_name has invalid initial status")
if initial_status < 0 && initial_power > 1e-3
error("unit $unit_name has invalid initial power")
end
@ -186,7 +199,7 @@ function _from_json(json; repair = true)
bus,
max_power,
min_power,
timeseries(dict["Must run?"], default = [false for t = 1:T]),
timeseries(dict["Must run?"], default = [false for t in 1:T]),
min_power_cost,
segments,
scalar(dict["Minimum uptime (h)"], default = 1) * time_multiplier,
@ -197,8 +210,14 @@ function _from_json(json; repair = true)
scalar(dict["Shutdown limit (MW)"], default = 1e6),
initial_status,
initial_power,
timeseries(dict["Provides spinning reserves?"], default = [true for t = 1:T]),
timeseries(dict["Provides flexible capacity?"], default = [true for t = 1:T]),
timeseries(
dict["Provides spinning reserves?"],
default = [true for t in 1:T],
),
timeseries(
dict["Provides flexible capacity?"],
default = [true for t in 1:T],
),
startup_categories,
)
push!(bus.units, unit)
@ -211,10 +230,14 @@ function _from_json(json; repair = true)
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))
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
@ -227,11 +250,17 @@ function _from_json(json; repair = true)
name_to_bus[dict["Target bus"]],
scalar(dict["Reactance (ohms)"]),
scalar(dict["Susceptance (S)"]),
timeseries(dict["Normal flow limit (MW)"], default = [1e8 for t = 1:T]),
timeseries(dict["Emergency flow limit (MW)"], default = [1e8 for t = 1:T]),
timeseries(
dict["Normal flow limit (MW)"],
default = [1e8 for t in 1:T],
),
timeseries(
dict["Emergency flow limit (MW)"],
default = [1e8 for t in 1:T],
),
timeseries(
dict["Flow limit penalty (\$/MW)"],
default = [5000.0 for t = 1:T],
default = [5000.0 for t in 1:T],
),
)
name_to_line[line_name] = line
@ -245,10 +274,12 @@ function _from_json(json; repair = true)
affected_units = Unit[]
affected_lines = TransmissionLine[]
if "Affected lines" in keys(dict)
affected_lines = [name_to_line[l] for l in dict["Affected lines"]]
affected_lines =
[name_to_line[l] for l in dict["Affected lines"]]
end
if "Affected units" in keys(dict)
affected_units = [name_to_unit[u] for u in dict["Affected units"]]
affected_units =
[name_to_unit[u] for u in dict["Affected units"]]
end
cont = Contingency(cont_name, affected_lines, affected_units)
push!(contingencies, cont)

@ -96,7 +96,10 @@ function Base.show(io::IO, instance::UnitCommitmentInstance)
print(io, "$(length(instance.buses)) buses, ")
print(io, "$(length(instance.lines)) lines, ")
print(io, "$(length(instance.contingencies)) contingencies, ")
print(io, "$(length(instance.price_sensitive_loads)) price sensitive loads, ")
print(
io,
"$(length(instance.price_sensitive_loads)) price sensitive loads, ",
)
print(io, "$(instance.time) time steps")
print(io, ")")
return

@ -34,14 +34,18 @@ function build_model(;
)::JuMP.Model
if formulation.ramping == WanHob2016.Ramping() &&
instance.reserves.spinning >= ones(instance.time) .* 1e-6
error("Spinning reserves are not supported by the WanHob2016 ramping formulation")
error(
"Spinning reserves are not supported by the WanHob2016 ramping formulation",
)
end
if formulation.ramping !== WanHob2016.Ramping() && (
instance.reserves.upflexiramp >= ones(instance.time) .* 1e-6 ||
instance.reserves.dwflexiramp >= ones(instance.time) .* 1e-6
)
error("Flexiramp is supported only by the WanHob2016 ramping formulation")
error(
"Flexiramp is supported only by the WanHob2016 ramping formulation",
)
end
@info "Building model..."

@ -32,7 +32,7 @@ function _add_ramp_eqs!(
switch_off = model[:switch_off]
switch_on = model[:switch_on]
for t = 1:model[:instance].time
for t in 1:model[:instance].time
# Ramp up limit
if t == 1
if is_initially_on
@ -49,8 +49,12 @@ function _add_ramp_eqs!(
max_prod_this_period =
g.min_power[t] * is_on[gn, t] +
prod_above[gn, t] +
(RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? reserve[gn, t] : 0.0)
min_prod_last_period = g.min_power[t-1] * is_on[gn, t-1] + prod_above[gn, t-1]
(
RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ?
reserve[gn, t] : 0.0
)
min_prod_last_period =
g.min_power[t-1] * is_on[gn, t-1] + prod_above[gn, t-1]
# Equation (24) in Kneuven et al. (2020)
eq_ramp_up[gn, t] = @constraint(
@ -77,10 +81,11 @@ function _add_ramp_eqs!(
g.min_power[t-1] * is_on[gn, t-1] +
prod_above[gn, t-1] +
(
RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ? reserve[gn, t-1] :
0.0
RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ?
reserve[gn, t-1] : 0.0
)
min_prod_this_period = g.min_power[t] * is_on[gn, t] + prod_above[gn, t]
min_prod_this_period =
g.min_power[t] * is_on[gn, t] + prod_above[gn, t]
# Equation (25) in Kneuven et al. (2020)
eq_ramp_down[gn, t] = @constraint(

@ -18,16 +18,18 @@ function _add_production_piecewise_linear_eqs!(
prod_above = model[:prod_above]
K = length(g.cost_segments)
for t = 1:model[:instance].time
for t in 1:model[:instance].time
gn = g.name
for k = 1:K
for k in 1:K
# Equation (45) in Kneuven et al. (2020)
# NB: when reading instance, UnitCommitment.jl already calculates
# difference between max power for segments k and k-1 so the
# value of cost_segments[k].mw[t] is the max production *for
# that segment*
eq_segprod_limit[gn, t, k] =
@constraint(model, segprod[gn, t, k] <= g.cost_segments[k].mw[t])
eq_segprod_limit[gn, t, k] = @constraint(
model,
segprod[gn, t, k] <= g.cost_segments[k].mw[t]
)
# Also add this as an explicit upper bound on segprod to make the
# solver's work a bit easier
@ -35,12 +37,18 @@ function _add_production_piecewise_linear_eqs!(
# Definition of production
# Equation (43) in Kneuven et al. (2020)
eq_prod_above_def[gn, t] =
@constraint(model, prod_above[gn, t] == sum(segprod[gn, t, k] for k = 1:K))
eq_prod_above_def[gn, t] = @constraint(
model,
prod_above[gn, t] == sum(segprod[gn, t, k] for k in 1:K)
)
# Objective function
# Equation (44) in Kneuven et al. (2020)
add_to_expression!(model[:obj], segprod[gn, t, k], g.cost_segments[k].cost[t])
add_to_expression!(
model[:obj],
segprod[gn, t, k],
g.cost_segments[k].cost[t],
)
end
end
end

@ -33,8 +33,9 @@ function _add_ramp_eqs!(
switch_off = model[:switch_off]
switch_on = model[:switch_on]
for t = 1:model[:instance].time
time_invariant = (t > 1) ? (abs(g.min_power[t] - g.min_power[t-1]) < 1e-7) : true
for t in 1:model[:instance].time
time_invariant =
(t > 1) ? (abs(g.min_power[t] - g.min_power[t-1]) < 1e-7) : true
# if t > 1 && !time_invariant
# @warn(
@ -47,8 +48,10 @@ function _add_ramp_eqs!(
# end
max_prod_this_period =
prod_above[gn, t] +
(RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? reserve[gn, t] : 0.0)
prod_above[gn, t] + (
RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ?
reserve[gn, t] : 0.0
)
min_prod_last_period = 0.0
if t > 1 && time_invariant
min_prod_last_period = prod_above[gn, t-1]
@ -58,7 +61,8 @@ function _add_ramp_eqs!(
eq_str_ramp_up[gn, t] = @constraint(
model,
max_prod_this_period - min_prod_last_period <=
(SU - g.min_power[t] - RU) * switch_on[gn, t] + RU * is_on[gn, t]
(SU - g.min_power[t] - RU) * switch_on[gn, t] +
RU * is_on[gn, t]
)
elseif (t == 1 && is_initially_on) || (t > 1 && !time_invariant)
if t > 1
@ -99,7 +103,8 @@ function _add_ramp_eqs!(
eq_str_ramp_down[gn, t] = @constraint(
model,
max_prod_last_period - min_prod_this_period <=
(SD - g.min_power[t] - RD) * switch_off[gn, t] + RD * on_last_period
(SD - g.min_power[t] - RD) * switch_off[gn, t] +
RD * on_last_period
)
elseif (t == 1 && is_initially_on) || (t > 1 && !time_invariant)
# Add back in min power

@ -9,8 +9,8 @@ function _add_production_vars!(
)::Nothing
prod_above = _init(model, :prod_above)
segprod = _init(model, :segprod)
for t = 1:model[:instance].time
for k = 1:length(g.cost_segments)
for t in 1:model[:instance].time
for k in 1:length(g.cost_segments)
segprod[g.name, t, k] = @variable(model, lower_bound = 0)
end
prod_above[g.name, t] = @variable(model, lower_bound = 0)
@ -28,7 +28,7 @@ function _add_production_limit_eqs!(
prod_above = model[:prod_above]
reserve = model[:reserve]
gn = g.name
for t = 1:model[:instance].time
for t in 1:model[:instance].time
# Objective function terms for production costs
# Part of (69) of Kneuven et al. (2020) as C^R_g * u_g(t) term
add_to_expression!(model[:obj], is_on[gn, t], g.min_power_cost[t])

@ -21,13 +21,15 @@ function _add_production_piecewise_linear_eqs!(
is_on = model[:is_on]
K = length(g.cost_segments)
for t = 1:model[:instance].time
for t in 1:model[:instance].time
# Definition of production
# Equation (43) in Kneuven et al. (2020)
eq_prod_above_def[gn, t] =
@constraint(model, prod_above[gn, t] == sum(segprod[gn, t, k] for k = 1:K))
eq_prod_above_def[gn, t] = @constraint(
model,
prod_above[gn, t] == sum(segprod[gn, t, k] for k in 1:K)
)
for k = 1:K
for k in 1:K
# Equation (42) in Kneuven et al. (2020)
# Without this, solvers will add a lot of implied bound cuts to
# have this same effect.
@ -46,7 +48,11 @@ function _add_production_piecewise_linear_eqs!(
# Objective function
# Equation (44) in Kneuven et al. (2020)
add_to_expression!(model[:obj], segprod[gn, t, k], g.cost_segments[k].cost[t])
add_to_expression!(
model[:obj],
segprod[gn, t, k],
g.cost_segments[k].cost[t],
)
end
end
return

@ -10,7 +10,7 @@ function _add_status_vars!(
is_on = _init(model, :is_on)
switch_on = _init(model, :switch_on)
switch_off = _init(model, :switch_off)
for t = 1:model[:instance].time
for t in 1:model[:instance].time
if g.must_run[t]
is_on[g.name, t] = 1.0
switch_on[g.name, t] = (t == 1 ? 1.0 - _is_initially_on(g) : 0.0)
@ -34,7 +34,7 @@ function _add_status_eqs!(
is_on = model[:is_on]
switch_off = model[:switch_off]
switch_on = model[:switch_on]
for t = 1:model[:instance].time
for t in 1:model[:instance].time
if !g.must_run[t]
# Link binary variables
if t == 1
@ -51,8 +51,10 @@ function _add_status_eqs!(
)
end
# Cannot switch on and off at the same time
eq_switch_on_off[g.name, t] =
@constraint(model, switch_on[g.name, t] + switch_off[g.name, t] <= 1)
eq_switch_on_off[g.name, t] = @constraint(
model,
switch_on[g.name, t] + switch_off[g.name, t] <= 1
)
end
end
return

@ -26,12 +26,12 @@ function _add_production_piecewise_linear_eqs!(
switch_on = model[:switch_on]
switch_off = model[:switch_off]
for t = 1:T
for k = 1:K
for t in 1:T
for k in 1:K
# Pbar^{k-1)
Pbar0 =
g.min_power[t] +
(k > 1 ? sum(g.cost_segments[ell].mw[t] for ell = 1:k-1) : 0.0)
(k > 1 ? sum(g.cost_segments[ell].mw[t] for ell in 1:k-1) : 0.0)
# Pbar^k
Pbar1 = g.cost_segments[k].mw[t] + Pbar0
@ -61,7 +61,8 @@ function _add_production_piecewise_linear_eqs!(
eq_segprod_limit_a[gn, t, k] = @constraint(
model,
segprod[gn, t, k] <=
g.cost_segments[k].mw[t] * is_on[gn, t] - Cv * switch_on[gn, t] -
g.cost_segments[k].mw[t] * is_on[gn, t] -
Cv * switch_on[gn, t] -
(t < T ? Cw * switch_off[gn, t+1] : 0.0)
)
else
@ -69,7 +70,8 @@ function _add_production_piecewise_linear_eqs!(
eq_segprod_limit_b[gn, t, k] = @constraint(
model,
segprod[gn, t, k] <=
g.cost_segments[k].mw[t] * is_on[gn, t] - Cv * switch_on[gn, t] -
g.cost_segments[k].mw[t] * is_on[gn, t] -
Cv * switch_on[gn, t] -
(t < T ? max(0, Cv - Cw) * switch_off[gn, t+1] : 0.0)
)
@ -85,12 +87,18 @@ function _add_production_piecewise_linear_eqs!(
# Definition of production
# Equation (43) in Kneuven et al. (2020)
eq_prod_above_def[gn, t] =
@constraint(model, prod_above[gn, t] == sum(segprod[gn, t, k] for k = 1:K))
eq_prod_above_def[gn, t] = @constraint(
model,
prod_above[gn, t] == sum(segprod[gn, t, k] for k in 1:K)
)
# Objective function
# Equation (44) in Kneuven et al. (2020)
add_to_expression!(model[:obj], segprod[gn, t, k], g.cost_segments[k].cost[t])
add_to_expression!(
model[:obj],
segprod[gn, t, k],
g.cost_segments[k].cost[t],
)
# Also add an explicit upper bound on segprod to make the solver's
# work a bit easier

@ -32,8 +32,9 @@ function _add_ramp_eqs!(
switch_off = model[:switch_off]
switch_on = model[:switch_on]
for t = 1:model[:instance].time
time_invariant = (t > 1) ? (abs(g.min_power[t] - g.min_power[t-1]) < 1e-7) : true
for t in 1:model[:instance].time
time_invariant =
(t > 1) ? (abs(g.min_power[t] - g.min_power[t-1]) < 1e-7) : true
# Ramp up limit
if t == 1
@ -59,8 +60,8 @@ function _add_ramp_eqs!(
g.min_power[t] * is_on[gn, t] +
prod_above[gn, t] +
(
RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? reserve[gn, t] :
0.0
RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ?
reserve[gn, t] : 0.0
)
min_prod_last_period =
g.min_power[t-1] * is_on[gn, t-1] + prod_above[gn, t-1]
@ -75,7 +76,8 @@ function _add_ramp_eqs!(
# prod_above[gn, t] when starting up, and creates diff with (24).
eq_ramp_up[gn, t] = @constraint(
model,
prod_above[gn, t] + (RESERVES_WHEN_RAMP_UP ? reserve[gn, t] : 0.0) -
prod_above[gn, t] +
(RESERVES_WHEN_RAMP_UP ? reserve[gn, t] : 0.0) -
prod_above[gn, t-1] <= RU
)
end
@ -105,7 +107,8 @@ function _add_ramp_eqs!(
RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ?
reserve[gn, t-1] : 0.0
)
min_prod_this_period = g.min_power[t] * is_on[gn, t] + prod_above[gn, t]
min_prod_this_period =
g.min_power[t] * is_on[gn, t] + prod_above[gn, t]
eq_ramp_down[gn, t] = @constraint(
model,
max_prod_last_period - min_prod_this_period <=

@ -11,27 +11,30 @@ function _add_startup_cost_eqs!(
eq_startup_restrict = _init(model, :eq_startup_restrict)
S = length(g.startup_categories)
startup = model[:startup]
for t = 1:model[:instance].time
for t in 1:model[:instance].time
# If unit is switching on, we must choose a startup category
eq_startup_choose[g.name, t] = @constraint(
model,
model[:switch_on][g.name, t] == sum(startup[g.name, t, s] for s = 1:S)
model[:switch_on][g.name, t] ==
sum(startup[g.name, t, s] for s in 1:S)
)
for s = 1:S
for s in 1:S
# If unit has not switched off in the last `delay` time periods, startup category is forbidden.
# The last startup category is always allowed.
if s < S
range_start = t - g.startup_categories[s+1].delay + 1
range_end = t - g.startup_categories[s].delay
range = (range_start:range_end)
initial_sum =
(g.initial_status < 0 && (g.initial_status + 1 in range) ? 1.0 : 0.0)
initial_sum = (
g.initial_status < 0 && (g.initial_status + 1 in range) ? 1.0 : 0.0
)
eq_startup_restrict[g.name, t, s] = @constraint(
model,
startup[g.name, t, s] <=
initial_sum +
sum(model[:switch_off][g.name, i] for i in range if i >= 1)
initial_sum + sum(
model[:switch_off][g.name, i] for i in range if i >= 1
)
)
end

@ -14,8 +14,10 @@ function _add_ramp_eqs!(
gn = g.name
reserve = model[:reserve]
eq_str_prod_limit = _init(model, :eq_str_prod_limit)
eq_prod_limit_ramp_up_extra_period = _init(model, :eq_prod_limit_ramp_up_extra_period)
eq_prod_limit_shutdown_trajectory = _init(model, :eq_prod_limit_shutdown_trajectory)
eq_prod_limit_ramp_up_extra_period =
_init(model, :eq_prod_limit_ramp_up_extra_period)
eq_prod_limit_shutdown_trajectory =
_init(model, :eq_prod_limit_shutdown_trajectory)
UT = g.min_uptime
SU = g.startup_limit # startup rate, i.e., max production right after startup
SD = g.shutdown_limit # shutdown rate, i.e., max production right before shutdown
@ -31,7 +33,7 @@ function _add_ramp_eqs!(
switch_off = model[:switch_off]
switch_on = model[:switch_on]
for t = 1:T
for t in 1:T
Pbar = g.max_power[t]
if Pbar < 1e-7
# Skip this time period if max power = 0
@ -52,10 +54,13 @@ function _add_ramp_eqs!(
# then switch_off[gn, t+1] = 0
eq_str_prod_limit[gn, t] = @constraint(
model,
prod_above[gn, t] + g.min_power[t] * is_on[gn, t] + reserve[gn, t] <=
Pbar * is_on[gn, t] - (t < T ? (Pbar - SD) * switch_off[gn, t+1] : 0.0) - sum(
prod_above[gn, t] +
g.min_power[t] * is_on[gn, t] +
reserve[gn, t] <=
Pbar * is_on[gn, t] -
(t < T ? (Pbar - SD) * switch_off[gn, t+1] : 0.0) - sum(
(Pbar - (SU + i * RU)) * switch_on[gn, t-i] for
i = 0:min(UT - 2, TRU, t - 1)
i in 0:min(UT - 2, TRU, t - 1)
)
)
@ -64,10 +69,12 @@ function _add_ramp_eqs!(
# Covers an additional time period of the ramp-up trajectory, compared to (38)
eq_prod_limit_ramp_up_extra_period[gn, t] = @constraint(
model,
prod_above[gn, t] + g.min_power[t] * is_on[gn, t] + reserve[gn, t] <=
prod_above[gn, t] +
g.min_power[t] * is_on[gn, t] +
reserve[gn, t] <=
Pbar * is_on[gn, t] - sum(
(Pbar - (SU + i * RU)) * switch_on[gn, t-i] for
i = 0:min(UT - 1, TRU, t - 1)
i in 0:min(UT - 1, TRU, t - 1)
)
)
end
@ -82,9 +89,13 @@ function _add_ramp_eqs!(
prod_above[gn, t] +
g.min_power[t] * is_on[gn, t] +
(RESERVES_WHEN_SHUT_DOWN ? reserve[gn, t] : 0.0) <=
Pbar * is_on[gn, t] -
sum((Pbar - (SD + i * RD)) * switch_off[gn, t+1+i] for i = 0:KSD) -
sum((Pbar - (SU + i * RU)) * switch_on[gn, t-i] for i = 0:KSU) - (
Pbar * is_on[gn, t] - sum(
(Pbar - (SD + i * RD)) * switch_off[gn, t+1+i] for
i in 0:KSD
) - sum(
(Pbar - (SU + i * RU)) * switch_on[gn, t-i] for
i in 0:KSU
) - (
(KSU >= TRU || KSU > t - 2) ? 0.0 :
max(0, (SU + (KSU + 1) * RU) - (SD + TRD * RD)) *
switch_on[gn, t-(KSU+1)]

@ -8,7 +8,7 @@ function _add_flexiramp_vars!(model::JuMP.Model, g::Unit)::Nothing
mfg = _init(model, :mfg)
dwflexiramp = _init(model, :dwflexiramp)
dwflexiramp_shortfall = _init(model, :dwflexiramp_shortfall)
for t = 1:model[:instance].time
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)
if g.provides_flexiramp_reserves[t]
@ -51,28 +51,37 @@ function _add_ramp_eqs!(
dwflexiramp = model[:dwflexiramp]
mfg = model[:mfg]
for t = 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)
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])
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]) <=
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])
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]) <=
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
@ -104,7 +113,8 @@ function _add_ramp_eqs!(
) # 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])) <=
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)
@ -113,7 +123,8 @@ function _add_ramp_eqs!(
@constraint(
model,
mfg[gn, t] <=
(SD * (is_on[gn, t] - is_on[gn, t+1])) + (maxp[t] * is_on[gn, t+1])
(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,
@ -141,13 +152,15 @@ function _add_ramp_eqs!(
) # 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]
-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])
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(

@ -5,12 +5,13 @@
function _add_bus!(model::JuMP.Model, b::Bus)::Nothing
net_injection = _init(model, :expr_net_injection)
curtail = _init(model, :curtail)
for t = 1:model[:instance].time
for t in 1:model[:instance].time
# Fixed load
net_injection[b.name, t] = AffExpr(-b.load[t])
# Load curtailment
curtail[b.name, t] = @variable(model, lower_bound = 0, upper_bound = b.load[t])
curtail[b.name, t] =
@variable(model, lower_bound = 0, upper_bound = b.load[t])
add_to_expression!(net_injection[b.name, t], curtail[b.name, t], 1.0)
add_to_expression!(

@ -8,9 +8,13 @@ function _add_transmission_line!(
f::ShiftFactorsFormulation,
)::Nothing
overflow = _init(model, :overflow)
for t = 1:model[:instance].time
for t in 1:model[:instance].time
overflow[lm.name, t] = @variable(model, lower_bound = 0)
add_to_expression!(model[:obj], overflow[lm.name, t], lm.flow_limit_penalty[t])
add_to_expression!(
model[:obj],
overflow[lm.name, t],
lm.flow_limit_penalty[t],
)
end
return
end

@ -2,18 +2,26 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
function _add_price_sensitive_load!(model::JuMP.Model, ps::PriceSensitiveLoad)::Nothing
function _add_price_sensitive_load!(
model::JuMP.Model,
ps::PriceSensitiveLoad,
)::Nothing
loads = _init(model, :loads)
net_injection = _init(model, :expr_net_injection)
for t = 1:model[:instance].time
for t in 1:model[:instance].time
# Decision variable
loads[ps.name, t] = @variable(model, lower_bound = 0, upper_bound = ps.demand[t])
loads[ps.name, t] =
@variable(model, lower_bound = 0, upper_bound = ps.demand[t])
# Objective function terms
add_to_expression!(model[:obj], loads[ps.name, t], -ps.revenue[t])
# Net injection
add_to_expression!(net_injection[ps.bus.name, t], loads[ps.name, t], -1.0)
add_to_expression!(
net_injection[ps.bus.name, t],
loads[ps.name, t],
-1.0,
)
end
return
end

@ -13,7 +13,10 @@ M[l.offset, b.offset] indicates the amount of power (in MW) that flows through
transmission line l when 1 MW of power is injected at the slack bus (the bus
that has offset zero) and withdrawn from b.
"""
function _injection_shift_factors(; buses::Array{Bus}, lines::Array{TransmissionLine})
function _injection_shift_factors(;
buses::Array{Bus},
lines::Array{TransmissionLine},
)
susceptance = _susceptance_matrix(lines)
incidence = _reduced_incidence_matrix(lines = lines, buses = buses)
laplacian = transpose(incidence) * susceptance * incidence
@ -30,7 +33,10 @@ is the number of buses and L is the number of lines. For each row, there is a 1
element and a -1 element, indicating the source and target buses, respectively,
for that line.
"""
function _reduced_incidence_matrix(; buses::Array{Bus}, lines::Array{TransmissionLine})
function _reduced_incidence_matrix(;
buses::Array{Bus},
lines::Array{TransmissionLine},
)
matrix = spzeros(Float64, length(lines), length(buses) - 1)
for line in lines
if line.source.offset > 0
@ -69,7 +75,7 @@ function _line_outage_factors(;
incidence = Array(_reduced_incidence_matrix(lines = lines, buses = buses))
lodf::Array{Float64,2} = isf * transpose(incidence)
_, n = size(lodf)
for i = 1:n
for i in 1:n
lodf[:, i] *= 1.0 / (1.0 - lodf[i, i])
lodf[i, i] = -1
end

@ -25,7 +25,14 @@ struct Formulation
status_vars::StatusVarsFormulation = Gar1962.StatusVars(),
transmission::TransmissionFormulation = ShiftFactorsFormulation(),
)
return new(prod_vars, pwl_costs, ramping, startup_costs, status_vars, transmission)
return new(
prod_vars,
pwl_costs,
ramping,
startup_costs,
status_vars,
transmission,
)
end
end

@ -14,12 +14,12 @@ function _add_net_injection_eqs!(model::JuMP.Model)::Nothing
net_injection = _init(model, :net_injection)
eq_net_injection = _init(model, :eq_net_injection)
eq_power_balance = _init(model, :eq_power_balance)
for t = 1:T, b in model[:instance].buses
for t in 1:T, b in model[:instance].buses
n = net_injection[b.name, t] = @variable(model)
eq_net_injection[b.name, t] =
@constraint(model, -n + model[:expr_net_injection][b.name, t] == 0)
end
for t = 1:T
for t in 1:T
eq_power_balance[t] = @constraint(
model,
sum(net_injection[b.name, t] for b in model[:instance].buses) == 0
@ -31,7 +31,7 @@ end
function _add_reserve_eqs!(model::JuMP.Model)::Nothing
eq_min_reserve = _init(model, :eq_min_reserve)
instance = model[:instance]
for t = 1:instance.time
for t in 1:instance.time
# Equation (68) in Kneuven et al. (2020)
# As in Morales-España et al. (2013a)
# Akin to the alternative formulation with max_power_avail
@ -46,7 +46,11 @@ function _add_reserve_eqs!(model::JuMP.Model)::Nothing
# Account for shortfall contribution to objective
if shortfall_penalty >= 0
add_to_expression!(model[:obj], shortfall_penalty, model[:reserve_shortfall][t])
add_to_expression!(
model[:obj],
shortfall_penalty,
model[:reserve_shortfall][t],
)
end
end
return
@ -61,20 +65,24 @@ function _add_flexiramp_eqs!(model::JuMP.Model)::Nothing
eq_min_upflexiramp = _init(model, :eq_min_upflexiramp)
eq_min_dwflexiramp = _init(model, :eq_min_dwflexiramp)
instance = model[:instance]
for t = 1:instance.time
for t in 1:instance.time
flexiramp_shortfall_penalty = instance.flexiramp_shortfall_penalty[t]
# Eq. (17) in Wang & Hobbs (2016)
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
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]
)
# 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
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]
)
@ -83,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

@ -46,7 +46,7 @@ _is_initially_on(g::Unit)::Float64 = (g.initial_status > 0 ? 1.0 : 0.0)
function _add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
reserve = _init(model, :reserve)
reserve_shortfall = _init(model, :reserve_shortfall)
for t = 1:model[:instance].time
for t in 1:model[:instance].time
if g.provides_spinning_reserves[t]
reserve[g.name, t] = @variable(model, lower_bound = 0)
else
@ -61,7 +61,7 @@ end
function _add_reserve_eqs!(model::JuMP.Model, g::Unit)::Nothing
reserve = model[:reserve]
for t = 1:model[:instance].time
for t in 1:model[:instance].time
add_to_expression!(expr_reserve[g.bus.name, t], reserve[g.name, t], 1.0)
end
return
@ -69,8 +69,8 @@ end
function _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing
startup = _init(model, :startup)
for t = 1:model[:instance].time
for s = 1:length(g.startup_categories)
for t in 1:model[:instance].time
for s in 1:length(g.startup_categories)
startup[g.name, t, s] = @variable(model, binary = true)
end
end
@ -86,7 +86,7 @@ function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
switch_off = model[:switch_off]
switch_on = model[:switch_on]
T = model[:instance].time
for t = 1:T
for t in 1:T
# Startup limit
eq_startup_limit[g.name, t] = @constraint(
model,
@ -96,14 +96,16 @@ function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
)
# Shutdown limit
if g.initial_power > g.shutdown_limit
eq_shutdown_limit[g.name, 0] = @constraint(model, switch_off[g.name, 1] <= 0)
eq_shutdown_limit[g.name, 0] =
@constraint(model, switch_off[g.name, 1] <= 0)
end
if t < T
eq_shutdown_limit[g.name, t] = @constraint(
model,
prod_above[g.name, t] <=
(g.max_power[t] - g.min_power[t]) * is_on[g.name, t] -
max(0, g.max_power[t] - g.shutdown_limit) * switch_off[g.name, t+1]
max(0, g.max_power[t] - g.shutdown_limit) *
switch_off[g.name, t+1]
)
end
end
@ -119,7 +121,7 @@ function _add_ramp_eqs!(
reserve = model[:reserve]
eq_ramp_up = _init(model, :eq_ramp_up)
eq_ramp_down = _init(model, :eq_ramp_down)
for t = 1:model[:instance].time
for t in 1:model[:instance].time
# Ramp up limit
if t == 1
if _is_initially_on(g) == 1
@ -149,7 +151,8 @@ function _add_ramp_eqs!(
else
eq_ramp_down[g.name, t] = @constraint(
model,
prod_above[g.name, t] >= prod_above[g.name, t-1] - g.ramp_down_limit
prod_above[g.name, t] >=
prod_above[g.name, t-1] - g.ramp_down_limit
)
end
end
@ -162,18 +165,18 @@ function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
eq_min_uptime = _init(model, :eq_min_uptime)
eq_min_downtime = _init(model, :eq_min_downtime)
T = model[:instance].time
for t = 1:T
for t in 1:T
# Minimum up-time
eq_min_uptime[g.name, t] = @constraint(
model,
sum(switch_on[g.name, i] for i in (t-g.min_uptime+1):t if i >= 1) <=
is_on[g.name, t]
sum(switch_on[g.name, i] for i in (t-g.min_uptime+1):t if i >= 1) <= is_on[g.name, t]
)
# Minimum down-time
eq_min_downtime[g.name, t] = @constraint(
model,
sum(switch_off[g.name, i] for i in (t-g.min_downtime+1):t if i >= 1) <=
1 - is_on[g.name, t]
sum(
switch_off[g.name, i] for i in (t-g.min_downtime+1):t if i >= 1
) <= 1 - is_on[g.name, t]
)
# Minimum up/down-time for initial periods
if t == 1
@ -200,7 +203,7 @@ end
function _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing
expr_net_injection = model[:expr_net_injection]
for t = 1:model[:instance].time
for t in 1:model[:instance].time
# Add to net injection expression
add_to_expression!(
expr_net_injection[g.bus.name, t],

@ -14,10 +14,12 @@ function fix!(model::JuMP.Model, solution::AbstractDict)::Nothing
prod_above = model[:prod_above]
reserve = model[:reserve]
for g in instance.units
for t = 1:T
for t in 1:T
is_on_value = round(solution["Is on"][g.name][t])
prod_value = round(solution["Production (MW)"][g.name][t], digits = 5)
reserve_value = round(solution["Reserve (MW)"][g.name][t], digits = 5)
prod_value =
round(solution["Production (MW)"][g.name][t], digits = 5)
reserve_value =
round(solution["Reserve (MW)"][g.name][t], digits = 5)
JuMP.fix(is_on[g.name, t], is_on_value, force = true)
JuMP.fix(
prod_above[g.name, t],

@ -2,7 +2,10 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
function _enforce_transmission(model::JuMP.Model, violations::Vector{_Violation})::Nothing
function _enforce_transmission(
model::JuMP.Model,
violations::Vector{_Violation},
)::Nothing
for v in violations
_enforce_transmission(
model = model,

@ -4,9 +4,11 @@
function _offer(filter::_ViolationFilter, v::_Violation)::Nothing
if v.monitored_line.offset keys(filter.queues)
filter.queues[v.monitored_line.offset] = PriorityQueue{_Violation,Float64}()
filter.queues[v.monitored_line.offset] =
PriorityQueue{_Violation,Float64}()
end
q::PriorityQueue{_Violation,Float64} = filter.queues[v.monitored_line.offset]
q::PriorityQueue{_Violation,Float64} =
filter.queues[v.monitored_line.offset]
if length(q) < filter.max_per_line
enqueue!(q, v => v.amount)
else

@ -4,7 +4,11 @@
import Base.Threads: @threads
function _find_violations(model::JuMP.Model; max_per_line::Int, max_per_period::Int)
function _find_violations(
model::JuMP.Model;
max_per_line::Int,
max_per_period::Int,
)
instance = model[:instance]
net_injection = model[:net_injection]
overflow = model[:overflow]
@ -14,10 +18,13 @@ function _find_violations(model::JuMP.Model; max_per_line::Int, max_per_period::
time_screening = @elapsed begin
non_slack_buses = [b for b in instance.buses if b.offset > 0]
net_injection_values = [
value(net_injection[b.name, t]) for b in non_slack_buses, t = 1:instance.time
value(net_injection[b.name, t]) for b in non_slack_buses,
t in 1:instance.time
]
overflow_values = [
value(overflow[lm.name, t]) for lm in instance.lines,
t in 1:instance.time
]
overflow_values =
[value(overflow[lm.name, t]) for lm in instance.lines, t = 1:instance.time]
violations = UnitCommitment._find_violations(
instance = instance,
net_injections = net_injection_values,
@ -28,7 +35,10 @@ function _find_violations(model::JuMP.Model; max_per_line::Int, max_per_period::
max_per_period = max_per_period,
)
end
@info @sprintf("Verified transmission limits in %.2f seconds", time_screening)
@info @sprintf(
"Verified transmission limits in %.2f seconds",
time_screening
)
return violations
end
@ -71,8 +81,10 @@ function _find_violations(;
size(lodf) == (L, L) || error("lodf has incorrect size")
filters = Dict(
t => _ViolationFilter(max_total = max_per_period, max_per_line = max_per_line)
for t = 1:T
t => _ViolationFilter(
max_total = max_per_period,
max_per_line = max_per_line,
) for t in 1:T
)
pre_flow::Array{Float64} = zeros(L, K) # pre_flow[lm, thread]
@ -80,30 +92,35 @@ function _find_violations(;
pre_v::Array{Float64} = zeros(L, K) # pre_v[lm, thread]
post_v::Array{Float64} = zeros(L, L, K) # post_v[lm, lc, thread]
normal_limits::Array{Float64,2} =
[l.normal_flow_limit[t] + overflow[l.offset, t] for l in instance.lines, t = 1:T]
normal_limits::Array{Float64,2} = [
l.normal_flow_limit[t] + overflow[l.offset, t] for
l in instance.lines, t in 1:T
]
emergency_limits::Array{Float64,2} =
[l.emergency_flow_limit[t] + overflow[l.offset, t] for l in instance.lines, t = 1:T]
emergency_limits::Array{Float64,2} = [
l.emergency_flow_limit[t] + overflow[l.offset, t] for
l in instance.lines, t in 1:T
]
is_vulnerable::Array{Bool} = zeros(Bool, L)
for c in instance.contingencies
is_vulnerable[c.lines[1].offset] = true
end
@threads for t = 1:T
@threads for t in 1:T
k = threadid()
# Pre-contingency flows
pre_flow[:, k] = isf * net_injections[:, t]
# Post-contingency flows
for lc = 1:L, lm = 1:L
post_flow[lm, lc, k] = pre_flow[lm, k] + pre_flow[lc, k] * lodf[lm, lc]
for lc in 1:L, lm in 1:L
post_flow[lm, lc, k] =
pre_flow[lm, k] + pre_flow[lc, k] * lodf[lm, lc]
end
# Pre-contingency violations
for lm = 1:L
for lm in 1:L
pre_v[lm, k] = max(
0.0,
pre_flow[lm, k] - normal_limits[lm, t],
@ -112,7 +129,7 @@ function _find_violations(;
end
# Post-contingency violations
for lc = 1:L, lm = 1:L
for lc in 1:L, lm in 1:L
post_v[lm, lc, k] = max(
0.0,
post_flow[lm, lc, k] - emergency_limits[lm, t],
@ -121,7 +138,7 @@ function _find_violations(;
end
# Offer pre-contingency violations
for lm = 1:L
for lm in 1:L
if pre_v[lm, k] > 1e-5
_offer(
filters[t],
@ -136,7 +153,7 @@ function _find_violations(;
end
# Offer post-contingency violations
for lm = 1:L, lc = 1:L
for lm in 1:L, lc in 1:L
if post_v[lm, lc, k] > 1e-5 && is_vulnerable[lc]
_offer(
filters[t],
@ -152,7 +169,7 @@ function _find_violations(;
end
violations = _Violation[]
for t = 1:instance.time
for t in 1:instance.time
append!(violations, _query(filters[t]))
end

@ -27,7 +27,10 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing
@info "Time limit exceeded"
break
end
@info @sprintf("Setting MILP time limit to %.2f seconds", time_remaining)
@info @sprintf(
"Setting MILP time limit to %.2f seconds",
time_remaining
)
JuMP.set_time_limit_sec(model, time_remaining)
@info "Solving MILP..."
JuMP.optimize!(model)

@ -6,40 +6,43 @@ function solution(model::JuMP.Model)::OrderedDict
instance, T = model[:instance], model[:instance].time
function timeseries(vars, collection)
return OrderedDict(
b.name => [round(value(vars[b.name, t]), digits = 5) for t = 1:T] for
b in collection
b.name => [round(value(vars[b.name, t]), digits = 5) for t in 1:T]
for b in collection
)
end
function production_cost(g)
return [
value(model[:is_on][g.name, t]) * g.min_power_cost[t] + sum(
Float64[
value(model[:segprod][g.name, t, k]) * g.cost_segments[k].cost[t] for
k = 1:length(g.cost_segments)
value(model[:segprod][g.name, t, k]) *
g.cost_segments[k].cost[t] for
k in 1:length(g.cost_segments)
],
) for t = 1:T
) for t in 1:T
]
end
function production(g)
return [
value(model[:is_on][g.name, t]) * g.min_power[t] + sum(
Float64[
value(model[:segprod][g.name, t, k]) for k = 1:length(g.cost_segments)
value(model[:segprod][g.name, t, k]) for
k in 1:length(g.cost_segments)
],
) for t = 1:T
) for t in 1:T
]
end
function startup_cost(g)
S = length(g.startup_categories)
return [
sum(
g.startup_categories[s].cost * value(model[:startup][g.name, t, s]) for
s = 1:S
) for t = 1:T
g.startup_categories[s].cost *
value(model[:startup][g.name, t, s]) for s in 1:S
) for t in 1:T
]
end
sol = OrderedDict()
sol["Production (MW)"] = OrderedDict(g.name => production(g) for g in instance.units)
sol["Production (MW)"] =
OrderedDict(g.name => production(g) for g in instance.units)
sol["Production cost (\$)"] =
OrderedDict(g.name => production_cost(g) for g in instance.units)
sol["Startup cost (\$)"] =
@ -51,19 +54,21 @@ function solution(model::JuMP.Model)::OrderedDict
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 = 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 = 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
@ -72,11 +77,12 @@ 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 = 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)
sol["Net injection (MW)"] =
timeseries(model[:net_injection], instance.buses)
sol["Load curtail (MW)"] = timeseries(model[:curtail], instance.buses)
if !isempty(instance.lines)
sol["Line overflow (MW)"] = timeseries(model[:overflow], instance.lines)

@ -6,10 +6,16 @@ function set_warm_start!(model::JuMP.Model, solution::AbstractDict)::Nothing
instance, T = model[:instance], model[:instance].time
is_on = model[:is_on]
for g in instance.units
for t = 1:T
for t in 1:T
JuMP.set_start_value(is_on[g.name, t], solution["Is on"][g.name][t])
JuMP.set_start_value(switch_on[g.name, t], solution["Switch on"][g.name][t])
JuMP.set_start_value(switch_off[g.name, t], solution["Switch off"][g.name][t])
JuMP.set_start_value(
switch_on[g.name, t],
solution["Switch on"][g.name][t],
)
JuMP.set_start_value(
switch_off[g.name, t],
solution["Switch off"][g.name][t],
)
end
end
return

@ -11,7 +11,10 @@ Generates feasible initial conditions for the given instance, by constructing
and solving a single-period mixed-integer optimization problem, using the given
optimizer. The instance is modified in-place.
"""
function generate_initial_conditions!(instance::UnitCommitmentInstance, optimizer)::Nothing
function generate_initial_conditions!(
instance::UnitCommitmentInstance,
optimizer,
)::Nothing
G = instance.units
B = instance.buses
t = 1
@ -28,7 +31,11 @@ function generate_initial_conditions!(instance::UnitCommitmentInstance, optimize
@constraint(mip, max_power[g in G], p[g] <= g.max_power[t] * x[g])
# Constraint: Production equals demand
@constraint(mip, power_balance, sum(b.load[t] for b in B) == sum(p[g] for g in G))
@constraint(
mip,
power_balance,
sum(b.load[t] for b in B) == sum(p[g] for g in G)
)
# Constraint: Must run
for g in G

@ -117,7 +117,10 @@ Base.@kwdef struct Randomization
randomize_load_share::Bool = true
end
function _randomize_costs(instance::UnitCommitmentInstance, distribution)::Nothing
function _randomize_costs(
instance::UnitCommitmentInstance,
distribution,
)::Nothing
for unit in instance.units
α = rand(distribution)
unit.min_power_cost *= α
@ -131,11 +134,17 @@ function _randomize_costs(instance::UnitCommitmentInstance, distribution)::Nothi
return
end
function _randomize_load_share(instance::UnitCommitmentInstance, distribution)::Nothing
function _randomize_load_share(
instance::UnitCommitmentInstance,
distribution,
)::Nothing
α = rand(distribution, length(instance.buses))
for t = 1:instance.time
for t in 1:instance.time
total = sum(bus.load[t] for bus in instance.buses)
den = sum(bus.load[t] / total * α[i] for (i, bus) in enumerate(instance.buses))
den = sum(
bus.load[t] / total * α[i] for
(i, bus) in enumerate(instance.buses)
)
for (i, bus) in enumerate(instance.buses)
bus.load[t] *= α[i] / den
end
@ -149,9 +158,11 @@ function _randomize_load_profile(
)::Nothing
# Generate new system load
system_load = [1.0]
for t = 2:instance.time
for t in 2:instance.time
idx = (t - 1) % length(params.load_profile_mu) + 1
gamma = rand(Normal(params.load_profile_mu[idx], params.load_profile_sigma[idx]))
gamma = rand(
Normal(params.load_profile_mu[idx], params.load_profile_sigma[idx]),
)
push!(system_load, system_load[t-1] * gamma)
end
capacity = sum(maximum(u.max_power) for u in instance.units)
@ -161,7 +172,7 @@ function _randomize_load_profile(
# Scale bus loads to match the new system load
prev_system_load = sum(b.load for b in instance.buses)
for b in instance.buses
for t = 1:instance.time
for t in 1:instance.time
b.load[t] *= system_load[t] / prev_system_load[t]
end
end

@ -93,8 +93,9 @@ function _run_benchmarks(;
trials,
)
combinations = [
(c, s.first, s.second, m.first, m.second, f.first, f.second, t) for c in cases
for s in optimizers for f in formulations for m in methods for t in trials
(c, s.first, s.second, m.first, m.second, f.first, f.second, t) for
c in cases for s in optimizers for f in formulations for
m in methods for t in trials
]
shuffle!(combinations)
if nworkers() > 1

@ -18,7 +18,7 @@ function repair!(instance::UnitCommitmentInstance)::Int
for g in instance.units
# Startup costs and delays must be increasing
for s = 2:length(g.startup_categories)
for s in 2:length(g.startup_categories)
if g.startup_categories[s].delay <= g.startup_categories[s-1].delay
prev_value = g.startup_categories[s].delay
new_value = g.startup_categories[s-1].delay + 1
@ -38,9 +38,9 @@ function repair!(instance::UnitCommitmentInstance)::Int
end
end
for t = 1:instance.time
for t in 1:instance.time
# Production cost curve should be convex
for k = 2:length(g.cost_segments)
for k in 2:length(g.cost_segments)
cost = g.cost_segments[k].cost[t]
min_cost = g.cost_segments[k-1].cost[t]
if cost < min_cost - 1e-5

@ -24,7 +24,10 @@ This function is implemented independently from the optimization model in
producing valid solutions. It can also be used to verify the solutions produced
by other optimization packages.
"""
function validate(instance::UnitCommitmentInstance, solution::Union{Dict,OrderedDict})::Bool
function validate(
instance::UnitCommitmentInstance,
solution::Union{Dict,OrderedDict},
)::Bool
err_count = 0
err_count += _validate_units(instance, solution)
err_count += _validate_reserve_and_demand(instance, solution)
@ -47,12 +50,13 @@ function _validate_units(instance, solution; tol = 0.01)
actual_startup_cost = solution["Startup cost (\$)"][unit.name]
is_on = bin(solution["Is on"][unit.name])
for t = 1:instance.time
for t in 1:instance.time
# Auxiliary variables
if t == 1
is_starting_up = (unit.initial_status < 0) && is_on[t]
is_shutting_down = (unit.initial_status > 0) && !is_on[t]
ramp_up = max(0, production[t] + reserve[t] - unit.initial_power)
ramp_up =
max(0, production[t] + reserve[t] - unit.initial_power)
ramp_down = max(0, unit.initial_power - production[t])
else
is_starting_up = !is_on[t-1] && is_on[t]
@ -86,7 +90,11 @@ function _validate_units(instance, solution; tol = 0.01)
# Verify must-run
if !is_on[t] && unit.must_run[t]
@error @sprintf("Must-run unit %s is offline at time %d", unit.name, t)
@error @sprintf(
"Must-run unit %s is offline at time %d",
unit.name,
t
)
err_count += 1
end
@ -113,7 +121,8 @@ function _validate_units(instance, solution; tol = 0.01)
end
# If unit is on, must produce at most its maximum power
if is_on[t] && (production[t] + reserve[t] > unit.max_power[t] + tol)
if is_on[t] &&
(production[t] + reserve[t] > unit.max_power[t] + tol)
@error @sprintf(
"Unit %s produces above its maximum limit at time %d (%.2f + %.2f> %.2f)",
unit.name,
@ -127,7 +136,11 @@ function _validate_units(instance, solution; tol = 0.01)
# If unit is off, must produce zero
if !is_on[t] && production[t] + reserve[t] > tol
@error @sprintf("Unit %s produces power at time %d while off", unit.name, t)
@error @sprintf(
"Unit %s produces power at time %d while off",
unit.name,
t
)
err_count += 1
end
@ -156,7 +169,9 @@ function _validate_units(instance, solution; tol = 0.01)
end
# Ramp-up limit
if !is_starting_up && !is_shutting_down && (ramp_up > unit.ramp_up_limit + tol)
if !is_starting_up &&
!is_shutting_down &&
(ramp_up > unit.ramp_up_limit + tol)
@error @sprintf(
"Unit %s exceeds ramp up limit at time %d (%.2f > %.2f)",
unit.name,
@ -186,7 +201,7 @@ function _validate_units(instance, solution; tol = 0.01)
# Calculate how much time the unit has been offline
time_down = 0
for k = 1:(t-1)
for k in 1:(t-1)
if !is_on[t-k]
time_down += 1
else
@ -220,7 +235,7 @@ function _validate_units(instance, solution; tol = 0.01)
# Calculate how much time the unit has been online
time_up = 0
for k = 1:(t-1)
for k in 1:(t-1)
if is_on[t-k]
time_up += 1
else
@ -273,7 +288,7 @@ end
function _validate_reserve_and_demand(instance, solution, tol = 0.01)
err_count = 0
for t = 1:instance.time
for t in 1:instance.time
load_curtail = 0
fixed_load = sum(b.load[t] for b in instance.buses)
ps_load = 0
@ -283,10 +298,13 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
ps in instance.price_sensitive_loads
)
end
production = sum(solution["Production (MW)"][g.name][t] for g in instance.units)
production =
sum(solution["Production (MW)"][g.name][t] for g in instance.units)
if "Load curtail (MW)" in keys(solution)
load_curtail =
sum(solution["Load curtail (MW)"][b.name][t] for b in instance.buses)
load_curtail = sum(
solution["Load curtail (MW)"][b.name][t] for
b in instance.buses
)
end
balance = fixed_load - load_curtail - production + ps_load
@ -303,18 +321,20 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
err_count += 1
end
# Verify flexiramp solutions only if either of the up-flexiramp and
# down-flexiramp requirements is not a default array of zeros
if instance.reserves.upflexiramp != zeros(instance.time) ||
instance.reserves.dwflexiramp != zeros(instance.time)
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,
@ -325,14 +345,16 @@ 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,
@ -345,7 +367,8 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
# Verify spinning reserve solutions only if both up-flexiramp and
# down-flexiramp requirements are arrays of zeros.
else
reserve = sum(solution["Reserve (MW)"][g.name][t] for g in instance.units)
reserve =
sum(solution["Reserve (MW)"][g.name][t] for g in instance.units)
reserve_shortfall =
(instance.shortfall_penalty[t] >= 0) ?
solution["Reserve shortfall (MW)"][t] : 0
@ -361,7 +384,6 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
err_count += 1
end
end
end
return err_count

@ -7,8 +7,9 @@ using UnitCommitment
basedir = @__DIR__
@testset "read_egret_solution" begin
solution =
UnitCommitment.read_egret_solution("$basedir/../fixtures/egret_output.json.gz")
solution = UnitCommitment.read_egret_solution(
"$basedir/../fixtures/egret_output.json.gz",
)
for attr in ["Is on", "Production (MW)", "Production cost (\$)"]
@test attr in keys(solution)
@test "115_STEAM_1" in keys(solution[attr])

@ -19,9 +19,9 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
@test instance.lines[5].target.name == "b5"
@test instance.lines[5].reactance 0.17388
@test instance.lines[5].susceptance 10.037550333
@test instance.lines[5].normal_flow_limit == [1e8 for t = 1:4]
@test instance.lines[5].emergency_flow_limit == [1e8 for t = 1:4]
@test instance.lines[5].flow_limit_penalty == [5e3 for t = 1:4]
@test instance.lines[5].normal_flow_limit == [1e8 for t in 1:4]
@test instance.lines[5].emergency_flow_limit == [1e8 for t in 1:4]
@test instance.lines[5].flow_limit_penalty == [5e3 for t in 1:4]
@test instance.lines_by_name["l5"].name == "l5"
@test instance.lines[1].name == "l1"
@ -29,9 +29,9 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
@test instance.lines[1].target.name == "b2"
@test instance.lines[1].reactance 0.059170
@test instance.lines[1].susceptance 29.496860773945
@test instance.lines[1].normal_flow_limit == [300.0 for t = 1:4]
@test instance.lines[1].emergency_flow_limit == [400.0 for t = 1:4]
@test instance.lines[1].flow_limit_penalty == [1e3 for t = 1:4]
@test instance.lines[1].normal_flow_limit == [300.0 for t in 1:4]
@test instance.lines[1].emergency_flow_limit == [400.0 for t in 1:4]
@test instance.lines[1].flow_limit_penalty == [1e3 for t in 1:4]
@test instance.buses[9].name == "b9"
@test instance.buses[9].load == [35.36638, 33.25495, 31.67138, 31.14353]
@ -44,12 +44,12 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
@test unit.ramp_down_limit == 1e6
@test unit.startup_limit == 1e6
@test unit.shutdown_limit == 1e6
@test unit.must_run == [false for t = 1:4]
@test unit.min_power_cost == [1400.0 for t = 1:4]
@test unit.must_run == [false for t in 1:4]
@test unit.min_power_cost == [1400.0 for t in 1:4]
@test unit.min_uptime == 1
@test unit.min_downtime == 1
@test unit.provides_spinning_reserves == [true for t = 1:4]
for t = 1:1
@test unit.provides_spinning_reserves == [true for t in 1:4]
for t in 1:1
@test unit.cost_segments[1].mw[t] == 10.0
@test unit.cost_segments[2].mw[t] == 20.0
@test unit.cost_segments[3].mw[t] == 5.0
@ -68,7 +68,7 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
unit = instance.units[2]
@test unit.name == "g2"
@test unit.must_run == [false for t = 1:4]
@test unit.must_run == [false for t in 1:4]
unit = instance.units[3]
@test unit.name == "g3"
@ -77,12 +77,12 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
@test unit.ramp_down_limit == 70.0
@test unit.startup_limit == 70.0
@test unit.shutdown_limit == 70.0
@test unit.must_run == [true for t = 1:4]
@test unit.min_power_cost == [0.0 for t = 1:4]
@test unit.must_run == [true for t in 1:4]
@test unit.min_power_cost == [0.0 for t in 1:4]
@test unit.min_uptime == 1
@test unit.min_downtime == 1
@test unit.provides_spinning_reserves == [true for t = 1:4]
for t = 1:4
@test unit.provides_spinning_reserves == [true for t in 1:4]
for t in 1:4
@test unit.cost_segments[1].mw[t] 33
@test unit.cost_segments[2].mw[t] 33
@test unit.cost_segments[3].mw[t] 34
@ -101,8 +101,8 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
load = instance.price_sensitive_loads[1]
@test load.name == "ps1"
@test load.bus.name == "b3"
@test load.revenue == [100.0 for t = 1:4]
@test load.demand == [50.0 for t = 1:4]
@test load.revenue == [100.0 for t in 1:4]
@test load.demand == [50.0 for t in 1:4]
@test instance.price_sensitive_loads_by_name["ps1"].name == "ps1"
end

@ -7,7 +7,7 @@ import UnitCommitment: _Violation, _offer, _query
@testset "find_violations" begin
instance = UnitCommitment.read_benchmark("test/case14")
for line in instance.lines, t = 1:instance.time
for line in instance.lines, t in 1:instance.time
line.normal_flow_limit[t] = 1.0
line.emergency_flow_limit[t] = 1.0
end
@ -20,8 +20,8 @@ import UnitCommitment: _Violation, _offer, _query
buses = instance.buses,
isf = isf,
)
inj = [1000.0 for b = 1:13, t = 1:instance.time]
overflow = [0.0 for l in instance.lines, t = 1:instance.time]
inj = [1000.0 for b in 1:13, t in 1:instance.time]
overflow = [0.0 for l in instance.lines, t in 1:instance.time]
violations = UnitCommitment._find_violations(
instance = instance,
net_injections = inj,

@ -8,7 +8,8 @@ basedir = @__DIR__
@testset "generate_initial_conditions!" begin
# Load instance
instance = UnitCommitment.read("$basedir/../fixtures/case118-initcond.json.gz")
instance =
UnitCommitment.read("$basedir/../fixtures/case118-initcond.json.gz")
optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
# All units should have unknown initial conditions

@ -28,7 +28,10 @@ test_approx(x, y) = @test isapprox(x, y, atol = 1e-3)
test_approx(bus.load[1] / prev_system_load[1], 0.012)
Random.seed!(42)
randomize!(instance, XavQiuAhm2021.Randomization(randomize_load_profile = false))
randomize!(
instance,
XavQiuAhm2021.Randomization(randomize_load_profile = false),
)
# Check randomized costs
test_approx(unit.min_power_cost[1], 831.977)

@ -6,7 +6,7 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON
@testset "build_model" begin
instance = UnitCommitment.read_benchmark("test/case14")
for line in instance.lines, t = 1:4
for line in instance.lines, t in 1:4
line.normal_flow_limit[t] = 10.0
end
optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)

Loading…
Cancel
Save