mirror of
https://github.com/ANL-CEEESA/UnitCommitment.jl.git
synced 2025-12-06 08:18:51 -06:00
new formatting
This commit is contained in:
BIN
benchmark/.DS_Store
vendored
Normal file
BIN
benchmark/.DS_Store
vendored
Normal file
Binary file not shown.
@@ -129,19 +129,15 @@ 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 in 1:T
|
||||
for t = 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,10 +23,7 @@ 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"
|
||||
@@ -65,9 +62,7 @@ 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
|
||||
@@ -97,8 +92,7 @@ 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
|
||||
|
||||
@@ -108,23 +102,23 @@ function _from_json(json; repair = true)
|
||||
|
||||
function timeseries(x; default = nothing)
|
||||
x !== nothing || return default
|
||||
x isa Array || return [x for t in 1:T]
|
||||
x isa Array || return [x for t = 1:T]
|
||||
return x
|
||||
end
|
||||
|
||||
# Read parameters
|
||||
power_balance_penalty = timeseries(
|
||||
json["Parameters"]["Power balance penalty (\$/MW)"],
|
||||
default = [1000.0 for t in 1:T],
|
||||
default = [1000.0 for t = 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 in 1:T],
|
||||
default = [500.0 for t = 1:T],
|
||||
)
|
||||
shortfall_penalty = timeseries(
|
||||
json["Parameters"]["Reserve shortfall penalty (\$/MW)"],
|
||||
default = [-1.0 for t in 1:T],
|
||||
default = [-1.0 for t = 1:T],
|
||||
)
|
||||
|
||||
# Read buses
|
||||
@@ -146,17 +140,14 @@ 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 in 1:K]...,
|
||||
)
|
||||
curve_cost = hcat(
|
||||
[timeseries(dict["Production cost curve (\$)"][k]) for k in 1:K]...,
|
||||
)
|
||||
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]...)
|
||||
min_power = curve_mw[:, 1]
|
||||
max_power = curve_mw[:, K]
|
||||
min_power_cost = curve_cost[:, 1]
|
||||
segments = CostSegment[]
|
||||
for k in 2:K
|
||||
for k = 2:K
|
||||
amount = curve_mw[:, k] - curve_mw[:, k-1]
|
||||
cost = (curve_cost[:, k] - curve_cost[:, k-1]) ./ amount
|
||||
replace!(cost, NaN => 0.0)
|
||||
@@ -167,13 +158,10 @@ 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 in 1:length(startup_delays)
|
||||
for k = 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
|
||||
|
||||
@@ -186,8 +174,7 @@ 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
|
||||
@@ -199,7 +186,7 @@ function _from_json(json; repair = true)
|
||||
bus,
|
||||
max_power,
|
||||
min_power,
|
||||
timeseries(dict["Must run?"], default = [false for t in 1:T]),
|
||||
timeseries(dict["Must run?"], default = [false for t = 1:T]),
|
||||
min_power_cost,
|
||||
segments,
|
||||
scalar(dict["Minimum uptime (h)"], default = 1) * time_multiplier,
|
||||
@@ -210,14 +197,8 @@ 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 in 1:T],
|
||||
),
|
||||
timeseries(
|
||||
dict["Provides flexible capacity?"],
|
||||
default = [true for t in 1:T],
|
||||
),
|
||||
timeseries(dict["Provides spinning reserves?"], default = [true for t = 1:T]),
|
||||
timeseries(dict["Provides flexible capacity?"], default = [true for t = 1:T]),
|
||||
startup_categories,
|
||||
)
|
||||
push!(bus.units, unit)
|
||||
@@ -230,14 +211,10 @@ 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
|
||||
@@ -250,17 +227,11 @@ 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 in 1:T],
|
||||
),
|
||||
timeseries(
|
||||
dict["Emergency flow limit (MW)"],
|
||||
default = [1e8 for t in 1:T],
|
||||
),
|
||||
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["Flow limit penalty (\$/MW)"],
|
||||
default = [5000.0 for t in 1:T],
|
||||
default = [5000.0 for t = 1:T],
|
||||
),
|
||||
)
|
||||
name_to_line[line_name] = line
|
||||
@@ -274,12 +245,10 @@ 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,10 +96,7 @@ 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
|
||||
|
||||
@@ -33,19 +33,15 @@ function build_model(;
|
||||
variable_names::Bool = false,
|
||||
)::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",
|
||||
)
|
||||
instance.reserves.spinning >= ones(instance.time) .* 1e-6
|
||||
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",
|
||||
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")
|
||||
end
|
||||
|
||||
@info "Building model..."
|
||||
|
||||
@@ -32,7 +32,7 @@ function _add_ramp_eqs!(
|
||||
switch_off = model[:switch_off]
|
||||
switch_on = model[:switch_on]
|
||||
|
||||
for t in 1:model[:instance].time
|
||||
for t = 1:model[:instance].time
|
||||
# Ramp up limit
|
||||
if t == 1
|
||||
if is_initially_on
|
||||
@@ -49,12 +49,8 @@ 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(
|
||||
@@ -81,11 +77,10 @@ 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,18 +18,16 @@ function _add_production_piecewise_linear_eqs!(
|
||||
prod_above = model[:prod_above]
|
||||
|
||||
K = length(g.cost_segments)
|
||||
for t in 1:model[:instance].time
|
||||
for t = 1:model[:instance].time
|
||||
gn = g.name
|
||||
for k in 1:K
|
||||
for k = 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
|
||||
@@ -37,18 +35,12 @@ 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 in 1:K)
|
||||
)
|
||||
eq_prod_above_def[gn, t] =
|
||||
@constraint(model, prod_above[gn, t] == sum(segprod[gn, t, k] for k = 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,9 +33,8 @@ function _add_ramp_eqs!(
|
||||
switch_off = model[:switch_off]
|
||||
switch_on = model[:switch_on]
|
||||
|
||||
for t in 1:model[:instance].time
|
||||
time_invariant =
|
||||
(t > 1) ? (abs(g.min_power[t] - g.min_power[t-1]) < 1e-7) : true
|
||||
for t = 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(
|
||||
@@ -48,10 +47,8 @@ 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]
|
||||
@@ -61,8 +58,7 @@ 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
|
||||
@@ -103,8 +99,7 @@ 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 in 1:model[:instance].time
|
||||
for k in 1:length(g.cost_segments)
|
||||
for t = 1:model[:instance].time
|
||||
for k = 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 in 1:model[:instance].time
|
||||
for t = 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,15 +21,13 @@ function _add_production_piecewise_linear_eqs!(
|
||||
is_on = model[:is_on]
|
||||
|
||||
K = length(g.cost_segments)
|
||||
for t in 1:model[:instance].time
|
||||
for t = 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 in 1:K)
|
||||
)
|
||||
eq_prod_above_def[gn, t] =
|
||||
@constraint(model, prod_above[gn, t] == sum(segprod[gn, t, k] for k = 1:K))
|
||||
|
||||
for k in 1:K
|
||||
for k = 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.
|
||||
@@ -48,11 +46,7 @@ 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 in 1:model[:instance].time
|
||||
for t = 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 in 1:model[:instance].time
|
||||
for t = 1:model[:instance].time
|
||||
if !g.must_run[t]
|
||||
# Link binary variables
|
||||
if t == 1
|
||||
@@ -51,10 +51,8 @@ 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 in 1:T
|
||||
for k in 1:K
|
||||
for t = 1:T
|
||||
for k = 1:K
|
||||
# Pbar^{k-1)
|
||||
Pbar0 =
|
||||
g.min_power[t] +
|
||||
(k > 1 ? sum(g.cost_segments[ell].mw[t] for ell in 1:k-1) : 0.0)
|
||||
(k > 1 ? sum(g.cost_segments[ell].mw[t] for ell = 1:k-1) : 0.0)
|
||||
# Pbar^k
|
||||
Pbar1 = g.cost_segments[k].mw[t] + Pbar0
|
||||
|
||||
@@ -61,8 +61,7 @@ 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
|
||||
@@ -70,8 +69,7 @@ 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)
|
||||
)
|
||||
|
||||
@@ -87,18 +85,12 @@ 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 in 1:K)
|
||||
)
|
||||
eq_prod_above_def[gn, t] =
|
||||
@constraint(model, prod_above[gn, t] == sum(segprod[gn, t, k] for k = 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,9 +32,8 @@ function _add_ramp_eqs!(
|
||||
switch_off = model[:switch_off]
|
||||
switch_on = model[:switch_on]
|
||||
|
||||
for t in 1:model[:instance].time
|
||||
time_invariant =
|
||||
(t > 1) ? (abs(g.min_power[t] - g.min_power[t-1]) < 1e-7) : true
|
||||
for t = 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
|
||||
@@ -60,8 +59,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]
|
||||
@@ -76,8 +75,7 @@ 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
|
||||
@@ -107,8 +105,7 @@ 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,30 +11,27 @@ function _add_startup_cost_eqs!(
|
||||
eq_startup_restrict = _init(model, :eq_startup_restrict)
|
||||
S = length(g.startup_categories)
|
||||
startup = model[:startup]
|
||||
for t in 1:model[:instance].time
|
||||
for t = 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 in 1:S)
|
||||
model[:switch_on][g.name, t] == sum(startup[g.name, t, s] for s = 1:S)
|
||||
)
|
||||
|
||||
for s in 1:S
|
||||
for s = 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,10 +14,8 @@ 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
|
||||
@@ -33,7 +31,7 @@ function _add_ramp_eqs!(
|
||||
switch_off = model[:switch_off]
|
||||
switch_on = model[:switch_on]
|
||||
|
||||
for t in 1:T
|
||||
for t = 1:T
|
||||
Pbar = g.max_power[t]
|
||||
if Pbar < 1e-7
|
||||
# Skip this time period if max power = 0
|
||||
@@ -54,13 +52,10 @@ 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 in 0:min(UT - 2, TRU, t - 1)
|
||||
i = 0:min(UT - 2, TRU, t - 1)
|
||||
)
|
||||
)
|
||||
|
||||
@@ -69,12 +64,10 @@ 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 in 0:min(UT - 1, TRU, t - 1)
|
||||
i = 0:min(UT - 1, TRU, t - 1)
|
||||
)
|
||||
)
|
||||
end
|
||||
@@ -89,13 +82,9 @@ 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 in 0:KSD
|
||||
) - sum(
|
||||
(Pbar - (SU + i * RU)) * switch_on[gn, t-i] for
|
||||
i in 0:KSU
|
||||
) - (
|
||||
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) - (
|
||||
(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 in 1:model[:instance].time
|
||||
for t = 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,37 +51,28 @@ function _add_ramp_eqs!(
|
||||
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)
|
||||
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)
|
||||
@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
|
||||
@@ -113,8 +104,7 @@ 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)
|
||||
@@ -123,8 +113,7 @@ 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,
|
||||
@@ -152,15 +141,13 @@ 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,13 +5,12 @@
|
||||
function _add_bus!(model::JuMP.Model, b::Bus)::Nothing
|
||||
net_injection = _init(model, :expr_net_injection)
|
||||
curtail = _init(model, :curtail)
|
||||
for t in 1:model[:instance].time
|
||||
for t = 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,13 +8,9 @@ function _add_transmission_line!(
|
||||
f::ShiftFactorsFormulation,
|
||||
)::Nothing
|
||||
overflow = _init(model, :overflow)
|
||||
for t in 1:model[:instance].time
|
||||
for t = 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,26 +2,18 @@
|
||||
# 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 in 1:model[:instance].time
|
||||
for t = 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,10 +13,7 @@ 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
|
||||
@@ -33,10 +30,7 @@ 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
|
||||
@@ -75,7 +69,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 in 1:n
|
||||
for i = 1:n
|
||||
lodf[:, i] *= 1.0 / (1.0 - lodf[i, i])
|
||||
lodf[i, i] = -1
|
||||
end
|
||||
|
||||
@@ -25,14 +25,7 @@ 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 in 1:T, b in model[:instance].buses
|
||||
for t = 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 in 1:T
|
||||
for t = 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 in 1:instance.time
|
||||
for t = 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,11 +46,7 @@ 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
|
||||
@@ -65,24 +61,20 @@ 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 in 1:instance.time
|
||||
for t = 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]
|
||||
)
|
||||
|
||||
@@ -91,10 +83,7 @@ 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 in 1:model[:instance].time
|
||||
for t = 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 in 1:model[:instance].time
|
||||
for t = 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 in 1:model[:instance].time
|
||||
for s in 1:length(g.startup_categories)
|
||||
for t = 1:model[:instance].time
|
||||
for s = 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 in 1:T
|
||||
for t = 1:T
|
||||
# Startup limit
|
||||
eq_startup_limit[g.name, t] = @constraint(
|
||||
model,
|
||||
@@ -96,16 +96,14 @@ 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
|
||||
@@ -121,7 +119,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 in 1:model[:instance].time
|
||||
for t = 1:model[:instance].time
|
||||
# Ramp up limit
|
||||
if t == 1
|
||||
if _is_initially_on(g) == 1
|
||||
@@ -151,8 +149,7 @@ 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
|
||||
@@ -165,18 +162,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 in 1:T
|
||||
for t = 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
|
||||
@@ -203,7 +200,7 @@ end
|
||||
|
||||
function _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing
|
||||
expr_net_injection = model[:expr_net_injection]
|
||||
for t in 1:model[:instance].time
|
||||
for t = 1:model[:instance].time
|
||||
# Add to net injection expression
|
||||
add_to_expression!(
|
||||
expr_net_injection[g.bus.name, t],
|
||||
|
||||
@@ -14,12 +14,10 @@ function fix!(model::JuMP.Model, solution::AbstractDict)::Nothing
|
||||
prod_above = model[:prod_above]
|
||||
reserve = model[:reserve]
|
||||
for g in instance.units
|
||||
for t in 1:T
|
||||
for t = 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,10 +2,7 @@
|
||||
# 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,11 +4,9 @@
|
||||
|
||||
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,11 +4,7 @@
|
||||
|
||||
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]
|
||||
@@ -18,13 +14,10 @@ function _find_violations(
|
||||
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 in 1:instance.time
|
||||
]
|
||||
overflow_values = [
|
||||
value(overflow[lm.name, t]) for lm in instance.lines,
|
||||
t in 1:instance.time
|
||||
value(net_injection[b.name, t]) for b in non_slack_buses, t = 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,
|
||||
@@ -35,10 +28,7 @@ function _find_violations(
|
||||
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
|
||||
|
||||
@@ -81,10 +71,8 @@ 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 in 1:T
|
||||
t => _ViolationFilter(max_total = max_per_period, max_per_line = max_per_line)
|
||||
for t = 1:T
|
||||
)
|
||||
|
||||
pre_flow::Array{Float64} = zeros(L, K) # pre_flow[lm, thread]
|
||||
@@ -92,35 +80,30 @@ 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 in 1:T
|
||||
]
|
||||
normal_limits::Array{Float64,2} =
|
||||
[l.normal_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
|
||||
]
|
||||
emergency_limits::Array{Float64,2} =
|
||||
[l.emergency_flow_limit[t] + overflow[l.offset, t] for l in instance.lines, t = 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 in 1:T
|
||||
@threads for t = 1:T
|
||||
k = threadid()
|
||||
|
||||
# Pre-contingency flows
|
||||
pre_flow[:, k] = isf * net_injections[:, t]
|
||||
|
||||
# Post-contingency flows
|
||||
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]
|
||||
for lc = 1:L, lm = 1:L
|
||||
post_flow[lm, lc, k] = pre_flow[lm, k] + pre_flow[lc, k] * lodf[lm, lc]
|
||||
end
|
||||
|
||||
# Pre-contingency violations
|
||||
for lm in 1:L
|
||||
for lm = 1:L
|
||||
pre_v[lm, k] = max(
|
||||
0.0,
|
||||
pre_flow[lm, k] - normal_limits[lm, t],
|
||||
@@ -129,7 +112,7 @@ function _find_violations(;
|
||||
end
|
||||
|
||||
# Post-contingency violations
|
||||
for lc in 1:L, lm in 1:L
|
||||
for lc = 1:L, lm = 1:L
|
||||
post_v[lm, lc, k] = max(
|
||||
0.0,
|
||||
post_flow[lm, lc, k] - emergency_limits[lm, t],
|
||||
@@ -138,7 +121,7 @@ function _find_violations(;
|
||||
end
|
||||
|
||||
# Offer pre-contingency violations
|
||||
for lm in 1:L
|
||||
for lm = 1:L
|
||||
if pre_v[lm, k] > 1e-5
|
||||
_offer(
|
||||
filters[t],
|
||||
@@ -153,7 +136,7 @@ function _find_violations(;
|
||||
end
|
||||
|
||||
# Offer post-contingency violations
|
||||
for lm in 1:L, lc in 1:L
|
||||
for lm = 1:L, lc = 1:L
|
||||
if post_v[lm, lc, k] > 1e-5 && is_vulnerable[lc]
|
||||
_offer(
|
||||
filters[t],
|
||||
@@ -169,7 +152,7 @@ function _find_violations(;
|
||||
end
|
||||
|
||||
violations = _Violation[]
|
||||
for t in 1:instance.time
|
||||
for t = 1:instance.time
|
||||
append!(violations, _query(filters[t]))
|
||||
end
|
||||
|
||||
|
||||
@@ -27,10 +27,7 @@ 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,43 +6,40 @@ 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 in 1:T]
|
||||
for b in collection
|
||||
b.name => [round(value(vars[b.name, t]), digits = 5) for t = 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 in 1:length(g.cost_segments)
|
||||
value(model[:segprod][g.name, t, k]) * g.cost_segments[k].cost[t] for
|
||||
k = 1:length(g.cost_segments)
|
||||
],
|
||||
) for t in 1:T
|
||||
) for t = 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 in 1:length(g.cost_segments)
|
||||
value(model[:segprod][g.name, t, k]) for k = 1:length(g.cost_segments)
|
||||
],
|
||||
) for t in 1:T
|
||||
) for t = 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 in 1:S
|
||||
) for t in 1:T
|
||||
g.startup_categories[s].cost * value(model[:startup][g.name, t, s]) for
|
||||
s = 1:S
|
||||
) for t = 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 (\$)"] =
|
||||
@@ -54,21 +51,19 @@ 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 in 1:instance.time
|
||||
round(value(model[:upflexiramp_shortfall][t]), digits = 5) : 0.0 for
|
||||
t = 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 = 1:instance.time
|
||||
)
|
||||
else
|
||||
# Report spinning reserve solutions only if both up-flexiramp and
|
||||
@@ -77,12 +72,11 @@ 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 = 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,16 +6,10 @@ 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 in 1:T
|
||||
for t = 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,10 +11,7 @@ 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
|
||||
@@ -31,11 +28,7 @@ function generate_initial_conditions!(
|
||||
@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,10 +117,7 @@ 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 *= α
|
||||
@@ -134,17 +131,11 @@ function _randomize_costs(
|
||||
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 in 1:instance.time
|
||||
for t = 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
|
||||
@@ -158,11 +149,9 @@ function _randomize_load_profile(
|
||||
)::Nothing
|
||||
# Generate new system load
|
||||
system_load = [1.0]
|
||||
for t in 2:instance.time
|
||||
for t = 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)
|
||||
@@ -172,7 +161,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 in 1:instance.time
|
||||
for t = 1:instance.time
|
||||
b.load[t] *= system_load[t] / prev_system_load[t]
|
||||
end
|
||||
end
|
||||
|
||||
@@ -93,9 +93,8 @@ 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 in 2:length(g.startup_categories)
|
||||
for s = 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 in 1:instance.time
|
||||
for t = 1:instance.time
|
||||
# Production cost curve should be convex
|
||||
for k in 2:length(g.cost_segments)
|
||||
for k = 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,10 +24,7 @@ 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)
|
||||
@@ -50,13 +47,12 @@ 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 in 1:instance.time
|
||||
for t = 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]
|
||||
@@ -90,11 +86,7 @@ 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
|
||||
|
||||
@@ -121,8 +113,7 @@ 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,
|
||||
@@ -136,11 +127,7 @@ 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
|
||||
|
||||
@@ -169,9 +156,7 @@ 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,
|
||||
@@ -201,7 +186,7 @@ function _validate_units(instance, solution; tol = 0.01)
|
||||
|
||||
# Calculate how much time the unit has been offline
|
||||
time_down = 0
|
||||
for k in 1:(t-1)
|
||||
for k = 1:(t-1)
|
||||
if !is_on[t-k]
|
||||
time_down += 1
|
||||
else
|
||||
@@ -235,7 +220,7 @@ function _validate_units(instance, solution; tol = 0.01)
|
||||
|
||||
# Calculate how much time the unit has been online
|
||||
time_up = 0
|
||||
for k in 1:(t-1)
|
||||
for k = 1:(t-1)
|
||||
if is_on[t-k]
|
||||
time_up += 1
|
||||
else
|
||||
@@ -288,7 +273,7 @@ end
|
||||
|
||||
function _validate_reserve_and_demand(instance, solution, tol = 0.01)
|
||||
err_count = 0
|
||||
for t in 1:instance.time
|
||||
for t = 1:instance.time
|
||||
load_curtail = 0
|
||||
fixed_load = sum(b.load[t] for b in instance.buses)
|
||||
ps_load = 0
|
||||
@@ -298,13 +283,10 @@ 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
|
||||
|
||||
@@ -324,7 +306,8 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
|
||||
|
||||
# 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)
|
||||
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_shortfall =
|
||||
@@ -362,8 +345,7 @@ 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
|
||||
|
||||
@@ -7,9 +7,8 @@ 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 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[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_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 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.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.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 in 1:4]
|
||||
@test unit.min_power_cost == [1400.0 for t in 1:4]
|
||||
@test unit.must_run == [false for t = 1:4]
|
||||
@test unit.min_power_cost == [1400.0 for t = 1:4]
|
||||
@test unit.min_uptime == 1
|
||||
@test unit.min_downtime == 1
|
||||
@test unit.provides_spinning_reserves == [true for t in 1:4]
|
||||
for t in 1:1
|
||||
@test unit.provides_spinning_reserves == [true for t = 1:4]
|
||||
for t = 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 in 1:4]
|
||||
@test unit.must_run == [false for t = 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 in 1:4]
|
||||
@test unit.min_power_cost == [0.0 for t in 1:4]
|
||||
@test unit.must_run == [true for t = 1:4]
|
||||
@test unit.min_power_cost == [0.0 for t = 1:4]
|
||||
@test unit.min_uptime == 1
|
||||
@test unit.min_downtime == 1
|
||||
@test unit.provides_spinning_reserves == [true for t in 1:4]
|
||||
for t in 1:4
|
||||
@test unit.provides_spinning_reserves == [true for t = 1:4]
|
||||
for t = 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 in 1:4]
|
||||
@test load.demand == [50.0 for t in 1:4]
|
||||
@test load.revenue == [100.0 for t = 1:4]
|
||||
@test load.demand == [50.0 for t = 1:4]
|
||||
@test instance.price_sensitive_loads_by_name["ps1"].name == "ps1"
|
||||
end
|
||||
|
||||
|
||||
@@ -23,6 +23,7 @@ function _small_test(formulation::Formulation)::Nothing
|
||||
instances = ["matpower/case118/2017-02-01", "test/case14"]
|
||||
for instance in instances
|
||||
# Should not crash
|
||||
@show "$(instance)"
|
||||
UnitCommitment.build_model(
|
||||
instance = UnitCommitment.read_benchmark(instance),
|
||||
formulation = formulation,
|
||||
@@ -58,17 +59,26 @@ function _test(formulation::Formulation)::Nothing
|
||||
end
|
||||
|
||||
@testset "formulations" begin
|
||||
@show "testset formulations"
|
||||
_test(Formulation())
|
||||
@show "ArrCon2000 ramping"
|
||||
_test(Formulation(ramping = ArrCon2000.Ramping()))
|
||||
|
||||
# _test(Formulation(ramping = DamKucRajAta2016.Ramping()))
|
||||
@show "MorLatRam2013 ramping"
|
||||
_test(
|
||||
Formulation(
|
||||
ramping = MorLatRam2013.Ramping(),
|
||||
startup_costs = MorLatRam2013.StartupCosts(),
|
||||
),
|
||||
)
|
||||
@show "PanGua2016 ramping"
|
||||
_test(Formulation(ramping = PanGua2016.Ramping()))
|
||||
@show "Gar1962 PwlCosts"
|
||||
_test(Formulation(pwl_costs = Gar1962.PwlCosts()))
|
||||
@show "CarArr2006 PwlCosts"
|
||||
_test(Formulation(pwl_costs = CarArr2006.PwlCosts()))
|
||||
@show "KnuOstWat2018 PwlCosts"
|
||||
_test(Formulation(pwl_costs = KnuOstWat2018.PwlCosts()))
|
||||
@show "formulations completed"
|
||||
end
|
||||
|
||||
@@ -11,6 +11,7 @@ UnitCommitment._setup_logger()
|
||||
const ENABLE_LARGE_TESTS = ("UCJL_LARGE_TESTS" in keys(ENV))
|
||||
|
||||
@testset "UnitCommitment" begin
|
||||
@show "running runtests.jl"
|
||||
include("usage.jl")
|
||||
@testset "import" begin
|
||||
include("import/egret_test.jl")
|
||||
@@ -27,6 +28,7 @@ const ENABLE_LARGE_TESTS = ("UCJL_LARGE_TESTS" in keys(ENV))
|
||||
include("solution/methods/XavQiuWanThi19/sensitivity_test.jl")
|
||||
end
|
||||
@testset "transform" begin
|
||||
@show "beginning transform"
|
||||
include("transform/initcond_test.jl")
|
||||
include("transform/slice_test.jl")
|
||||
@testset "randomize" begin
|
||||
|
||||
@@ -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 in 1:instance.time
|
||||
for line in instance.lines, t = 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 in 1:13, t in 1:instance.time]
|
||||
overflow = [0.0 for l in instance.lines, t in 1:instance.time]
|
||||
inj = [1000.0 for b = 1:13, t = 1:instance.time]
|
||||
overflow = [0.0 for l in instance.lines, t = 1:instance.time]
|
||||
violations = UnitCommitment._find_violations(
|
||||
instance = instance,
|
||||
net_injections = inj,
|
||||
|
||||
@@ -8,8 +8,7 @@ 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,10 +28,7 @@ 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)
|
||||
|
||||
@@ -5,6 +5,7 @@
|
||||
using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
|
||||
|
||||
@testset "slice" begin
|
||||
@show "beginning slice"
|
||||
instance = UnitCommitment.read_benchmark("test/case14")
|
||||
modified = UnitCommitment.slice(instance, 1:2)
|
||||
|
||||
@@ -35,7 +36,8 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
|
||||
@test length(ps.demand) == 2
|
||||
@test length(ps.revenue) == 2
|
||||
end
|
||||
|
||||
@show "beginning building model under slice"
|
||||
@show instance.reserves
|
||||
# Should be able to build model without errors
|
||||
optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
|
||||
model = UnitCommitment.build_model(
|
||||
|
||||
@@ -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 in 1:4
|
||||
for line in instance.lines, t = 1:4
|
||||
line.normal_flow_limit[t] = 10.0
|
||||
end
|
||||
optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
|
||||
|
||||
Reference in New Issue
Block a user