ANL-CEEESA
/
UnitCommitment.jl
mirror of https://github.com/ANL-CEEESA/UnitCommitment.jl.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

stochastic extension

Browse Source
pull/25/head
oyurdakul 3 years ago
parent c95b01dadf
commit 7e8a2ee026
39 changed files with 1197 additions and 1092 deletions
Show Stats Download Patch File Download Diff File

1
deps/formatter/format.jl vendored
Unescape View File

@ -1,4 +1,5 @@
using JuliaFormatter using JuliaFormatter
print(pwd())
format( format(
[ [
"../../src", "../../src",

2
src/format.jl
Unescape View File

@ -0,0 +1,2 @@
using JuliaFormatter
format(["../../src", "../../test", "../../benchmark/run.jl"], verbose = true)

101
src/instance/read.jl
Unescape View File

@ -44,25 +44,6 @@ function read_benchmark(
return UnitCommitment.read(filename) return UnitCommitment.read(filename)
end end
function read_scenarios(path::AbstractString)::UnitCommitmentInstance
scenario_paths = glob("*.json", path)
scenarios = Vector{UnitCommitmentScenario}()
total_number_of_scenarios = length(scenario_paths)
for (scenario_index, scenario_path) in enumerate(scenario_paths)
scenario = read_scenario(scenario_path,
total_number_of_scenarios = total_number_of_scenarios,
scenario_index = scenario_index)
push!(scenarios, scenario)
end
instance = UnitCommitmentInstance(
time = scenarios[1].time,
scenarios = scenarios
)
abs(sum(scenario.probability for scenario in instance.scenarios) - 1.0) <= 0.01 ||
error("scenario probabilities do not add up to one")
return instance
end
""" """
read(path::AbstractString)::UnitCommitmentInstance read(path::AbstractString)::UnitCommitmentInstance
@ -74,33 +55,54 @@ Read instance from a file. The file may be gzipped.
instance = UnitCommitment.read("/path/to/input.json.gz") instance = UnitCommitment.read("/path/to/input.json.gz")
``` ```
""" """
function read_scenario(path::AbstractString; total_number_of_scenarios = 1, scenario_index = 1)::UnitCommitmentScenario
if endswith(path, ".gz") function _repair_scenario_name_and_probability(
return _read(gzopen(path),total_number_of_scenarios, scenario_index) sc::UnitCommitmentScenario,
else path::String,
return _read(open(path), total_number_of_scenarios, scenario_index) number_of_scenarios::Int,
end )::UnitCommitmentScenario
sc.name !== nothing || (sc.name = first(split(last(split(path, "/")), ".")))
sc.probability !== nothing || (sc.probability = (1 / number_of_scenarios))
return sc
end end
function read(path::AbstractString; total_number_of_scenarios = 1, scenario_index = 1)::UnitCommitmentInstance function read(path)::UnitCommitmentInstance
if endswith(path, ".gz") scenarios = Vector{UnitCommitmentScenario}()
scenario = _read(gzopen(path),total_number_of_scenarios, scenario_index) if (endswith(path, ".gz") || endswith(path, ".json"))
endswith(path, ".gz") ? (scenario = _read(gzopen(path))) :
(scenario = _read(open(path)))
scenario = _repair_scenario_name_and_probability(scenario, "s1", 1)
scenarios = [scenario]
elseif typeof(path) == Vector{String}
number_of_scenarios = length(paths)
for scenario_path in path
if endswith(scenario_path, ".gz")
scenario = _read(gzopen(scenario_path))
elseif endswith(scenario_path, ".json")
scenario = _read(open(scenario_path))
else
error("Unsupported input format")
end
scenario = _repair_scenario_name_and_probability(
scenario,
scenario_path,
number_of_scenarios,
)
push!(scenarios, scenario)
end
else else
scenario = _read(open(path), total_number_of_scenarios, scenario_index) error("Unsupported input format")
end end
instance = UnitCommitmentInstance(
time = scenario.time, instance =
scenarios = [scenario] UnitCommitmentInstance(time = scenarios[1].time, scenarios = scenarios)
)
return instance return instance
end end
function _read(file::IO, total_number_of_scenarios::Int, scenario_index::Int)::UnitCommitmentScenario function _read(file::IO)::UnitCommitmentScenario
return _from_json( return _from_json(
JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)), JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)),
total_number_of_scenarios, )
scenario_index
)
end end
function _read_json(path::String)::OrderedDict function _read_json(path::String)::OrderedDict
@ -112,7 +114,7 @@ function _read_json(path::String)::OrderedDict
return JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)) return JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing))
end end
function _from_json(json, total_number_of_scenarios::Int, scenario_index::Int; repair = true)::UnitCommitmentScenario function _from_json(json; repair = true)::UnitCommitmentScenario
_migrate(json) _migrate(json)
units = Unit[] units = Unit[]
buses = Bus[] buses = Bus[]
@ -136,18 +138,8 @@ function _from_json(json, total_number_of_scenarios::Int, scenario_index::Int; r
error("Time step $time_step is not a divisor of 60") error("Time step $time_step is not a divisor of 60")
time_multiplier = 60 ÷ time_step time_multiplier = 60 ÷ time_step
T = time_horizon * time_multiplier T = time_horizon * time_multiplier
#####
probability = json["Parameters"]["Scenario probability"] probability = json["Parameters"]["Scenario probability"]
if probability === nothing
probability = (1 / total_number_of_scenarios)
end
scenario_name = json["Parameters"]["Scenario name"] scenario_name = json["Parameters"]["Scenario name"]
if scenario_name === nothing
scenario_name = "s$(scenario_index)"
end
######
name_to_bus = Dict{String,Bus}() name_to_bus = Dict{String,Bus}()
name_to_line = Dict{String,TransmissionLine}() name_to_line = Dict{String,TransmissionLine}()
name_to_unit = Dict{String,Unit}() name_to_unit = Dict{String,Unit}()
@ -164,15 +156,6 @@ function _from_json(json, total_number_of_scenarios::Int, scenario_index::Int; r
json["Parameters"]["Power balance penalty (\$/MW)"], json["Parameters"]["Power balance penalty (\$/MW)"],
default = [1000.0 for t in 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 in 1:T],
)
shortfall_penalty = timeseries(
json["Parameters"]["Reserve shortfall penalty (\$/MW)"],
default = [-1.0 for t in 1:T],
)
# Read buses # Read buses
for (bus_name, dict) in json["Buses"] for (bus_name, dict) in json["Buses"]
@ -368,13 +351,11 @@ function _from_json(json, total_number_of_scenarios::Int, scenario_index::Int; r
price_sensitive_loads = loads, price_sensitive_loads = loads,
reserves = reserves, reserves = reserves,
reserves_by_name = name_to_reserve, reserves_by_name = name_to_reserve,
# shortfall_penalty = shortfall_penalty,
# flexiramp_shortfall_penalty = flexiramp_shortfall_penalty,
time = T, time = T,
units_by_name = Dict(g.name => g for g in units), units_by_name = Dict(g.name => g for g in units),
units = units, units = units,
isf = spzeros(Float64, length(lines), length(buses) - 1), isf = spzeros(Float64, length(lines), length(buses) - 1),
lodf = spzeros(Float64, length(lines), length(lines)) lodf = spzeros(Float64, length(lines), length(lines)),
) )
if repair if repair
UnitCommitment.repair!(scenario) UnitCommitment.repair!(scenario)

23
src/instance/structs.jl
Unescape View File

@ -74,8 +74,8 @@ mutable struct PriceSensitiveLoad
end end
Base.@kwdef mutable struct UnitCommitmentScenario Base.@kwdef mutable struct UnitCommitmentScenario
name::String name::Any
probability::Float64 probability::Any
buses_by_name::Dict{AbstractString,Bus} buses_by_name::Dict{AbstractString,Bus}
buses::Vector{Bus} buses::Vector{Bus}
contingencies_by_name::Dict{AbstractString,Contingency} contingencies_by_name::Dict{AbstractString,Contingency}
@ -87,8 +87,6 @@ Base.@kwdef mutable struct UnitCommitmentScenario
price_sensitive_loads::Vector{PriceSensitiveLoad} price_sensitive_loads::Vector{PriceSensitiveLoad}
reserves::Vector{Reserve} reserves::Vector{Reserve}
reserves_by_name::Dict{AbstractString,Reserve} reserves_by_name::Dict{AbstractString,Reserve}
# shortfall_penalty::Vector{Float64}
# flexiramp_shortfall_penalty::Vector{Float64}
units_by_name::Dict{AbstractString,Unit} units_by_name::Dict{AbstractString,Unit}
units::Vector{Unit} units::Vector{Unit}
time::Int time::Int
@ -104,16 +102,13 @@ end
function Base.show(io::IO, instance::UnitCommitmentInstance) function Base.show(io::IO, instance::UnitCommitmentInstance)
print(io, "UnitCommitmentInstance(") print(io, "UnitCommitmentInstance(")
print(io, "$(length(instance.scenarios)) scenarios: ") print(io, "$(length(instance.scenarios)) scenarios: ")
for scenario in instance.scenarios for sc in instance.scenarios
print(io, "Scenario $(scenario.name): ") print(io, "Scenario $(sc.name): ")
print(io, "$(length(scenario.units)) units, ") print(io, "$(length(sc.units)) units, ")
print(io, "$(length(scenario.buses)) buses, ") print(io, "$(length(sc.buses)) buses, ")
print(io, "$(length(scenario.lines)) lines, ") print(io, "$(length(sc.lines)) lines, ")
print(io, "$(length(scenario.contingencies)) contingencies, ") print(io, "$(length(sc.contingencies)) contingencies, ")
print( print(io, "$(length(sc.price_sensitive_loads)) price sensitive loads, ")
io,
"$(length(scenario.price_sensitive_loads)) price sensitive loads, ",
)
end end
print(io, "$(instance.time) time steps") print(io, "$(instance.time) time steps")
print(io, ")") print(io, ")")

29
src/model/build.jl
Unescape View File

@ -78,26 +78,25 @@ function build_model(;
model[:obj] = AffExpr() model[:obj] = AffExpr()
model[:instance] = instance model[:instance] = instance
for g in instance.scenarios[1].units for g in instance.scenarios[1].units
_add_unit_first_stage!(model, g, formulation) _add_unit_commitment!(model, g, formulation)
end end
for scenario in instance.scenarios for sc in instance.scenarios
@info "Building scenario $(scenario.name) with" * @info "Building scenario $(sc.name) with" *
"probability $(scenario.probability)" "probability $(sc.probability)"
_setup_transmission(model, formulation.transmission, scenario) _setup_transmission(formulation.transmission, sc)
for l in scenario.lines for l in sc.lines
_add_transmission_line!(model, l, formulation.transmission, _add_transmission_line!(model, l, formulation.transmission, sc)
scenario)
end end
for b in scenario.buses for b in sc.buses
_add_bus!(model, b, scenario) _add_bus!(model, b, sc)
end end
for ps in scenario.price_sensitive_loads for ps in sc.price_sensitive_loads
_add_price_sensitive_load!(model, ps, scenario) _add_price_sensitive_load!(model, ps, sc)
end end
for g in scenario.units for g in sc.units
_add_unit_second_stage!(model, g, formulation, scenario) _add_unit_dispatch!(model, g, formulation, sc)
end end
_add_system_wide_eqs!(model, scenario) _add_system_wide_eqs!(model, sc)
end end
@objective(model, Min, model[:obj]) @objective(model, Min, model[:obj])
end end

5
src/model/formulations/ArrCon2000/ramp.jl
Unescape View File

@ -8,7 +8,7 @@ function _add_ramp_eqs!(
formulation_prod_vars::Gar1962.ProdVars, formulation_prod_vars::Gar1962.ProdVars,
formulation_ramping::ArrCon2000.Ramping, formulation_ramping::ArrCon2000.Ramping,
formulation_status_vars::Gar1962.StatusVars, formulation_status_vars::Gar1962.StatusVars,
sc::UnitCommitmentScenario sc::UnitCommitmentScenario,
)::Nothing )::Nothing
# TODO: Move upper case constants to model[:instance] # TODO: Move upper case constants to model[:instance]
RESERVES_WHEN_START_UP = true RESERVES_WHEN_START_UP = true
@ -74,7 +74,8 @@ function _add_ramp_eqs!(
# but the constraint below will force the unit to produce power # but the constraint below will force the unit to produce power
eq_ramp_down[sc.name, gn, t] = @constraint( eq_ramp_down[sc.name, gn, t] = @constraint(
model, model,
g.initial_power - (g.min_power[t] + prod_above[sc.name, gn, t]) <= RD g.initial_power -
(g.min_power[t] + prod_above[sc.name, gn, t]) <= RD
) )
end end
else else

10
src/model/formulations/CarArr2006/pwlcosts.jl
Unescape View File

@ -8,7 +8,7 @@ function _add_production_piecewise_linear_eqs!(
formulation_prod_vars::Gar1962.ProdVars, formulation_prod_vars::Gar1962.ProdVars,
formulation_pwl_costs::CarArr2006.PwlCosts, formulation_pwl_costs::CarArr2006.PwlCosts,
formulation_status_vars::StatusVarsFormulation, formulation_status_vars::StatusVarsFormulation,
sc::UnitCommitmentScenario sc::UnitCommitmentScenario,
)::Nothing )::Nothing
eq_prod_above_def = _init(model, :eq_prod_above_def) eq_prod_above_def = _init(model, :eq_prod_above_def)
eq_segprod_limit = _init(model, :eq_segprod_limit) eq_segprod_limit = _init(model, :eq_segprod_limit)
@ -34,13 +34,17 @@ function _add_production_piecewise_linear_eqs!(
# Also add this as an explicit upper bound on segprod to make the # Also add this as an explicit upper bound on segprod to make the
# solver's work a bit easier # solver's work a bit easier
set_upper_bound(segprod[sc.name, gn, t, k], g.cost_segments[k].mw[t]) set_upper_bound(
segprod[sc.name, gn, t, k],
g.cost_segments[k].mw[t],
)
# Definition of production # Definition of production
# Equation (43) in Kneuven et al. (2020) # Equation (43) in Kneuven et al. (2020)
eq_prod_above_def[sc.name, gn, t] = @constraint( eq_prod_above_def[sc.name, gn, t] = @constraint(
model, model,
prod_above[sc.name, gn, t] == sum(segprod[sc.name, gn, t, k] for k in 1:K) prod_above[sc.name, gn, t] ==
sum(segprod[sc.name, gn, t, k] for k in 1:K)
) )
# Objective function # Objective function

5
src/model/formulations/DamKucRajAta2016/ramp.jl
Unescape View File

@ -8,7 +8,7 @@ function _add_ramp_eqs!(
formulation_prod_vars::Gar1962.ProdVars, formulation_prod_vars::Gar1962.ProdVars,
formulation_ramping::DamKucRajAta2016.Ramping, formulation_ramping::DamKucRajAta2016.Ramping,
formulation_status_vars::Gar1962.StatusVars, formulation_status_vars::Gar1962.StatusVars,
sc::UnitCommitmentScenario sc::UnitCommitmentScenario,
)::Nothing )::Nothing
# TODO: Move upper case constants to model[:instance] # TODO: Move upper case constants to model[:instance]
RESERVES_WHEN_START_UP = true RESERVES_WHEN_START_UP = true
@ -66,7 +66,8 @@ function _add_ramp_eqs!(
elseif (t == 1 && is_initially_on) || (t > 1 && !time_invariant) elseif (t == 1 && is_initially_on) || (t > 1 && !time_invariant)
if t > 1 if t > 1
min_prod_last_period = min_prod_last_period =
prod_above[sc.name, gn, t-1] + g.min_power[t-1] * is_on[gn, t-1] prod_above[sc.name, gn, t-1] +
g.min_power[t-1] * is_on[gn, t-1]
else else
min_prod_last_period = max(g.initial_power, 0.0) min_prod_last_period = max(g.initial_power, 0.0)
end end

11
src/model/formulations/Gar1962/prod.jl
Unescape View File

@ -6,7 +6,7 @@ function _add_production_vars!(
model::JuMP.Model, model::JuMP.Model,
g::Unit, g::Unit,
formulation_prod_vars::Gar1962.ProdVars, formulation_prod_vars::Gar1962.ProdVars,
sc::UnitCommitmentScenario sc::UnitCommitmentScenario,
)::Nothing )::Nothing
prod_above = _init(model, :prod_above) prod_above = _init(model, :prod_above)
segprod = _init(model, :segprod) segprod = _init(model, :segprod)
@ -23,7 +23,7 @@ function _add_production_limit_eqs!(
model::JuMP.Model, model::JuMP.Model,
g::Unit, g::Unit,
formulation_prod_vars::Gar1962.ProdVars, formulation_prod_vars::Gar1962.ProdVars,
sc::UnitCommitmentScenario sc::UnitCommitmentScenario,
)::Nothing )::Nothing
eq_prod_limit = _init(model, :eq_prod_limit) eq_prod_limit = _init(model, :eq_prod_limit)
is_on = model[:is_on] is_on = model[:is_on]
@ -34,10 +34,6 @@ function _add_production_limit_eqs!(
# Objective function terms for production costs # Objective function terms for production costs
# Part of (69) of Kneuven et al. (2020) as C^R_g * u_g(t) term # Part of (69) of Kneuven et al. (2020) as C^R_g * u_g(t) term
### Moving this term to another function
# add_to_expression!(model[:obj], is_on[gn, t], g.min_power_cost[t])
###
# Production limit # Production limit
# Equation (18) in Kneuven et al. (2020) # Equation (18) in Kneuven et al. (2020)
# as \bar{p}_g(t) \le \bar{P}_g u_g(t) # as \bar{p}_g(t) \le \bar{P}_g u_g(t)
@ -49,7 +45,8 @@ function _add_production_limit_eqs!(
end end
eq_prod_limit[sc.name, gn, t] = @constraint( eq_prod_limit[sc.name, gn, t] = @constraint(
model, model,
prod_above[sc.name, gn, t] + reserve[t] <= power_diff * is_on[gn, t] prod_above[sc.name, gn, t] + reserve[t] <=
power_diff * is_on[gn, t]
) )
end end
end end

13
src/model/formulations/Gar1962/pwlcosts.jl
Unescape View File

@ -8,7 +8,7 @@ function _add_production_piecewise_linear_eqs!(
formulation_prod_vars::Gar1962.ProdVars, formulation_prod_vars::Gar1962.ProdVars,
formulation_pwl_costs::Gar1962.PwlCosts, formulation_pwl_costs::Gar1962.PwlCosts,
formulation_status_vars::Gar1962.StatusVars, formulation_status_vars::Gar1962.StatusVars,
sc::UnitCommitmentScenario sc::UnitCommitmentScenario,
)::Nothing )::Nothing
eq_prod_above_def = _init(model, :eq_prod_above_def) eq_prod_above_def = _init(model, :eq_prod_above_def)
eq_segprod_limit = _init(model, :eq_segprod_limit) eq_segprod_limit = _init(model, :eq_segprod_limit)
@ -27,7 +27,8 @@ function _add_production_piecewise_linear_eqs!(
# Equation (43) in Kneuven et al. (2020) # Equation (43) in Kneuven et al. (2020)
eq_prod_above_def[sc.name, gn, t] = @constraint( eq_prod_above_def[sc.name, gn, t] = @constraint(
model, model,
prod_above[sc.name, gn, t] == sum(segprod[sc.name, gn, t, k] for k in 1:K) prod_above[sc.name, gn, t] ==
sum(segprod[sc.name, gn, t, k] for k in 1:K)
) )
for k in 1:K for k in 1:K
@ -40,12 +41,16 @@ function _add_production_piecewise_linear_eqs!(
# that segment* # that segment*
eq_segprod_limit[sc.name, gn, t, k] = @constraint( eq_segprod_limit[sc.name, gn, t, k] = @constraint(
model, model,
segprod[sc.name, gn, t, k] <= g.cost_segments[k].mw[t] * is_on[gn, t] segprod[sc.name, gn, t, k] <=
g.cost_segments[k].mw[t] * is_on[gn, t]
) )
# Also add this as an explicit upper bound on segprod to make the # Also add this as an explicit upper bound on segprod to make the
# solver's work a bit easier # solver's work a bit easier
set_upper_bound(segprod[sc.name, gn, t, k], g.cost_segments[k].mw[t]) set_upper_bound(
segprod[sc.name, gn, t, k],
g.cost_segments[k].mw[t],
)
# Objective function # Objective function
# Equation (44) in Kneuven et al. (2020) # Equation (44) in Kneuven et al. (2020)

10
src/model/formulations/KnuOstWat2018/pwlcosts.jl
Unescape View File

@ -8,7 +8,7 @@ function _add_production_piecewise_linear_eqs!(
formulation_prod_vars::Gar1962.ProdVars, formulation_prod_vars::Gar1962.ProdVars,
formulation_pwl_costs::KnuOstWat2018.PwlCosts, formulation_pwl_costs::KnuOstWat2018.PwlCosts,
formulation_status_vars::Gar1962.StatusVars, formulation_status_vars::Gar1962.StatusVars,
sc::UnitCommitmentScenario sc::UnitCommitmentScenario,
)::Nothing )::Nothing
eq_prod_above_def = _init(model, :eq_prod_above_def) eq_prod_above_def = _init(model, :eq_prod_above_def)
eq_segprod_limit_a = _init(model, :eq_segprod_limit_a) eq_segprod_limit_a = _init(model, :eq_segprod_limit_a)
@ -90,7 +90,8 @@ function _add_production_piecewise_linear_eqs!(
# Equation (43) in Kneuven et al. (2020) # Equation (43) in Kneuven et al. (2020)
eq_prod_above_def[sc.name, gn, t] = @constraint( eq_prod_above_def[sc.name, gn, t] = @constraint(
model, model,
prod_above[sc.name, gn, t] == sum(segprod[sc.name, gn, t, k] for k in 1:K) prod_above[sc.name, gn, t] ==
sum(segprod[sc.name, gn, t, k] for k in 1:K)
) )
# Objective function # Objective function
@ -103,7 +104,10 @@ function _add_production_piecewise_linear_eqs!(
# Also add an explicit upper bound on segprod to make the solver's # Also add an explicit upper bound on segprod to make the solver's
# work a bit easier # work a bit easier
set_upper_bound(segprod[sc.name, gn, t, k], g.cost_segments[k].mw[t]) set_upper_bound(
segprod[sc.name, gn, t, k],
g.cost_segments[k].mw[t],
)
end end
end end
end end

8
src/model/formulations/MorLatRam2013/ramp.jl
Unescape View File

@ -8,7 +8,7 @@ function _add_ramp_eqs!(
formulation_prod_vars::Gar1962.ProdVars, formulation_prod_vars::Gar1962.ProdVars,
formulation_ramping::MorLatRam2013.Ramping, formulation_ramping::MorLatRam2013.Ramping,
formulation_status_vars::Gar1962.StatusVars, formulation_status_vars::Gar1962.StatusVars,
sc::UnitCommitmentScenario sc::UnitCommitmentScenario,
)::Nothing )::Nothing
# TODO: Move upper case constants to model[:instance] # TODO: Move upper case constants to model[:instance]
RESERVES_WHEN_START_UP = true RESERVES_WHEN_START_UP = true
@ -65,7 +65,8 @@ function _add_ramp_eqs!(
reserve[t] : 0.0 reserve[t] : 0.0
) )
min_prod_last_period = min_prod_last_period =
g.min_power[t-1] * is_on[gn, t-1] + prod_above[sc.name, gn, t-1] g.min_power[t-1] * is_on[gn, t-1] +
prod_above[sc.name, gn, t-1]
eq_ramp_up[gn, t] = @constraint( eq_ramp_up[gn, t] = @constraint(
model, model,
max_prod_this_period - min_prod_last_period <= max_prod_this_period - min_prod_last_period <=
@ -93,7 +94,8 @@ function _add_ramp_eqs!(
# but the constraint below will force the unit to produce power # but the constraint below will force the unit to produce power
eq_ramp_down[sc.name, gn, t] = @constraint( eq_ramp_down[sc.name, gn, t] = @constraint(
model, model,
g.initial_power - (g.min_power[t] + prod_above[sc.name, gn, t]) <= RD g.initial_power -
(g.min_power[t] + prod_above[sc.name, gn, t]) <= RD
) )
end end
else else

21
src/model/formulations/PanGua2016/ramp.jl
Unescape View File

@ -8,7 +8,7 @@ function _add_ramp_eqs!(
formulation_prod_vars::Gar1962.ProdVars, formulation_prod_vars::Gar1962.ProdVars,
formulation_ramping::PanGua2016.Ramping, formulation_ramping::PanGua2016.Ramping,
formulation_status_vars::Gar1962.StatusVars, formulation_status_vars::Gar1962.StatusVars,
sc::UnitCommitmentScenario sc::UnitCommitmentScenario,
)::Nothing )::Nothing
# TODO: Move upper case constants to model[:instance] # TODO: Move upper case constants to model[:instance]
RESERVES_WHEN_SHUT_DOWN = true RESERVES_WHEN_SHUT_DOWN = true
@ -68,16 +68,17 @@ function _add_ramp_eqs!(
if UT - 2 < TRU if UT - 2 < TRU
# Equation (40) in Kneuven et al. (2020) # Equation (40) in Kneuven et al. (2020)
# Covers an additional time period of the ramp-up trajectory, compared to (38) # Covers an additional time period of the ramp-up trajectory, compared to (38)
eq_prod_limit_ramp_up_extra_period[sc.name, gn, t] = @constraint( eq_prod_limit_ramp_up_extra_period[sc.name, gn, t] =
model, @constraint(
prod_above[sc.name, gn, t] + model,
g.min_power[t] * is_on[gn, t] + prod_above[sc.name, gn, t] +
reserve[t] <= g.min_power[t] * is_on[gn, t] +
Pbar * is_on[gn, t] - sum( reserve[t] <=
(Pbar - (SU + i * RU)) * switch_on[gn, t-i] for Pbar * is_on[gn, t] - sum(
i in 0:min(UT - 1, TRU, t - 1) (Pbar - (SU + i * RU)) * switch_on[gn, t-i] for
i in 0:min(UT - 1, TRU, t - 1)
)
) )
)
end end
# Add in shutdown trajectory if KSD >= 0 (else this is dominated by (38)) # Add in shutdown trajectory if KSD >= 0 (else this is dominated by (38))

57
src/model/formulations/WanHob2016/ramp.jl
Unescape View File

@ -8,7 +8,7 @@ function _add_ramp_eqs!(
::Gar1962.ProdVars, ::Gar1962.ProdVars,
::WanHob2016.Ramping, ::WanHob2016.Ramping,
::Gar1962.StatusVars, ::Gar1962.StatusVars,
sc::UnitCommitmentScenario sc::UnitCommitmentScenario,
)::Nothing )::Nothing
is_initially_on = (g.initial_status > 0) is_initially_on = (g.initial_status > 0)
SU = g.startup_limit SU = g.startup_limit
@ -30,7 +30,7 @@ function _add_ramp_eqs!(
error("Each generator may only provide one flexiramp reserve") error("Each generator may only provide one flexiramp reserve")
end end
for r in g.reserves for r in g.reserves
if r.type !== "up-frp" && r.type !== "down-frp" if r.type !== "flexiramp"
error( error(
"This formulation only supports flexiramp reserves, not $(r.type)", "This formulation only supports flexiramp reserves, not $(r.type)",
) )
@ -39,21 +39,23 @@ function _add_ramp_eqs!(
for t in 1:model[:instance].time for t in 1:model[:instance].time
@constraint( @constraint(
model, model,
prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t]) <= mfg[sc.name, rn, gn, t] prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t]) <=
mfg[sc.name, gn, t]
) # Eq. (19) in Wang & Hobbs (2016) ) # Eq. (19) in Wang & Hobbs (2016)
@constraint(model, mfg[sc.name, rn, gn, t] <= is_on[gn, t] * maxp[t]) # Eq. (22) in Wang & Hobbs (2016) @constraint(model, mfg[sc.name, gn, t] <= is_on[gn, t] * maxp[t]) # Eq. (22) in Wang & Hobbs (2016)
if t != model[:instance].time if t != model[:instance].time
@constraint( @constraint(
model, model,
minp[t] * (is_on[gn, t+1] + is_on[gn, t] - 1) <= minp[t] * (is_on[gn, t+1] + is_on[gn, t] - 1) <=
prod_above[sc.name, gn, t] - dwflexiramp[sc.name, rn, gn, t] + prod_above[sc.name, gn, t] -
(is_on[gn, t] * minp[t]) dwflexiramp[sc.name, rn, gn, t] + (is_on[gn, t] * minp[t])
) # first inequality of Eq. (20) in Wang & Hobbs (2016) ) # first inequality of Eq. (20) in Wang & Hobbs (2016)
@constraint( @constraint(
model, model,
prod_above[sc.name, gn, t] - dwflexiramp[sc.name, rn, gn, t] + prod_above[sc.name, gn, t] -
dwflexiramp[sc.name, rn, gn, t] +
(is_on[gn, t] * minp[t]) <= (is_on[gn, t] * minp[t]) <=
mfg[sc.name, rn, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1])) mfg[sc.name, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1]))
) # second inequality of Eq. (20) in Wang & Hobbs (2016) ) # second inequality of Eq. (20) in Wang & Hobbs (2016)
@constraint( @constraint(
model, model,
@ -67,12 +69,12 @@ function _add_ramp_eqs!(
prod_above[sc.name, gn, t] + prod_above[sc.name, gn, t] +
upflexiramp[sc.name, rn, gn, t] + upflexiramp[sc.name, rn, gn, t] +
(is_on[gn, t] * minp[t]) <= (is_on[gn, t] * minp[t]) <=
mfg[sc.name, rn, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1])) mfg[sc.name, gn, t+1] + (maxp[t] * (1 - is_on[gn, t+1]))
) # second inequality of Eq. (21) in Wang & Hobbs (2016) ) # second inequality of Eq. (21) in Wang & Hobbs (2016)
if t != 1 if t != 1
@constraint( @constraint(
model, model,
mfg[sc.name, rn, gn, t] <= mfg[sc.name, gn, t] <=
prod_above[sc.name, gn, t-1] + prod_above[sc.name, gn, t-1] +
(is_on[gn, t-1] * minp[t]) + (is_on[gn, t-1] * minp[t]) +
(RU * is_on[gn, t-1]) + (RU * is_on[gn, t-1]) +
@ -81,8 +83,13 @@ function _add_ramp_eqs!(
) # Eq. (23) in Wang & Hobbs (2016) ) # Eq. (23) in Wang & Hobbs (2016)
@constraint( @constraint(
model, model,
(prod_above[sc.name, gn, t-1] + (is_on[gn, t-1] * minp[t])) - (
(prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t])) <= prod_above[sc.name, gn, t-1] +
(is_on[gn, t-1] * minp[t])
) - (
prod_above[sc.name, gn, t] +
(is_on[gn, t] * minp[t])
) <=
RD * is_on[gn, t] + RD * is_on[gn, t] +
SD * (is_on[gn, t-1] - is_on[gn, t]) + SD * (is_on[gn, t-1] - is_on[gn, t]) +
maxp[t] * (1 - is_on[gn, t-1]) maxp[t] * (1 - is_on[gn, t-1])
@ -90,7 +97,7 @@ function _add_ramp_eqs!(
else else
@constraint( @constraint(
model, model,
mfg[sc.name, rn, gn, t] <= mfg[sc.name, gn, t] <=
initial_power + initial_power +
(RU * is_initially_on) + (RU * is_initially_on) +
(SU * (is_on[gn, t] - is_initially_on)) + (SU * (is_on[gn, t] - is_initially_on)) +
@ -98,8 +105,10 @@ function _add_ramp_eqs!(
) # Eq. (23) in Wang & Hobbs (2016) for the first time period ) # Eq. (23) in Wang & Hobbs (2016) for the first time period
@constraint( @constraint(
model, model,
initial_power - initial_power - (
(prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t])) <= prod_above[sc.name, gn, t] +
(is_on[gn, t] * minp[t])
) <=
RD * is_on[gn, t] + RD * is_on[gn, t] +
SD * (is_initially_on - is_on[gn, t]) + SD * (is_initially_on - is_on[gn, t]) +
maxp[t] * (1 - is_initially_on) maxp[t] * (1 - is_initially_on)
@ -107,7 +116,7 @@ function _add_ramp_eqs!(
end end
@constraint( @constraint(
model, model,
mfg[sc.name, rn, gn, t] <= mfg[sc.name, gn, t] <=
(SD * (is_on[gn, t] - is_on[gn, t+1])) + (SD * (is_on[gn, t] - is_on[gn, t+1])) +
(maxp[t] * is_on[gn, t+1]) (maxp[t] * is_on[gn, t+1])
) # Eq. (24) in Wang & Hobbs (2016) ) # Eq. (24) in Wang & Hobbs (2016)
@ -115,7 +124,8 @@ function _add_ramp_eqs!(
model, model,
-RD * is_on[gn, t+1] - -RD * is_on[gn, t+1] -
SD * (is_on[gn, t] - is_on[gn, t+1]) - SD * (is_on[gn, t] - is_on[gn, t+1]) -
maxp[t] * (1 - is_on[gn, t]) <= upflexiramp[sc.name, rn, gn, t] maxp[t] * (1 - is_on[gn, t]) <=
upflexiramp[sc.name, rn, gn, t]
) # first inequality of Eq. (26) in Wang & Hobbs (2016) ) # first inequality of Eq. (26) in Wang & Hobbs (2016)
@constraint( @constraint(
model, model,
@ -127,7 +137,8 @@ function _add_ramp_eqs!(
@constraint( @constraint(
model, model,
-RU * is_on[gn, t] - SU * (is_on[gn, t+1] - is_on[gn, t]) - -RU * is_on[gn, t] - SU * (is_on[gn, t+1] - is_on[gn, t]) -
maxp[t] * (1 - is_on[gn, t+1]) <= dwflexiramp[sc.name, rn, gn, t] maxp[t] * (1 - is_on[gn, t+1]) <=
dwflexiramp[sc.name, rn, gn, t]
) # first inequality of Eq. (27) in Wang & Hobbs (2016) ) # first inequality of Eq. (27) in Wang & Hobbs (2016)
@constraint( @constraint(
model, model,
@ -147,7 +158,8 @@ function _add_ramp_eqs!(
) # second inequality of Eq. (28) in Wang & Hobbs (2016) ) # second inequality of Eq. (28) in Wang & Hobbs (2016)
@constraint( @constraint(
model, model,
-maxp[t] * is_on[gn, t+1] <= dwflexiramp[sc.name, rn, gn, t] -maxp[t] * is_on[gn, t+1] <=
dwflexiramp[sc.name, rn, gn, t]
) # first inequality of Eq. (29) in Wang & Hobbs (2016) ) # first inequality of Eq. (29) in Wang & Hobbs (2016)
@constraint( @constraint(
model, model,
@ -157,7 +169,7 @@ function _add_ramp_eqs!(
else else
@constraint( @constraint(
model, model,
mfg[sc.name, rn, gn, t] <= mfg[sc.name, gn, t] <=
prod_above[sc.name, gn, t-1] + prod_above[sc.name, gn, t-1] +
(is_on[gn, t-1] * minp[t]) + (is_on[gn, t-1] * minp[t]) +
(RU * is_on[gn, t-1]) + (RU * is_on[gn, t-1]) +
@ -166,7 +178,10 @@ function _add_ramp_eqs!(
) # Eq. (23) in Wang & Hobbs (2016) for the last time period ) # Eq. (23) in Wang & Hobbs (2016) for the last time period
@constraint( @constraint(
model, model,
(prod_above[sc.name, gn, t-1] + (is_on[gn, t-1] * minp[t])) - (
prod_above[sc.name, gn, t-1] +
(is_on[gn, t-1] * minp[t])
) -
(prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t])) <= (prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t])) <=
RD * is_on[gn, t] + RD * is_on[gn, t] +
SD * (is_on[gn, t-1] - is_on[gn, t]) + SD * (is_on[gn, t-1] - is_on[gn, t]) +

12
src/model/formulations/base/bus.jl
Unescape View File

@ -2,7 +2,11 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details. # Released under the modified BSD license. See COPYING.md for more details.
function _add_bus!(model::JuMP.Model, b::Bus, sc::UnitCommitmentScenario)::Nothing function _add_bus!(
model::JuMP.Model,
b::Bus,
sc::UnitCommitmentScenario,
)::Nothing
net_injection = _init(model, :expr_net_injection) net_injection = _init(model, :expr_net_injection)
curtail = _init(model, :curtail) curtail = _init(model, :curtail)
for t in 1:model[:instance].time for t in 1:model[:instance].time
@ -13,7 +17,11 @@ function _add_bus!(model::JuMP.Model, b::Bus, sc::UnitCommitmentScenario)::Nothi
curtail[sc.name, b.name, t] = curtail[sc.name, b.name, t] =
@variable(model, lower_bound = 0, upper_bound = b.load[t]) @variable(model, lower_bound = 0, upper_bound = b.load[t])
add_to_expression!(net_injection[sc.name, b.name, t], curtail[sc.name, b.name, t], 1.0) add_to_expression!(
net_injection[sc.name, b.name, t],
curtail[sc.name, b.name, t],
1.0,
)
add_to_expression!( add_to_expression!(
model[:obj], model[:obj],
curtail[sc.name, b.name, t], curtail[sc.name, b.name, t],

20
src/model/formulations/base/line.jl
Unescape View File

@ -6,7 +6,7 @@ function _add_transmission_line!(
model::JuMP.Model, model::JuMP.Model,
lm::TransmissionLine, lm::TransmissionLine,
f::ShiftFactorsFormulation, f::ShiftFactorsFormulation,
sc::UnitCommitmentScenario sc::UnitCommitmentScenario,
)::Nothing )::Nothing
overflow = _init(model, :overflow) overflow = _init(model, :overflow)
for t in 1:model[:instance].time for t in 1:model[:instance].time
@ -21,30 +21,28 @@ function _add_transmission_line!(
end end
function _setup_transmission( function _setup_transmission(
model::JuMP.Model,
formulation::ShiftFactorsFormulation, formulation::ShiftFactorsFormulation,
scenario::UnitCommitmentScenario sc::UnitCommitmentScenario,
)::Nothing )::Nothing
instance = model[:instance]
isf = formulation.precomputed_isf isf = formulation.precomputed_isf
lodf = formulation.precomputed_lodf lodf = formulation.precomputed_lodf
if length(scenario.buses) == 1 if length(sc.buses) == 1
isf = zeros(0, 0) isf = zeros(0, 0)
lodf = zeros(0, 0) lodf = zeros(0, 0)
elseif isf === nothing elseif isf === nothing
@info "Computing injection shift factors..." @info "Computing injection shift factors..."
time_isf = @elapsed begin time_isf = @elapsed begin
isf = UnitCommitment._injection_shift_factors( isf = UnitCommitment._injection_shift_factors(
lines = scenario.lines, buses = sc.buses,
buses = scenario.buses, lines = sc.lines,
) )
end end
@info @sprintf("Computed ISF in %.2f seconds", time_isf) @info @sprintf("Computed ISF in %.2f seconds", time_isf)
@info "Computing line outage factors..." @info "Computing line outage factors..."
time_lodf = @elapsed begin time_lodf = @elapsed begin
lodf = UnitCommitment._line_outage_factors( lodf = UnitCommitment._line_outage_factors(
lines = scenario.lines, buses = sc.buses,
buses = scenario.buses, lines = sc.lines,
isf = isf, isf = isf,
) )
end end
@ -57,7 +55,7 @@ function _setup_transmission(
isf[abs.(isf).<formulation.isf_cutoff] .= 0 isf[abs.(isf).<formulation.isf_cutoff] .= 0
lodf[abs.(lodf).<formulation.lodf_cutoff] .= 0 lodf[abs.(lodf).<formulation.lodf_cutoff] .= 0
end end
scenario.isf = isf sc.isf = isf
scenario.lodf = lodf sc.lodf = lodf
return return
end end

9
src/model/formulations/base/psload.jl
Unescape View File

@ -5,7 +5,7 @@
function _add_price_sensitive_load!( function _add_price_sensitive_load!(
model::JuMP.Model, model::JuMP.Model,
ps::PriceSensitiveLoad, ps::PriceSensitiveLoad,
sc::UnitCommitmentScenario sc::UnitCommitmentScenario,
)::Nothing )::Nothing
loads = _init(model, :loads) loads = _init(model, :loads)
net_injection = _init(model, :expr_net_injection) net_injection = _init(model, :expr_net_injection)
@ -15,8 +15,11 @@ function _add_price_sensitive_load!(
@variable(model, lower_bound = 0, upper_bound = ps.demand[t]) @variable(model, lower_bound = 0, upper_bound = ps.demand[t])
# Objective function terms # Objective function terms
add_to_expression!(model[:obj], loads[ps.name, t], add_to_expression!(
-ps.revenue[t] * sc.probability) model[:obj],
loads[sc.name, ps.name, t],
-ps.revenue[t] * sc.probability,
)
# Net injection # Net injection
add_to_expression!( add_to_expression!(

95
src/model/formulations/base/system.jl
Unescape View File

@ -2,22 +2,30 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details. # Released under the modified BSD license. See COPYING.md for more details.
function _add_system_wide_eqs!(model::JuMP.Model, sc::UnitCommitmentScenario)::Nothing function _add_system_wide_eqs!(
model::JuMP.Model,
sc::UnitCommitmentScenario,
)::Nothing
_add_net_injection_eqs!(model, sc) _add_net_injection_eqs!(model, sc)
_add_spinning_reserve_eqs!(model, sc) _add_spinning_reserve_eqs!(model, sc)
_add_flexiramp_reserve_eqs!(model, sc) _add_flexiramp_reserve_eqs!(model, sc)
return return
end end
function _add_net_injection_eqs!(model::JuMP.Model, sc::UnitCommitmentScenario)::Nothing function _add_net_injection_eqs!(
model::JuMP.Model,
sc::UnitCommitmentScenario,
)::Nothing
T = model[:instance].time T = model[:instance].time
net_injection = _init(model, :net_injection) net_injection = _init(model, :net_injection)
eq_net_injection = _init(model, :eq_net_injection) eq_net_injection = _init(model, :eq_net_injection)
eq_power_balance = _init(model, :eq_power_balance) eq_power_balance = _init(model, :eq_power_balance)
for t in 1:T, b in sc.buses for t in 1:T, b in sc.buses
n = net_injection[sc.name, b.name, t] = @variable(model) n = net_injection[sc.name, b.name, t] = @variable(model)
eq_net_injection[sc.name, b.name, t] = eq_net_injection[sc.name, b.name, t] = @constraint(
@constraint(model, -n + model[:expr_net_injection][sc.name, b.name, t] == 0) model,
-n + model[:expr_net_injection][sc.name, b.name, t] == 0
)
end end
for t in 1:T for t in 1:T
eq_power_balance[sc.name, t] = @constraint( eq_power_balance[sc.name, t] = @constraint(
@ -28,7 +36,10 @@ function _add_net_injection_eqs!(model::JuMP.Model, sc::UnitCommitmentScenario):
return return
end end
function _add_spinning_reserve_eqs!(model::JuMP.Model, sc::UnitCommitmentScenario)::Nothing function _add_spinning_reserve_eqs!(
model::JuMP.Model,
sc::UnitCommitmentScenario,
)::Nothing
T = model[:instance].time T = model[:instance].time
eq_min_spinning_reserve = _init(model, :eq_min_spinning_reserve) eq_min_spinning_reserve = _init(model, :eq_min_spinning_reserve)
for r in sc.reserves for r in sc.reserves
@ -40,8 +51,11 @@ function _add_spinning_reserve_eqs!(model::JuMP.Model, sc::UnitCommitmentScenari
# from Carrión and Arroyo (2006) and Ostrowski et al. (2012) # from Carrión and Arroyo (2006) and Ostrowski et al. (2012)
eq_min_spinning_reserve[sc.name, r.name, t] = @constraint( eq_min_spinning_reserve[sc.name, r.name, t] = @constraint(
model, model,
sum(model[:reserve][sc.name, r.name, g.name, t] for g in r.units) + sum(
model[:reserve_shortfall][sc.name, r.name, t] >= r.amount[t] model[:reserve][sc.name, r.name, g.name, t] for
g in r.units
) + model[:reserve_shortfall][sc.name, r.name, t] >=
r.amount[t]
) )
# Account for shortfall contribution to objective # Account for shortfall contribution to objective
@ -57,7 +71,10 @@ function _add_spinning_reserve_eqs!(model::JuMP.Model, sc::UnitCommitmentScenari
return return
end end
function _add_flexiramp_reserve_eqs!(model::JuMP.Model, sc::UnitCommitmentScenario)::Nothing function _add_flexiramp_reserve_eqs!(
model::JuMP.Model,
sc::UnitCommitmentScenario,
)::Nothing
# Note: The flexpramp requirements in Wang & Hobbs (2016) are imposed as hard constraints # Note: The flexpramp requirements in Wang & Hobbs (2016) are imposed as hard constraints
# through Eq. (17) and Eq. (18). The constraints eq_min_upflexiramp and eq_min_dwflexiramp # through Eq. (17) and Eq. (18). The constraints eq_min_upflexiramp and eq_min_dwflexiramp
# provided below are modified versions of Eq. (17) and Eq. (18), respectively, in that # provided below are modified versions of Eq. (17) and Eq. (18), respectively, in that
@ -67,39 +84,37 @@ function _add_flexiramp_reserve_eqs!(model::JuMP.Model, sc::UnitCommitmentScenar
eq_min_dwflexiramp = _init(model, :eq_min_dwflexiramp) eq_min_dwflexiramp = _init(model, :eq_min_dwflexiramp)
T = model[:instance].time T = model[:instance].time
for r in sc.reserves for r in sc.reserves
if r.type == "up-frp" r.type == "flexiramp" || continue
for t in 1:T for t in 1:T
# Eq. (17) in Wang & Hobbs (2016) # Eq. (17) in Wang & Hobbs (2016)
eq_min_upflexiramp[sc.name, r.name, t] = @constraint( eq_min_upflexiramp[sc.name, r.name, t] = @constraint(
model, model,
sum(model[:upflexiramp][sc.name, r.name, g.name, t] for g in r.units) + sum(
model[:upflexiramp_shortfall][sc.name, r.name, t] >= r.amount[t] model[:upflexiramp][sc.name, r.name, g.name, t] for
) g in r.units
# Account for flexiramp shortfall contribution to objective ) + model[:upflexiramp_shortfall][sc.name, r.name, t] >=
if r.shortfall_penalty >= 0 r.amount[t]
add_to_expression!( )
model[:obj], # Eq. (18) in Wang & Hobbs (2016)
r.shortfall_penalty * sc.probability, eq_min_dwflexiramp[sc.name, r.name, t] = @constraint(
model[:upflexiramp_shortfall][sc.name, r.name, t] model,
) sum(
end model[:dwflexiramp][sc.name, r.name, g.name, t] for
end g in r.units
elseif r.type == "down-frp" ) + model[:dwflexiramp_shortfall][sc.name, r.name, t] >=
for t in 1:T r.amount[t]
# Eq. (18) in Wang & Hobbs (2016) )
eq_min_dwflexiramp[sc.name, r.name, t] = @constraint(
model, # Account for flexiramp shortfall contribution to objective
sum(model[:dwflexiramp][sc.name, r.name, g.name, t] for g in r.units) + if r.shortfall_penalty >= 0
model[:dwflexiramp_shortfall][sc.name, r.name, t] >= r.amount[t] add_to_expression!(
) model[:obj],
# Account for flexiramp shortfall contribution to objective r.shortfall_penalty * sc.probability,
if r.shortfall_penalty >= 0 (
add_to_expression!( model[:upflexiramp_shortfall][sc.name, r.name, t] +
model[:obj],
r.shortfall_penalty * sc.probability,
model[:dwflexiramp_shortfall][sc.name, r.name, t] model[:dwflexiramp_shortfall][sc.name, r.name, t]
) ),
end )
end end
end end
end end

154
src/model/formulations/base/unit.jl
Unescape View File

@ -2,7 +2,13 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details. # Released under the modified BSD license. See COPYING.md for more details.
function _add_unit_first_stage!(model::JuMP.Model, g::Unit, formulation::Formulation) # Function for adding variables, constraints, and objective function terms
# related to the binary commitment, startup and shutdown decisions of units
function _add_unit_commitment!(
model::JuMP.Model,
g::Unit,
formulation::Formulation,
)
if !all(g.must_run) && any(g.must_run) if !all(g.must_run) && any(g.must_run)
error("Partially must-run units are not currently supported") error("Partially must-run units are not currently supported")
end end
@ -21,24 +27,30 @@ function _add_unit_first_stage!(model::JuMP.Model, g::Unit, formulation::Formula
return return
end end
function _add_unit_second_stage!(model::JuMP.Model, g::Unit, formulation::Formulation, # Function for adding variables, constraints, and objective function terms
scenario::UnitCommitmentScenario) # related to the continuous dispatch decisions of units
function _add_unit_dispatch!(
model::JuMP.Model,
g::Unit,
formulation::Formulation,
sc::UnitCommitmentScenario,
)
# Variables # Variables
_add_production_vars!(model, g, formulation.prod_vars, scenario) _add_production_vars!(model, g, formulation.prod_vars, sc)
_add_spinning_reserve_vars!(model, g, scenario) _add_spinning_reserve_vars!(model, g, sc)
_add_flexiramp_reserve_vars!(model, g, scenario) _add_flexiramp_reserve_vars!(model, g, sc)
# Constraints and objective function # Constraints and objective function
_add_net_injection_eqs!(model, g, scenario) _add_net_injection_eqs!(model, g, sc)
_add_production_limit_eqs!(model, g, formulation.prod_vars, scenario) _add_production_limit_eqs!(model, g, formulation.prod_vars, sc)
_add_production_piecewise_linear_eqs!( _add_production_piecewise_linear_eqs!(
model, model,
g, g,
formulation.prod_vars, formulation.prod_vars,
formulation.pwl_costs, formulation.pwl_costs,
formulation.status_vars, formulation.status_vars,
scenario sc,
) )
_add_ramp_eqs!( _add_ramp_eqs!(
model, model,
@ -46,62 +58,29 @@ function _add_unit_second_stage!(model::JuMP.Model, g::Unit, formulation::Formul
formulation.prod_vars, formulation.prod_vars,
formulation.ramping, formulation.ramping,
formulation.status_vars, formulation.status_vars,
scenario sc,
) )
_add_startup_shutdown_limit_eqs!(model, g, scenario) _add_startup_shutdown_limit_eqs!(model, g, sc)
return return
end end
# function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
# if !all(g.must_run) && any(g.must_run)
# error("Partially must-run units are not currently supported")
# end
# if g.initial_power === nothing || g.initial_status === nothing
# error("Initial conditions for $(g.name) must be provided")
# end
# # Variables
# _add_production_vars!(model, g, formulation.prod_vars)
# _add_spinning_reserve_vars!(model, g)
# _add_flexiramp_reserve_vars!(model, g)
# _add_startup_shutdown_vars!(model, g)
# _add_status_vars!(model, g, formulation.status_vars)
# # Constraints and objective function
# _add_min_uptime_downtime_eqs!(model, g)
# _add_net_injection_eqs!(model, g)
# _add_production_limit_eqs!(model, g, formulation.prod_vars)
# _add_production_piecewise_linear_eqs!(
# model,
# g,
# formulation.prod_vars,
# formulation.pwl_costs,
# formulation.status_vars,
# )
# _add_ramp_eqs!(
# model,
# g,
# formulation.prod_vars,
# formulation.ramping,
# formulation.status_vars,
# )
# _add_startup_cost_eqs!(model, g, formulation.startup_costs)
# _add_startup_shutdown_limit_eqs!(model, g)
# _add_status_eqs!(model, g, formulation.status_vars)
# return
# end
_is_initially_on(g::Unit)::Float64 = (g.initial_status > 0 ? 1.0 : 0.0) _is_initially_on(g::Unit)::Float64 = (g.initial_status > 0 ? 1.0 : 0.0)
function _add_spinning_reserve_vars!(model::JuMP.Model, g::Unit, sc::UnitCommitmentScenario)::Nothing function _add_spinning_reserve_vars!(
model::JuMP.Model,
g::Unit,
sc::UnitCommitmentScenario,
)::Nothing
reserve = _init(model, :reserve) reserve = _init(model, :reserve)
reserve_shortfall = _init(model, :reserve_shortfall) reserve_shortfall = _init(model, :reserve_shortfall)
for r in g.reserves for r in g.reserves
r.type == "spinning" || continue r.type == "spinning" || continue
for t in 1:model[:instance].time for t in 1:model[:instance].time
reserve[sc.name, r.name, g.name, t] = @variable(model, lower_bound = 0) reserve[sc.name, r.name, g.name, t] =
@variable(model, lower_bound = 0)
if (sc.name, r.name, t) ∉ keys(reserve_shortfall) if (sc.name, r.name, t) ∉ keys(reserve_shortfall)
reserve_shortfall[sc.name, r.name, t] = @variable(model, lower_bound = 0) reserve_shortfall[sc.name, r.name, t] =
@variable(model, lower_bound = 0)
if r.shortfall_penalty < 0 if r.shortfall_penalty < 0
set_upper_bound(reserve_shortfall[sc.name, r.name, t], 0.0) set_upper_bound(reserve_shortfall[sc.name, r.name, t], 0.0)
end end
@ -111,35 +90,37 @@ function _add_spinning_reserve_vars!(model::JuMP.Model, g::Unit, sc::UnitCommitm
return return
end end
function _add_flexiramp_reserve_vars!(model::JuMP.Model, g::Unit, sc::UnitCommitmentScenario)::Nothing function _add_flexiramp_reserve_vars!(
model::JuMP.Model,
g::Unit,
sc::UnitCommitmentScenario,
)::Nothing
upflexiramp = _init(model, :upflexiramp) upflexiramp = _init(model, :upflexiramp)
upflexiramp_shortfall = _init(model, :upflexiramp_shortfall) upflexiramp_shortfall = _init(model, :upflexiramp_shortfall)
mfg = _init(model, :mfg) mfg = _init(model, :mfg)
dwflexiramp = _init(model, :dwflexiramp) dwflexiramp = _init(model, :dwflexiramp)
dwflexiramp_shortfall = _init(model, :dwflexiramp_shortfall) dwflexiramp_shortfall = _init(model, :dwflexiramp_shortfall)
for r in g.reserves for t in 1:model[:instance].time
if r.type == "up-frp" # maximum feasible generation, \bar{g_{its}} in Wang & Hobbs (2016)
for t in 1:model[:instance].time mfg[sc.name, g.name, t] = @variable(model, lower_bound = 0)
# maximum feasible generation, \bar{g_{its}} in Wang & Hobbs (2016) for r in g.reserves
mfg[sc.name, r.name, g.name, t] = @variable(model, lower_bound = 0) r.type == "flexiramp" || continue
upflexiramp[sc.name, r.name, g.name, t] = @variable(model) # up-flexiramp, ur_{it} in Wang & Hobbs (2016) upflexiramp[sc.name, r.name, g.name, t] = @variable(model) # up-flexiramp, ur_{it} in Wang & Hobbs (2016)
if (sc.name, r.name, t) ∉ keys(upflexiramp_shortfall) dwflexiramp[sc.name, r.name, g.name, t] = @variable(model) # down-flexiramp, dr_{it} in Wang & Hobbs (2016)
upflexiramp_shortfall[sc.name, r.name, t] = if (sc.name, r.name, t) ∉ keys(upflexiramp_shortfall)
@variable(model, lower_bound = 0) upflexiramp_shortfall[sc.name, r.name, t] =
if r.shortfall_penalty < 0 @variable(model, lower_bound = 0)
set_upper_bound(upflexiramp_shortfall[sc.name, r.name, t], 0.0) dwflexiramp_shortfall[sc.name, r.name, t] =
end @variable(model, lower_bound = 0)
end if r.shortfall_penalty < 0
end set_upper_bound(
elseif r.type == "down-frp" upflexiramp_shortfall[sc.name, r.name, t],
for t in 1:model[:instance].time 0.0,
dwflexiramp[sc.name, r.name, g.name, t] = @variable(model) # down-flexiramp, dr_{it} in Wang & Hobbs (2016) )
if (sc.name, r.name, t) ∉ keys(dwflexiramp_shortfall) set_upper_bound(
dwflexiramp_shortfall[sc.name, r.name, t] = dwflexiramp_shortfall[sc.name, r.name, t],
@variable(model, lower_bound = 0) 0.0,
if r.shortfall_penalty < 0 )
set_upper_bound(dwflexiramp_shortfall[sc.name, r.name, t], 0.0)
end
end end
end end
end end
@ -157,7 +138,11 @@ function _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing
return return
end end
function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit, sc::UnitCommitmentScenario)::Nothing function _add_startup_shutdown_limit_eqs!(
model::JuMP.Model,
g::Unit,
sc::UnitCommitmentScenario,
)::Nothing
eq_shutdown_limit = _init(model, :eq_shutdown_limit) eq_shutdown_limit = _init(model, :eq_shutdown_limit)
eq_startup_limit = _init(model, :eq_startup_limit) eq_startup_limit = _init(model, :eq_startup_limit)
is_on = model[:is_on] is_on = model[:is_on]
@ -196,7 +181,7 @@ function _add_ramp_eqs!(
model::JuMP.Model, model::JuMP.Model,
g::Unit, g::Unit,
formulation::RampingFormulation, formulation::RampingFormulation,
sc::UnitCommitmentScenario sc::UnitCommitmentScenario,
)::Nothing )::Nothing
prod_above = model[:prod_above] prod_above = model[:prod_above]
reserve = _total_reserves(model, g, sc) reserve = _total_reserves(model, g, sc)
@ -282,7 +267,11 @@ function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
end end
end end
function _add_net_injection_eqs!(model::JuMP.Model, g::Unit, sc::UnitCommitmentScenario)::Nothing function _add_net_injection_eqs!(
model::JuMP.Model,
g::Unit,
sc::UnitCommitmentScenario,
)::Nothing
expr_net_injection = model[:expr_net_injection] expr_net_injection = model[:expr_net_injection]
for t in 1:model[:instance].time for t in 1:model[:instance].time
# Add to net injection expression # Add to net injection expression
@ -305,7 +294,10 @@ function _total_reserves(model, g, sc)::Vector
spinning_reserves = [r for r in g.reserves if r.type == "spinning"] spinning_reserves = [r for r in g.reserves if r.type == "spinning"]
if !isempty(spinning_reserves) if !isempty(spinning_reserves)
reserve += [ reserve += [
sum(model[:reserve][sc.name, r.name, g.name, t] for r in spinning_reserves) for t in 1:model[:instance].time sum(
model[:reserve][sc.name, r.name, g.name, t] for
r in spinning_reserves
) for t in 1:model[:instance].time
] ]
end end
return reserve return reserve

46
src/solution/fix.jl
Unescape View File

@ -10,37 +10,43 @@ solution. Useful for computing LMPs.
""" """
function fix!(model::JuMP.Model, solution::AbstractDict)::Nothing function fix!(model::JuMP.Model, solution::AbstractDict)::Nothing
instance, T = model[:instance], model[:instance].time instance, T = model[:instance], model[:instance].time
"Production (MW)" ∈ keys(solution) ? solution = Dict("s1" => solution) :
nothing
is_on = model[:is_on] is_on = model[:is_on]
prod_above = model[:prod_above] prod_above = model[:prod_above]
reserve = model[:reserve] reserve = model[:reserve]
for g in instance.units for sc in instance.scenarios
for t in 1:T for g in sc.units
is_on_value = round(solution["Is on"][g.name][t])
prod_value =
round(solution["Production (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],
prod_value - is_on_value * g.min_power[t],
force = true,
)
end
end
for r in instance.reserves
r.type == "spinning" || continue
for g in r.units
for t in 1:T for t in 1:T
reserve_value = round( is_on_value = round(solution[sc.name]["Is on"][g.name][t])
solution["Spinning reserve (MW)"][r.name][g.name][t], prod_value = round(
solution[sc.name]["Production (MW)"][g.name][t],
digits = 5, digits = 5,
) )
JuMP.fix(is_on[g.name, t], is_on_value, force = true)
JuMP.fix( JuMP.fix(
reserve[r.name, g.name, t], prod_above[sc.name, g.name, t],
reserve_value, prod_value - is_on_value * g.min_power[t],
force = true, force = true,
) )
end end
end end
for r in sc.reserves
r.type == "spinning" || continue
for g in r.units
for t in 1:T
reserve_value = round(
solution[sc.name]["Spinning reserve (MW)"][r.name][g.name][t],
digits = 5,
)
JuMP.fix(
reserve[sc.name, r.name, g.name, t],
reserve_value,
force = true,
)
end
end
end
end end
return return
end end

2
src/solution/methods/XavQiuWanThi2019/enforce.jl
Unescape View File

@ -5,7 +5,7 @@
function _enforce_transmission( function _enforce_transmission(
model::JuMP.Model, model::JuMP.Model,
violations::Vector{_Violation}, violations::Vector{_Violation},
sc::UnitCommitmentScenario sc::UnitCommitmentScenario,
)::Nothing )::Nothing
for v in violations for v in violations
_enforce_transmission( _enforce_transmission(

16
src/solution/methods/XavQiuWanThi2019/find.jl
Unescape View File

@ -19,8 +19,8 @@ function _find_violations(
time_screening = @elapsed begin time_screening = @elapsed begin
non_slack_buses = [b for b in sc.buses if b.offset > 0] non_slack_buses = [b for b in sc.buses if b.offset > 0]
net_injection_values = [ net_injection_values = [
value(net_injection[sc.name, b.name, t]) for b in non_slack_buses, value(net_injection[sc.name, b.name, t]) for
t in 1:instance.time b in non_slack_buses, t in 1:instance.time
] ]
overflow_values = [ overflow_values = [
value(overflow[sc.name, lm.name, t]) for lm in sc.lines, value(overflow[sc.name, lm.name, t]) for lm in sc.lines,
@ -72,7 +72,7 @@ function _find_violations(;
isf::Array{Float64,2}, isf::Array{Float64,2},
lodf::Array{Float64,2}, lodf::Array{Float64,2},
max_per_line::Int, max_per_line::Int,
max_per_period::Int max_per_period::Int,
)::Array{_Violation,1} )::Array{_Violation,1}
B = length(sc.buses) - 1 B = length(sc.buses) - 1
L = length(sc.lines) L = length(sc.lines)
@ -96,13 +96,13 @@ function _find_violations(;
post_v::Array{Float64} = zeros(L, L, K) # post_v[lm, lc, thread] post_v::Array{Float64} = zeros(L, L, K) # post_v[lm, lc, thread]
normal_limits::Array{Float64,2} = [ normal_limits::Array{Float64,2} = [
l.normal_flow_limit[t] + overflow[l.offset, t] for l.normal_flow_limit[t] + overflow[l.offset, t] for l in sc.lines,
l in sc.lines, t in 1:T t in 1:T
] ]
emergency_limits::Array{Float64,2} = [ emergency_limits::Array{Float64,2} = [
l.emergency_flow_limit[t] + overflow[l.offset, t] for l.emergency_flow_limit[t] + overflow[l.offset, t] for l in sc.lines,
l in sc.lines, t in 1:T t in 1:T
] ]
is_vulnerable::Array{Bool} = zeros(Bool, L) is_vulnerable::Array{Bool} = zeros(Bool, L)
@ -114,7 +114,7 @@ function _find_violations(;
k = threadid() k = threadid()
# Pre-contingency flows # Pre-contingency flows
pre_flow[:, k] = isf * net_injections[ :, t] pre_flow[:, k] = isf * net_injections[:, t]
# Post-contingency flows # Post-contingency flows
for lc in 1:L, lm in 1:L for lc in 1:L, lm in 1:L

8
src/solution/methods/XavQiuWanThi2019/optimize.jl
Unescape View File

@ -10,9 +10,9 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing
JuMP.set_optimizer_attribute(model, "MIPGap", gap) JuMP.set_optimizer_attribute(model, "MIPGap", gap)
@info @sprintf("MIP gap tolerance set to %f", gap) @info @sprintf("MIP gap tolerance set to %f", gap)
end end
for scenario in model[:instance].scenarios for sc in model[:instance].scenarios
large_gap = false large_gap = false
has_transmission = (length(scenario.isf) > 0) has_transmission = (length(sc.isf) > 0)
if has_transmission && method.two_phase_gap if has_transmission && method.two_phase_gap
set_gap(1e-2) set_gap(1e-2)
large_gap = true large_gap = true
@ -36,7 +36,7 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing
has_transmission || break has_transmission || break
violations = _find_violations( violations = _find_violations(
model, model,
scenario, sc,
max_per_line = method.max_violations_per_line, max_per_line = method.max_violations_per_line,
max_per_period = method.max_violations_per_period, max_per_period = method.max_violations_per_period,
) )
@ -49,7 +49,7 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing
break break
end end
else else
_enforce_transmission(model, violations, scenario) _enforce_transmission(model, violations, sc)
end end
end end
end end

72
src/solution/solution.jl
Unescape View File

@ -16,17 +16,21 @@ solution = UnitCommitment.solution(model)
""" """
function solution(model::JuMP.Model)::OrderedDict function solution(model::JuMP.Model)::OrderedDict
instance, T = model[:instance], model[:instance].time instance, T = model[:instance], model[:instance].time
function timeseries_first_stage(vars, collection) function timeseries(vars, collection; sc = nothing)
return OrderedDict( if sc === nothing
b.name => [round(value(vars[b.name, t]), digits = 5) for t in 1:T] return OrderedDict(
for b in collection b.name =>
) [round(value(vars[b.name, t]), digits = 5) for t in 1:T] for
end b in collection
function timeseries_second_stage(vars, collection, sc) )
return OrderedDict( else
b.name => [round(value(vars[sc.name, b.name, t]), digits = 5) for t in 1:T] return OrderedDict(
for b in collection b.name => [
) round(value(vars[sc.name, b.name, t]), digits = 5) for
t in 1:T
] for b in collection
)
end
end end
function production_cost(g, sc) function production_cost(g, sc)
return [ return [
@ -67,24 +71,25 @@ function solution(model::JuMP.Model)::OrderedDict
OrderedDict(g.name => production_cost(g, sc) for g in sc.units) OrderedDict(g.name => production_cost(g, sc) for g in sc.units)
sol[sc.name]["Startup cost (\$)"] = sol[sc.name]["Startup cost (\$)"] =
OrderedDict(g.name => startup_cost(g, sc) for g in sc.units) OrderedDict(g.name => startup_cost(g, sc) for g in sc.units)
sol[sc.name]["Is on"] = timeseries_first_stage(model[:is_on], sc.units) sol[sc.name]["Is on"] = timeseries(model[:is_on], sc.units)
sol[sc.name]["Switch on"] = timeseries_first_stage(model[:switch_on], sc.units) sol[sc.name]["Switch on"] = timeseries(model[:switch_on], sc.units)
sol[sc.name]["Switch off"] = timeseries_first_stage(model[:switch_off], sc.units) sol[sc.name]["Switch off"] = timeseries(model[:switch_off], sc.units)
sol[sc.name]["Net injection (MW)"] = sol[sc.name]["Net injection (MW)"] =
timeseries_second_stage(model[:net_injection], sc.buses, sc) timeseries(model[:net_injection], sc.buses, sc = sc)
sol[sc.name]["Load curtail (MW)"] = timeseries_second_stage(model[:curtail], sc.buses, sc) sol[sc.name]["Load curtail (MW)"] =
timeseries(model[:curtail], sc.buses, sc = sc)
if !isempty(sc.lines) if !isempty(sc.lines)
sol[sc.name]["Line overflow (MW)"] = timeseries_second_stage(model[:overflow], sc.lines, sc) sol[sc.name]["Line overflow (MW)"] =
timeseries(model[:overflow], sc.lines, sc = sc)
end end
if !isempty(sc.price_sensitive_loads) if !isempty(sc.price_sensitive_loads)
sol[sc.name]["Price-sensitive loads (MW)"] = sol[sc.name]["Price-sensitive loads (MW)"] =
timeseries_second_stage(model[:loads], sc.price_sensitive_loads, sc) timeseries(model[:loads], sc.price_sensitive_loads, sc = sc)
end end
sol[sc.name]["Spinning reserve (MW)"] = OrderedDict( sol[sc.name]["Spinning reserve (MW)"] = OrderedDict(
r.name => OrderedDict( r.name => OrderedDict(
g.name => [ g.name => [
value(model[:reserve][sc.name, r.name, g.name, t]) for value(model[:reserve][sc.name, r.name, g.name, t]) for t in 1:instance.time
t in 1:instance.time
] for g in r.units ] for g in r.units
) for r in sc.reserves if r.type == "spinning" ) for r in sc.reserves if r.type == "spinning"
) )
@ -97,32 +102,31 @@ function solution(model::JuMP.Model)::OrderedDict
sol[sc.name]["Up-flexiramp (MW)"] = OrderedDict( sol[sc.name]["Up-flexiramp (MW)"] = OrderedDict(
r.name => OrderedDict( r.name => OrderedDict(
g.name => [ g.name => [
value(model[:upflexiramp][sc.name, r.name, g.name, t]) for value(model[:upflexiramp][sc.name, r.name, g.name, t]) for t in 1:instance.time
t in 1:instance.time
] for g in r.units ] for g in r.units
) for r in sc.reserves if r.type == "up-frp" ) for r in sc.reserves if r.type == "flexiramp"
) )
sol[sc.name]["Up-flexiramp shortfall (MW)"] = OrderedDict( sol[sc.name]["Up-flexiramp shortfall (MW)"] = OrderedDict(
r.name => [ r.name => [
value(model[:upflexiramp_shortfall][sc.name, r.name, t]) for value(model[:upflexiramp_shortfall][sc.name, r.name, t]) for t in 1:instance.time
t in 1:instance.time ] for r in sc.reserves if r.type == "flexiramp"
] for r in sc.reserves if r.type == "up-frp"
) )
sol[sc.name]["Down-flexiramp (MW)"] = OrderedDict( sol[sc.name]["Down-flexiramp (MW)"] = OrderedDict(
r.name => OrderedDict( r.name => OrderedDict(
g.name => [ g.name => [
value(model[:dwflexiramp][sc.name, r.name, g.name, t]) for value(model[:dwflexiramp][sc.name, r.name, g.name, t]) for t in 1:instance.time
t in 1:instance.time
] for g in r.units ] for g in r.units
) for r in sc.reserves if r.type == "down-frp" ) for r in sc.reserves if r.type == "flexiramp"
) )
sol[sc.name]["Down-flexiramp shortfall (MW)"] = OrderedDict( sol[sc.name]["Down-flexiramp shortfall (MW)"] = OrderedDict(
r.name => [ r.name => [
value(model[:dwflexiramp_shortfall][sc.name, r.name, t]) for value(model[:dwflexiramp_shortfall][sc.name, r.name, t]) for t in 1:instance.time
t in 1:instance.time ] for r in sc.reserves if r.type == "flexiramp"
] for r in sc.reserves if r.type == "down-frp"
) )
end end
length(keys(sol)) > 1 ? nothing : sol = Dict(sol_key => sol_val for scen_key in keys(sol) for (sol_key, sol_val) in sol[scen_key]) if length(instance.scenarios) == 1
return sol return first(values(sol))
else
return sol
end
end end

12
src/transform/initcond.jl
Unescape View File

@ -5,18 +5,18 @@
using JuMP using JuMP
""" """
generate_initial_conditions!(instance, optimizer) generate_initial_conditions!(sc, optimizer)
Generates feasible initial conditions for the given instance, by constructing Generates feasible initial conditions for the given scenario, by constructing
and solving a single-period mixed-integer optimization problem, using the given and solving a single-period mixed-integer optimization problem, using the given
optimizer. The instance is modified in-place. optimizer. The scenario is modified in-place.
""" """
function generate_initial_conditions!( function generate_initial_conditions!(
instance::UnitCommitmentInstance, sc::UnitCommitmentScenario,
optimizer, optimizer,
)::Nothing )::Nothing
G = instance.units G = sc.units
B = instance.buses B = sc.buses
t = 1 t = 1
mip = JuMP.Model(optimizer) mip = JuMP.Model(optimizer)

50
src/transform/randomize/XavQiuAhm2021.jl
Unescape View File

@ -6,6 +6,7 @@ module XavQiuAhm2021
using Distributions using Distributions
import ..UnitCommitmentInstance import ..UnitCommitmentInstance
import ..UnitCommitmentScenario
""" """
struct Randomization struct Randomization
@ -119,10 +120,10 @@ end
function _randomize_costs( function _randomize_costs(
rng, rng,
instance::UnitCommitmentInstance, sc::UnitCommitmentScenario,
distribution, distribution,
)::Nothing )::Nothing
for unit in instance.units for unit in sc.units
α = rand(rng, distribution) α = rand(rng, distribution)
unit.min_power_cost *= α unit.min_power_cost *= α
for k in unit.cost_segments for k in unit.cost_segments
@ -137,17 +138,15 @@ end
function _randomize_load_share( function _randomize_load_share(
rng, rng,
instance::UnitCommitmentInstance, sc::UnitCommitmentScenario,
distribution, distribution,
)::Nothing )::Nothing
α = rand(rng, distribution, length(instance.buses)) α = rand(rng, distribution, length(sc.buses))
for t in 1:instance.time for t in 1:sc.time
total = sum(bus.load[t] for bus in instance.buses) total = sum(bus.load[t] for bus in sc.buses)
den = sum( den =
bus.load[t] / total * α[i] for sum(bus.load[t] / total * α[i] for (i, bus) in enumerate(sc.buses))
(i, bus) in enumerate(instance.buses) for (i, bus) in enumerate(sc.buses)
)
for (i, bus) in enumerate(instance.buses)
bus.load[t] *= α[i] / den bus.load[t] *= α[i] / den
end end
end end
@ -156,12 +155,12 @@ end
function _randomize_load_profile( function _randomize_load_profile(
rng, rng,
instance::UnitCommitmentInstance, sc::UnitCommitmentScenario,
params::Randomization, params::Randomization,
)::Nothing )::Nothing
# Generate new system load # Generate new system load
system_load = [1.0] system_load = [1.0]
for t in 2:instance.time for t in 2:sc.time
idx = (t - 1) % length(params.load_profile_mu) + 1 idx = (t - 1) % length(params.load_profile_mu) + 1
gamma = rand( gamma = rand(
rng, rng,
@ -169,14 +168,14 @@ function _randomize_load_profile(
) )
push!(system_load, system_load[t-1] * gamma) push!(system_load, system_load[t-1] * gamma)
end end
capacity = sum(maximum(u.max_power) for u in instance.units) capacity = sum(maximum(u.max_power) for u in sc.units)
peak_load = rand(rng, params.peak_load) * capacity peak_load = rand(rng, params.peak_load) * capacity
system_load = system_load ./ maximum(system_load) .* peak_load system_load = system_load ./ maximum(system_load) .* peak_load
# Scale bus loads to match the new system load # Scale bus loads to match the new system load
prev_system_load = sum(b.load for b in instance.buses) prev_system_load = sum(b.load for b in sc.buses)
for b in instance.buses for b in sc.buses
for t in 1:instance.time for t in 1:sc.time
b.load[t] *= system_load[t] / prev_system_load[t] b.load[t] *= system_load[t] / prev_system_load[t]
end end
end end
@ -199,15 +198,26 @@ function randomize!(
instance::UnitCommitment.UnitCommitmentInstance, instance::UnitCommitment.UnitCommitmentInstance,
method::XavQiuAhm2021.Randomization; method::XavQiuAhm2021.Randomization;
rng = MersenneTwister(), rng = MersenneTwister(),
)::Nothing
for sc in instance.scenarios
randomize!(sc; method, rng)
end
return
end
function randomize!(
sc::UnitCommitment.UnitCommitmentScenario;
method::XavQiuAhm2021.Randomization,
rng = MersenneTwister(),
)::Nothing )::Nothing
if method.randomize_costs if method.randomize_costs
XavQiuAhm2021._randomize_costs(rng, instance, method.cost) XavQiuAhm2021._randomize_costs(rng, sc, method.cost)
end end
if method.randomize_load_share if method.randomize_load_share
XavQiuAhm2021._randomize_load_share(rng, instance, method.load_share) XavQiuAhm2021._randomize_load_share(rng, sc, method.load_share)
end end
if method.randomize_load_profile if method.randomize_load_profile
XavQiuAhm2021._randomize_load_profile(rng, instance, method) XavQiuAhm2021._randomize_load_profile(rng, sc, method)
end end
return return
end end

50
src/transform/slice.jl
Unescape View File

@ -24,31 +24,33 @@ function slice(
)::UnitCommitmentInstance )::UnitCommitmentInstance
modified = deepcopy(instance) modified = deepcopy(instance)
modified.time = length(range) modified.time = length(range)
modified.power_balance_penalty = modified.power_balance_penalty[range] for sc in modified.scenarios
for r in modified.reserves sc.power_balance_penalty = sc.power_balance_penalty[range]
r.amount = r.amount[range] for r in sc.reserves
end r.amount = r.amount[range]
for u in modified.units end
u.max_power = u.max_power[range] for u in sc.units
u.min_power = u.min_power[range] u.max_power = u.max_power[range]
u.must_run = u.must_run[range] u.min_power = u.min_power[range]
u.min_power_cost = u.min_power_cost[range] u.must_run = u.must_run[range]
for s in u.cost_segments u.min_power_cost = u.min_power_cost[range]
s.mw = s.mw[range] for s in u.cost_segments
s.cost = s.cost[range] s.mw = s.mw[range]
s.cost = s.cost[range]
end
end
for b in sc.buses
b.load = b.load[range]
end
for l in sc.lines
l.normal_flow_limit = l.normal_flow_limit[range]
l.emergency_flow_limit = l.emergency_flow_limit[range]
l.flow_limit_penalty = l.flow_limit_penalty[range]
end
for ps in sc.price_sensitive_loads
ps.demand = ps.demand[range]
ps.revenue = ps.revenue[range]
end end
end
for b in modified.buses
b.load = b.load[range]
end
for l in modified.lines
l.normal_flow_limit = l.normal_flow_limit[range]
l.emergency_flow_limit = l.emergency_flow_limit[range]
l.flow_limit_penalty = l.flow_limit_penalty[range]
end
for ps in modified.price_sensitive_loads
ps.demand = ps.demand[range]
ps.revenue = ps.revenue[range]
end end
return modified return modified
end end

611
src/validation/validate.jl
Unescape View File

@ -28,6 +28,8 @@ function validate(
instance::UnitCommitmentInstance, instance::UnitCommitmentInstance,
solution::Union{Dict,OrderedDict}, solution::Union{Dict,OrderedDict},
)::Bool )::Bool
"Production (MW)" ∈ keys(solution) ? solution = Dict("s1" => solution) :
nothing
err_count = 0 err_count = 0
err_count += _validate_units(instance, solution) err_count += _validate_units(instance, solution)
err_count += _validate_reserve_and_demand(instance, solution) err_count += _validate_reserve_and_demand(instance, solution)
@ -42,358 +44,369 @@ end
function _validate_units(instance::UnitCommitmentInstance, solution; tol = 0.01) function _validate_units(instance::UnitCommitmentInstance, solution; tol = 0.01)
err_count = 0 err_count = 0
for sc in instance.scenarios
for unit in instance.units for unit in sc.units
production = solution["Production (MW)"][unit.name] production = solution[sc.name]["Production (MW)"][unit.name]
reserve = [0.0 for _ in 1:instance.time] reserve = [0.0 for _ in 1:instance.time]
spinning_reserves = [r for r in unit.reserves if r.type == "spinning"] spinning_reserves =
if !isempty(spinning_reserves) [r for r in unit.reserves if r.type == "spinning"]
reserve += sum( if !isempty(spinning_reserves)
solution["Spinning reserve (MW)"][r.name][unit.name] for reserve += sum(
r in spinning_reserves solution[sc.name]["Spinning reserve (MW)"][r.name][unit.name]
) for r in spinning_reserves
end )
actual_production_cost = solution["Production cost (\$)"][unit.name]
actual_startup_cost = solution["Startup cost (\$)"][unit.name]
is_on = bin(solution["Is on"][unit.name])
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_down = max(0, unit.initial_power - production[t])
else
is_starting_up = !is_on[t-1] && is_on[t]
is_shutting_down = is_on[t-1] && !is_on[t]
ramp_up = max(0, production[t] + reserve[t] - production[t-1])
ramp_down = max(0, production[t-1] - production[t])
end end
actual_production_cost =
solution[sc.name]["Production cost (\$)"][unit.name]
actual_startup_cost =
solution[sc.name]["Startup cost (\$)"][unit.name]
is_on = bin(solution[sc.name]["Is on"][unit.name])
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_down = max(0, unit.initial_power - production[t])
else
is_starting_up = !is_on[t-1] && is_on[t]
is_shutting_down = is_on[t-1] && !is_on[t]
ramp_up =
max(0, production[t] + reserve[t] - production[t-1])
ramp_down = max(0, production[t-1] - production[t])
end
# Compute production costs # Compute production costs
production_cost, startup_cost = 0, 0 production_cost, startup_cost = 0, 0
if is_on[t] if is_on[t]
production_cost += unit.min_power_cost[t] production_cost += unit.min_power_cost[t]
residual = max(0, production[t] - unit.min_power[t]) residual = max(0, production[t] - unit.min_power[t])
for s in unit.cost_segments for s in unit.cost_segments
cleared = min(residual, s.mw[t]) cleared = min(residual, s.mw[t])
production_cost += cleared * s.cost[t] production_cost += cleared * s.cost[t]
residual = max(0, residual - s.mw[t]) residual = max(0, residual - s.mw[t])
end
end end
end
# Production should be non-negative # Production should be non-negative
if production[t] < -tol if production[t] < -tol
@error @sprintf( @error @sprintf(
"Unit %s produces negative amount of power at time %d (%.2f)", "Unit %s produces negative amount of power at time %d (%.2f)",
unit.name, unit.name,
t, t,
production[t] production[t]
) )
err_count += 1 err_count += 1
end end
# Verify must-run # Verify must-run
if !is_on[t] && unit.must_run[t] if !is_on[t] && unit.must_run[t]
@error @sprintf( @error @sprintf(
"Must-run unit %s is offline at time %d", "Must-run unit %s is offline at time %d",
unit.name, unit.name,
t t
) )
err_count += 1 err_count += 1
end end
# Verify reserve eligibility # Verify reserve eligibility
for r in instance.reserves for r in sc.reserves
if r.type == "spinning" if r.type == "spinning"
if unit ∉ r.units && if unit ∉ r.units && (
(unit in keys(solution["Spinning reserve (MW)"][r.name])) unit in keys(
@error @sprintf( solution[sc.name]["Spinning reserve (MW)"][r.name],
"Unit %s is not eligible to provide reserve %s", )
unit.name,
r.name,
) )
err_count += 1 @error @sprintf(
"Unit %s is not eligible to provide reserve %s",
unit.name,
r.name,
)
err_count += 1
end
end end
end end
end
# If unit is on, must produce at least its minimum power
if is_on[t] && (production[t] < unit.min_power[t] - tol)
@error @sprintf(
"Unit %s produces below its minimum limit at time %d (%.2f < %.2f)",
unit.name,
t,
production[t],
unit.min_power[t]
)
err_count += 1
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)
@error @sprintf(
"Unit %s produces above its maximum limit at time %d (%.2f + %.2f> %.2f)",
unit.name,
t,
production[t],
reserve[t],
unit.max_power[t]
)
err_count += 1
end
# If unit is off, must produce zero # If unit is on, must produce at least its minimum power
if !is_on[t] && production[t] + reserve[t] > tol if is_on[t] && (production[t] < unit.min_power[t] - tol)
@error @sprintf( @error @sprintf(
"Unit %s produces power at time %d while off (%.2f + %.2f > 0)", "Unit %s produces below its minimum limit at time %d (%.2f < %.2f)",
unit.name, unit.name,
t, t,
production[t], production[t],
reserve[t], unit.min_power[t]
) )
err_count += 1 err_count += 1
end end
# Startup limit
if is_starting_up && (ramp_up > unit.startup_limit + tol)
@error @sprintf(
"Unit %s exceeds startup limit at time %d (%.2f > %.2f)",
unit.name,
t,
ramp_up,
unit.startup_limit
)
err_count += 1
end
# Shutdown limit
if is_shutting_down && (ramp_down > unit.shutdown_limit + tol)
@error @sprintf(
"Unit %s exceeds shutdown limit at time %d (%.2f > %.2f)",
unit.name,
t,
ramp_down,
unit.shutdown_limit
)
err_count += 1
end
# Ramp-up limit
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,
t,
ramp_up,
unit.ramp_up_limit
)
err_count += 1
end
# Ramp-down limit # If unit is on, must produce at most its maximum power
if !is_starting_up && if is_on[t] &&
!is_shutting_down && (production[t] + reserve[t] > unit.max_power[t] + tol)
(ramp_down > unit.ramp_down_limit + tol) @error @sprintf(
@error @sprintf( "Unit %s produces above its maximum limit at time %d (%.2f + %.2f> %.2f)",
"Unit %s exceeds ramp down limit at time %d (%.2f > %.2f)", unit.name,
unit.name, t,
t, production[t],
ramp_down, reserve[t],
unit.ramp_down_limit unit.max_power[t]
) )
err_count += 1 err_count += 1
end end
# Verify startup costs & minimum downtime # If unit is off, must produce zero
if is_starting_up if !is_on[t] && production[t] + reserve[t] > tol
@error @sprintf(
"Unit %s produces power at time %d while off (%.2f + %.2f > 0)",
unit.name,
t,
production[t],
reserve[t],
)
err_count += 1
end
# Calculate how much time the unit has been offline # Startup limit
time_down = 0 if is_starting_up && (ramp_up > unit.startup_limit + tol)
for k in 1:(t-1) @error @sprintf(
if !is_on[t-k] "Unit %s exceeds startup limit at time %d (%.2f > %.2f)",
time_down += 1 unit.name,
else t,
break ramp_up,
end unit.startup_limit
)
err_count += 1
end end
if (t == time_down + 1) && (unit.initial_status < 0)
time_down -= unit.initial_status # Shutdown limit
if is_shutting_down && (ramp_down > unit.shutdown_limit + tol)
@error @sprintf(
"Unit %s exceeds shutdown limit at time %d (%.2f > %.2f)",
unit.name,
t,
ramp_down,
unit.shutdown_limit
)
err_count += 1
end end
# Calculate startup costs # Ramp-up limit
for c in unit.startup_categories if !is_starting_up &&
if time_down >= c.delay !is_shutting_down &&
startup_cost = c.cost (ramp_up > unit.ramp_up_limit + tol)
end @error @sprintf(
"Unit %s exceeds ramp up limit at time %d (%.2f > %.2f)",
unit.name,
t,
ramp_up,
unit.ramp_up_limit
)
err_count += 1
end end
# Check minimum downtime # Ramp-down limit
if time_down < unit.min_downtime if !is_starting_up &&
!is_shutting_down &&
(ramp_down > unit.ramp_down_limit + tol)
@error @sprintf( @error @sprintf(
"Unit %s violates minimum downtime at time %d", "Unit %s exceeds ramp down limit at time %d (%.2f > %.2f)",
unit.name, unit.name,
t t,
ramp_down,
unit.ramp_down_limit
) )
err_count += 1 err_count += 1
end end
end
# Verify minimum uptime # Verify startup costs & minimum downtime
if is_shutting_down if is_starting_up
# Calculate how much time the unit has been offline
time_down = 0
for k in 1:(t-1)
if !is_on[t-k]
time_down += 1
else
break
end
end
if (t == time_down + 1) && (unit.initial_status < 0)
time_down -= unit.initial_status
end
# Calculate startup costs
for c in unit.startup_categories
if time_down >= c.delay
startup_cost = c.cost
end
end
# Calculate how much time the unit has been online # Check minimum downtime
time_up = 0 if time_down < unit.min_downtime
for k in 1:(t-1) @error @sprintf(
if is_on[t-k] "Unit %s violates minimum downtime at time %d",
time_up += 1 unit.name,
else t
break )
err_count += 1
end end
end end
if (t == time_up + 1) && (unit.initial_status > 0)
time_up += unit.initial_status # Verify minimum uptime
if is_shutting_down
# Calculate how much time the unit has been online
time_up = 0
for k in 1:(t-1)
if is_on[t-k]
time_up += 1
else
break
end
end
if (t == time_up + 1) && (unit.initial_status > 0)
time_up += unit.initial_status
end
# Check minimum uptime
if time_up < unit.min_uptime
@error @sprintf(
"Unit %s violates minimum uptime at time %d",
unit.name,
t
)
err_count += 1
end
end end
# Check minimum uptime # Verify production costs
if time_up < unit.min_uptime if abs(actual_production_cost[t] - production_cost) > 1.00
@error @sprintf( @error @sprintf(
"Unit %s violates minimum uptime at time %d", "Unit %s has unexpected production cost at time %d (%.2f should be %.2f)",
unit.name, unit.name,
t t,
actual_production_cost[t],
production_cost
) )
err_count += 1 err_count += 1
end end
end
# Verify production costs
if abs(actual_production_cost[t] - production_cost) > 1.00
@error @sprintf(
"Unit %s has unexpected production cost at time %d (%.2f should be %.2f)",
unit.name,
t,
actual_production_cost[t],
production_cost
)
err_count += 1
end
# Verify startup costs # Verify startup costs
if abs(actual_startup_cost[t] - startup_cost) > 1.00 if abs(actual_startup_cost[t] - startup_cost) > 1.00
@error @sprintf( @error @sprintf(
"Unit %s has unexpected startup cost at time %d (%.2f should be %.2f)", "Unit %s has unexpected startup cost at time %d (%.2f should be %.2f)",
unit.name, unit.name,
t, t,
actual_startup_cost[t], actual_startup_cost[t],
startup_cost startup_cost
) )
err_count += 1 err_count += 1
end
end end
end end
end end
return err_count return err_count
end end
function _validate_reserve_and_demand(instance, solution, tol = 0.01) function _validate_reserve_and_demand(instance, solution, tol = 0.01)
err_count = 0 err_count = 0
for t in 1:instance.time for sc in instance.scenarios
load_curtail = 0 for t in 1:instance.time
fixed_load = sum(b.load[t] for b in instance.buses) load_curtail = 0
ps_load = 0 fixed_load = sum(b.load[t] for b in sc.buses)
if length(instance.price_sensitive_loads) > 0 ps_load = 0
ps_load = sum( if length(sc.price_sensitive_loads) > 0
solution["Price-sensitive loads (MW)"][ps.name][t] for ps_load = sum(
ps in instance.price_sensitive_loads solution[sc.name]["Price-sensitive loads (MW)"][ps.name][t]
) for ps in sc.price_sensitive_loads
end )
production = end
sum(solution["Production (MW)"][g.name][t] for g in instance.units) production = sum(
if "Load curtail (MW)" in keys(solution) solution[sc.name]["Production (MW)"][g.name][t] for
load_curtail = sum( g in sc.units
solution["Load curtail (MW)"][b.name][t] for
b in instance.buses
)
end
balance = fixed_load - load_curtail - production + ps_load
# Verify that production equals demand
if abs(balance) > tol
@error @sprintf(
"Non-zero power balance at time %d (%.2f + %.2f - %.2f - %.2f != 0)",
t,
fixed_load,
ps_load,
load_curtail,
production,
) )
err_count += 1 if "Load curtail (MW)" in keys(solution)
end load_curtail = sum(
solution[sc.name]["Load curtail (MW)"][b.name][t] for
b in sc.buses
)
end
balance = fixed_load - load_curtail - production + ps_load
# Verify reserves # Verify that production equals demand
for r in instance.reserves if abs(balance) > tol
if r.type == "spinning" @error @sprintf(
provided = sum( "Non-zero power balance at time %d (%.2f + %.2f - %.2f - %.2f != 0)",
solution["Spinning reserve (MW)"][r.name][g.name][t] for t,
g in r.units fixed_load,
ps_load,
load_curtail,
production,
) )
shortfall = err_count += 1
solution["Spinning reserve shortfall (MW)"][r.name][t] end
required = r.amount[t]
if provided + shortfall < required - tol # Verify reserves
@error @sprintf( for r in sc.reserves
"Insufficient reserve %s at time %d (%.2f + %.2f < %.2f)", if r.type == "spinning"
r.name, provided = sum(
t, solution[sc.name]["Spinning reserve (MW)"][r.name][g.name][t]
provided, for g in r.units
shortfall,
required,
) )
end shortfall =
elseif r.type == "up-frp" solution[sc.name]["Spinning reserve shortfall (MW)"][r.name][t]
upflexiramp = sum( required = r.amount[t]
solution["Up-flexiramp (MW)"][r.name][g.name][t] for
g in r.units
)
upflexiramp_shortfall =
solution["Up-flexiramp shortfall (MW)"][r.name][t]
if upflexiramp + upflexiramp_shortfall < r.amount[t] - tol if provided + shortfall < required - tol
@error @sprintf( @error @sprintf(
"Insufficient up-flexiramp at time %d (%.2f + %.2f < %.2f)", "Insufficient reserve %s at time %d (%.2f + %.2f < %.2f)",
t, r.name,
upflexiramp, t,
upflexiramp_shortfall, provided,
r.amount[t], shortfall,
required,
)
end
elseif r.type == "flexiramp"
upflexiramp = sum(
solution[sc.name]["Up-flexiramp (MW)"][r.name][g.name][t]
for g in r.units
) )
err_count += 1 upflexiramp_shortfall =
end solution[sc.name]["Up-flexiramp shortfall (MW)"][r.name][t]
elseif r.type == "down-frp"
dwflexiramp = sum(
solution["Down-flexiramp (MW)"][r.name][g.name][t] for
g in r.units
)
dwflexiramp_shortfall =
solution["Down-flexiramp shortfall (MW)"][r.name][t]
if dwflexiramp + dwflexiramp_shortfall < r.amount[t] - tol if upflexiramp + upflexiramp_shortfall < r.amount[t] - tol
@error @sprintf( @error @sprintf(
"Insufficient down-flexiramp at time %d (%.2f + %.2f < %.2f)", "Insufficient up-flexiramp at time %d (%.2f + %.2f < %.2f)",
t, t,
dwflexiramp, upflexiramp,
dwflexiramp_shortfall, upflexiramp_shortfall,
r.amount[t], r.amount[t],
)
err_count += 1
end
dwflexiramp = sum(
solution[sc.name]["Down-flexiramp (MW)"][r.name][g.name][t]
for g in r.units
) )
err_count += 1 dwflexiramp_shortfall =
solution[sc.name]["Down-flexiramp shortfall (MW)"][r.name][t]
if dwflexiramp + dwflexiramp_shortfall < r.amount[t] - tol
@error @sprintf(
"Insufficient down-flexiramp at time %d (%.2f + %.2f < %.2f)",
t,
dwflexiramp,
dwflexiramp_shortfall,
r.amount[t],
)
err_count += 1
end
else
error("Unknown reserve type: $(r.type)")
end end
else
error("Unknown reserve type: $(r.type)")
end end
end end
end end

1
test/Project.toml
Unescape View File

@ -2,6 +2,7 @@
Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76" Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63" GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63"
JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572" JuMP = "4076af6c-e467-56ae-b986-b466b2749572"

14
test/instance/migrate_test.jl
Unescape View File

@ -6,13 +6,17 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
@testset "read v0.2" begin @testset "read v0.2" begin
instance = UnitCommitment.read("$FIXTURES/ucjl-0.2.json.gz") instance = UnitCommitment.read("$FIXTURES/ucjl-0.2.json.gz")
@test length(instance.reserves_by_name["r1"].amount) == 4 for sc in instance.scenarios
@test instance.units_by_name["g2"].reserves[1].name == "r1" @test length(sc.reserves_by_name["r1"].amount) == 4
@test sc.units_by_name["g2"].reserves[1].name == "r1"
end
end end
@testset "read v0.3" begin @testset "read v0.3" begin
instance = UnitCommitment.read("$FIXTURES/ucjl-0.3.json.gz") instance = UnitCommitment.read("$FIXTURES/ucjl-0.3.json.gz")
@test length(instance.units) == 6 for sc in instance.scenarios
@test length(instance.buses) == 14 @test length(sc.units) == 6
@test length(instance.lines) == 20 @test length(sc.buses) == 14
@test length(sc.lines) == 20
end
end end

191
test/instance/read_test.jl
Unescape View File

@ -6,115 +6,116 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
@testset "read_benchmark" begin @testset "read_benchmark" begin
instance = UnitCommitment.read("$FIXTURES/case14.json.gz") instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
for sc in instance.scenarios
@test length(sc.lines) == 20
@test length(sc.buses) == 14
@test length(sc.units) == 6
@test length(sc.contingencies) == 19
@test length(sc.price_sensitive_loads) == 1
@test instance.time == 4
@test length(instance.lines) == 20 @test sc.lines[5].name == "l5"
@test length(instance.buses) == 14 @test sc.lines[5].source.name == "b2"
@test length(instance.units) == 6 @test sc.lines[5].target.name == "b5"
@test length(instance.contingencies) == 19 @test sc.lines[5].reactance ≈ 0.17388
@test length(instance.price_sensitive_loads) == 1 @test sc.lines[5].susceptance ≈ 10.037550333
@test instance.time == 4 @test sc.lines[5].normal_flow_limit == [1e8 for t in 1:4]
@test sc.lines[5].emergency_flow_limit == [1e8 for t in 1:4]
@test sc.lines[5].flow_limit_penalty == [5e3 for t in 1:4]
@test sc.lines_by_name["l5"].name == "l5"
@test instance.lines[5].name == "l5" @test sc.lines[1].name == "l1"
@test instance.lines[5].source.name == "b2" @test sc.lines[1].source.name == "b1"
@test instance.lines[5].target.name == "b5" @test sc.lines[1].target.name == "b2"
@test instance.lines[5].reactance ≈ 0.17388 @test sc.lines[1].reactance ≈ 0.059170
@test instance.lines[5].susceptance ≈ 10.037550333 @test sc.lines[1].susceptance ≈ 29.496860773945
@test instance.lines[5].normal_flow_limit == [1e8 for t in 1:4] @test sc.lines[1].normal_flow_limit == [300.0 for t in 1:4]
@test instance.lines[5].emergency_flow_limit == [1e8 for t in 1:4] @test sc.lines[1].emergency_flow_limit == [400.0 for t in 1:4]
@test instance.lines[5].flow_limit_penalty == [5e3 for t in 1:4] @test sc.lines[1].flow_limit_penalty == [1e3 for t in 1:4]
@test instance.lines_by_name["l5"].name == "l5"
@test instance.lines[1].name == "l1" @test sc.buses[9].name == "b9"
@test instance.lines[1].source.name == "b1" @test sc.buses[9].load == [35.36638, 33.25495, 31.67138, 31.14353]
@test instance.lines[1].target.name == "b2" @test sc.buses_by_name["b9"].name == "b9"
@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.buses[9].name == "b9" @test sc.reserves[1].name == "r1"
@test instance.buses[9].load == [35.36638, 33.25495, 31.67138, 31.14353] @test sc.reserves[1].type == "spinning"
@test instance.buses_by_name["b9"].name == "b9" @test sc.reserves[1].amount == [100.0, 100.0, 100.0, 100.0]
@test sc.reserves_by_name["r1"].name == "r1"
@test instance.reserves[1].name == "r1" unit = sc.units[1]
@test instance.reserves[1].type == "spinning" @test unit.name == "g1"
@test instance.reserves[1].amount == [100.0, 100.0, 100.0, 100.0] @test unit.bus.name == "b1"
@test instance.reserves_by_name["r1"].name == "r1" @test unit.ramp_up_limit == 1e6
@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.min_uptime == 1
@test unit.min_downtime == 1
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
@test unit.cost_segments[1].cost[t] ≈ 20.0
@test unit.cost_segments[2].cost[t] ≈ 30.0
@test unit.cost_segments[3].cost[t] ≈ 40.0
end
@test length(unit.startup_categories) == 3
@test unit.startup_categories[1].delay == 1
@test unit.startup_categories[2].delay == 2
@test unit.startup_categories[3].delay == 3
@test unit.startup_categories[1].cost == 1000.0
@test unit.startup_categories[2].cost == 1500.0
@test unit.startup_categories[3].cost == 2000.0
@test length(unit.reserves) == 0
@test sc.units_by_name["g1"].name == "g1"
unit = instance.units[1] unit = sc.units[2]
@test unit.name == "g1" @test unit.name == "g2"
@test unit.bus.name == "b1" @test unit.must_run == [false for t in 1:4]
@test unit.ramp_up_limit == 1e6 @test length(unit.reserves) == 1
@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.min_uptime == 1
@test unit.min_downtime == 1
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
@test unit.cost_segments[1].cost[t] ≈ 20.0
@test unit.cost_segments[2].cost[t] ≈ 30.0
@test unit.cost_segments[3].cost[t] ≈ 40.0
end
@test length(unit.startup_categories) == 3
@test unit.startup_categories[1].delay == 1
@test unit.startup_categories[2].delay == 2
@test unit.startup_categories[3].delay == 3
@test unit.startup_categories[1].cost == 1000.0
@test unit.startup_categories[2].cost == 1500.0
@test unit.startup_categories[3].cost == 2000.0
@test length(unit.reserves) == 0
@test instance.units_by_name["g1"].name == "g1"
unit = instance.units[2] unit = sc.units[3]
@test unit.name == "g2" @test unit.name == "g3"
@test unit.must_run == [false for t in 1:4] @test unit.bus.name == "b3"
@test length(unit.reserves) == 1 @test unit.ramp_up_limit == 70.0
@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.min_uptime == 1
@test unit.min_downtime == 1
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
@test unit.cost_segments[1].cost[t] ≈ 33.75
@test unit.cost_segments[2].cost[t] ≈ 38.04
@test unit.cost_segments[3].cost[t] ≈ 44.77853
end
@test length(unit.reserves) == 1
@test unit.reserves[1].name == "r1"
unit = instance.units[3] @test sc.contingencies[1].lines == [sc.lines[1]]
@test unit.name == "g3" @test sc.contingencies[1].units == []
@test unit.bus.name == "b3" @test sc.contingencies[1].name == "c1"
@test unit.ramp_up_limit == 70.0 @test sc.contingencies_by_name["c1"].name == "c1"
@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.min_uptime == 1
@test unit.min_downtime == 1
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
@test unit.cost_segments[1].cost[t] ≈ 33.75
@test unit.cost_segments[2].cost[t] ≈ 38.04
@test unit.cost_segments[3].cost[t] ≈ 44.77853
end
@test length(unit.reserves) == 1
@test unit.reserves[1].name == "r1"
@test instance.contingencies[1].lines == [instance.lines[1]]
@test instance.contingencies[1].units == []
@test instance.contingencies[1].name == "c1"
@test instance.contingencies_by_name["c1"].name == "c1"
load = instance.price_sensitive_loads[1] load = sc.price_sensitive_loads[1]
@test load.name == "ps1" @test load.name == "ps1"
@test load.bus.name == "b3" @test load.bus.name == "b3"
@test load.revenue == [100.0 for t in 1:4] @test load.revenue == [100.0 for t in 1:4]
@test load.demand == [50.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" @test sc.price_sensitive_loads_by_name["ps1"].name == "ps1"
end
end end
@testset "read_benchmark sub-hourly" begin @testset "read_benchmark sub-hourly" begin
instance = UnitCommitment.read("$FIXTURES/case14-sub-hourly.json.gz") instance = UnitCommitment.read("$FIXTURES/case14-sub-hourly.json.gz")
@test instance.time == 4 @test instance.time == 4
unit = instance.units[1] unit = instance.scenarios[1].units[1]
@test unit.name == "g1" @test unit.name == "g1"
@test unit.min_uptime == 2 @test unit.min_uptime == 2
@test unit.min_downtime == 2 @test unit.min_downtime == 2

143
test/solution/methods/XavQiuWanThi19/filter_test.jl
Unescape View File

@ -8,76 +8,77 @@ import UnitCommitment: _Violation, _offer, _query
@testset "_ViolationFilter" begin @testset "_ViolationFilter" begin
instance = UnitCommitment.read("$FIXTURES/case14.json.gz") instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
filter = UnitCommitment._ViolationFilter(max_per_line = 1, max_total = 2) filter = UnitCommitment._ViolationFilter(max_per_line = 1, max_total = 2)
for sc in instance.scenarios
_offer(
filter,
_Violation(
time = 1,
monitored_line = sc.lines[1],
outage_line = nothing,
amount = 100.0,
),
)
_offer(
filter,
_Violation(
time = 1,
monitored_line = sc.lines[1],
outage_line = sc.lines[1],
amount = 300.0,
),
)
_offer(
filter,
_Violation(
time = 1,
monitored_line = sc.lines[1],
outage_line = sc.lines[5],
amount = 500.0,
),
)
_offer(
filter,
_Violation(
time = 1,
monitored_line = sc.lines[1],
outage_line = sc.lines[4],
amount = 400.0,
),
)
_offer(
filter,
_Violation(
time = 1,
monitored_line = sc.lines[2],
outage_line = sc.lines[1],
amount = 200.0,
),
)
_offer(
filter,
_Violation(
time = 1,
monitored_line = sc.lines[2],
outage_line = sc.lines[8],
amount = 100.0,
),
)
_offer( actual = _query(filter)
filter, expected = [
_Violation( _Violation(
time = 1, time = 1,
monitored_line = instance.lines[1], monitored_line = sc.lines[2],
outage_line = nothing, outage_line = sc.lines[1],
amount = 100.0, amount = 200.0,
), ),
) _Violation(
_offer( time = 1,
filter, monitored_line = sc.lines[1],
_Violation( outage_line = sc.lines[5],
time = 1, amount = 500.0,
monitored_line = instance.lines[1], ),
outage_line = instance.lines[1], ]
amount = 300.0, @test actual == expected
), end
)
_offer(
filter,
_Violation(
time = 1,
monitored_line = instance.lines[1],
outage_line = instance.lines[5],
amount = 500.0,
),
)
_offer(
filter,
_Violation(
time = 1,
monitored_line = instance.lines[1],
outage_line = instance.lines[4],
amount = 400.0,
),
)
_offer(
filter,
_Violation(
time = 1,
monitored_line = instance.lines[2],
outage_line = instance.lines[1],
amount = 200.0,
),
)
_offer(
filter,
_Violation(
time = 1,
monitored_line = instance.lines[2],
outage_line = instance.lines[8],
amount = 100.0,
),
)
actual = _query(filter)
expected = [
_Violation(
time = 1,
monitored_line = instance.lines[2],
outage_line = instance.lines[1],
amount = 200.0,
),
_Violation(
time = 1,
monitored_line = instance.lines[1],
outage_line = instance.lines[5],
amount = 500.0,
),
]
@test actual == expected
end end

51
test/solution/methods/XavQiuWanThi19/find_test.jl
Unescape View File

@ -7,29 +7,32 @@ import UnitCommitment: _Violation, _offer, _query
@testset "find_violations" begin @testset "find_violations" begin
instance = UnitCommitment.read("$FIXTURES/case14.json.gz") instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
for line in instance.lines, t in 1:instance.time for sc in instance.scenarios
line.normal_flow_limit[t] = 1.0 for line in sc.lines, t in 1:instance.time
line.emergency_flow_limit[t] = 1.0 line.normal_flow_limit[t] = 1.0
line.emergency_flow_limit[t] = 1.0
end
isf = UnitCommitment._injection_shift_factors(
lines = sc.lines,
buses = sc.buses,
)
lodf = UnitCommitment._line_outage_factors(
lines = sc.lines,
buses = sc.buses,
isf = isf,
)
inj = [1000.0 for b in 1:13, t in 1:instance.time]
overflow = [0.0 for l in sc.lines, t in 1:instance.time]
violations = UnitCommitment._find_violations(
instance = instance,
sc = sc,
net_injections = inj,
overflow = overflow,
isf = isf,
lodf = lodf,
max_per_line = 1,
max_per_period = 5,
)
@test length(violations) == 20
end end
isf = UnitCommitment._injection_shift_factors(
lines = instance.lines,
buses = instance.buses,
)
lodf = UnitCommitment._line_outage_factors(
lines = instance.lines,
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]
violations = UnitCommitment._find_violations(
instance = instance,
net_injections = inj,
overflow = overflow,
isf = isf,
lodf = lodf,
max_per_line = 1,
max_per_period = 5,
)
@test length(violations) == 20
end end

246
test/solution/methods/XavQiuWanThi19/sensitivity_test.jl
Unescape View File

@ -6,137 +6,145 @@ using UnitCommitment, Test, LinearAlgebra
@testset "_susceptance_matrix" begin @testset "_susceptance_matrix" begin
instance = UnitCommitment.read("$FIXTURES/case14.json.gz") instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
actual = UnitCommitment._susceptance_matrix(instance.lines) for sc in instance.scenarios
@test size(actual) == (20, 20) actual = UnitCommitment._susceptance_matrix(sc.lines)
expected = Diagonal([ @test size(actual) == (20, 20)
29.5, expected = Diagonal([
7.83, 29.5,
8.82, 7.83,
9.9, 8.82,
10.04, 9.9,
10.2, 10.04,
41.45, 10.2,
8.35, 41.45,
3.14, 8.35,
6.93, 3.14,
8.77, 6.93,
6.82, 8.77,
13.4, 6.82,
9.91, 13.4,
15.87, 9.91,
20.65, 15.87,
6.46, 20.65,
9.09, 6.46,
8.73, 9.09,
5.02, 8.73,
]) 5.02,
@test round.(actual, digits = 2) == expected ])
@test round.(actual, digits = 2) == expected
end
end end
@testset "_reduced_incidence_matrix" begin @testset "_reduced_incidence_matrix" begin
instance = UnitCommitment.read("$FIXTURES/case14.json.gz") instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
actual = UnitCommitment._reduced_incidence_matrix( for sc in instance.scenarios
lines = instance.lines, actual = UnitCommitment._reduced_incidence_matrix(
buses = instance.buses, lines = sc.lines,
) buses = sc.buses,
@test size(actual) == (20, 13) )
@test actual[1, 1] == -1.0 @test size(actual) == (20, 13)
@test actual[3, 1] == 1.0 @test actual[1, 1] == -1.0
@test actual[4, 1] == 1.0 @test actual[3, 1] == 1.0
@test actual[5, 1] == 1.0 @test actual[4, 1] == 1.0
@test actual[3, 2] == -1.0 @test actual[5, 1] == 1.0
@test actual[6, 2] == 1.0 @test actual[3, 2] == -1.0
@test actual[4, 3] == -1.0 @test actual[6, 2] == 1.0
@test actual[6, 3] == -1.0 @test actual[4, 3] == -1.0
@test actual[7, 3] == 1.0 @test actual[6, 3] == -1.0
@test actual[8, 3] == 1.0 @test actual[7, 3] == 1.0
@test actual[9, 3] == 1.0 @test actual[8, 3] == 1.0
@test actual[2, 4] == -1.0 @test actual[9, 3] == 1.0
@test actual[5, 4] == -1.0 @test actual[2, 4] == -1.0
@test actual[7, 4] == -1.0 @test actual[5, 4] == -1.0
@test actual[10, 4] == 1.0 @test actual[7, 4] == -1.0
@test actual[10, 5] == -1.0 @test actual[10, 4] == 1.0
@test actual[11, 5] == 1.0 @test actual[10, 5] == -1.0
@test actual[12, 5] == 1.0 @test actual[11, 5] == 1.0
@test actual[13, 5] == 1.0 @test actual[12, 5] == 1.0
@test actual[8, 6] == -1.0 @test actual[13, 5] == 1.0
@test actual[14, 6] == 1.0 @test actual[8, 6] == -1.0
@test actual[15, 6] == 1.0 @test actual[14, 6] == 1.0
@test actual[14, 7] == -1.0 @test actual[15, 6] == 1.0
@test actual[9, 8] == -1.0 @test actual[14, 7] == -1.0
@test actual[15, 8] == -1.0 @test actual[9, 8] == -1.0
@test actual[16, 8] == 1.0 @test actual[15, 8] == -1.0
@test actual[17, 8] == 1.0 @test actual[16, 8] == 1.0
@test actual[16, 9] == -1.0 @test actual[17, 8] == 1.0
@test actual[18, 9] == 1.0 @test actual[16, 9] == -1.0
@test actual[11, 10] == -1.0 @test actual[18, 9] == 1.0
@test actual[18, 10] == -1.0 @test actual[11, 10] == -1.0
@test actual[12, 11] == -1.0 @test actual[18, 10] == -1.0
@test actual[19, 11] == 1.0 @test actual[12, 11] == -1.0
@test actual[13, 12] == -1.0 @test actual[19, 11] == 1.0
@test actual[19, 12] == -1.0 @test actual[13, 12] == -1.0
@test actual[20, 12] == 1.0 @test actual[19, 12] == -1.0
@test actual[17, 13] == -1.0 @test actual[20, 12] == 1.0
@test actual[20, 13] == -1.0 @test actual[17, 13] == -1.0
@test actual[20, 13] == -1.0
end
end end
@testset "_injection_shift_factors" begin @testset "_injection_shift_factors" begin
instance = UnitCommitment.read("$FIXTURES/case14.json.gz") instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
actual = UnitCommitment._injection_shift_factors( for sc in instance.scenarios
lines = instance.lines, actual = UnitCommitment._injection_shift_factors(
buses = instance.buses, lines = sc.lines,
) buses = sc.buses,
@test size(actual) == (20, 13) )
@test round.(actual, digits = 2) == [ @test size(actual) == (20, 13)
-0.84 -0.75 -0.67 -0.61 -0.63 -0.66 -0.66 -0.65 -0.65 -0.64 -0.63 -0.63 -0.64 @test round.(actual, digits = 2) == [
-0.16 -0.25 -0.33 -0.39 -0.37 -0.34 -0.34 -0.35 -0.35 -0.36 -0.37 -0.37 -0.36 -0.84 -0.75 -0.67 -0.61 -0.63 -0.66 -0.66 -0.65 -0.65 -0.64 -0.63 -0.63 -0.64
0.03 -0.53 -0.15 -0.1 -0.12 -0.14 -0.14 -0.14 -0.13 -0.13 -0.12 -0.12 -0.13 -0.16 -0.25 -0.33 -0.39 -0.37 -0.34 -0.34 -0.35 -0.35 -0.36 -0.37 -0.37 -0.36
0.06 -0.14 -0.32 -0.22 -0.25 -0.3 -0.3 -0.29 -0.28 -0.27 -0.25 -0.26 -0.27 0.03 -0.53 -0.15 -0.1 -0.12 -0.14 -0.14 -0.14 -0.13 -0.13 -0.12 -0.12 -0.13
0.08 -0.07 -0.2 -0.29 -0.26 -0.22 -0.22 -0.22 -0.23 -0.25 -0.26 -0.26 -0.24 0.06 -0.14 -0.32 -0.22 -0.25 -0.3 -0.3 -0.29 -0.28 -0.27 -0.25 -0.26 -0.27
0.03 0.47 -0.15 -0.1 -0.12 -0.14 -0.14 -0.14 -0.13 -0.13 -0.12 -0.12 -0.13 0.08 -0.07 -0.2 -0.29 -0.26 -0.22 -0.22 -0.22 -0.23 -0.25 -0.26 -0.26 -0.24
0.08 0.31 0.5 -0.3 -0.03 0.36 0.36 0.28 0.23 0.1 -0.0 0.02 0.17 0.03 0.47 -0.15 -0.1 -0.12 -0.14 -0.14 -0.14 -0.13 -0.13 -0.12 -0.12 -0.13
0.0 0.01 0.02 -0.01 -0.22 -0.63 -0.63 -0.45 -0.41 -0.32 -0.24 -0.25 -0.36 0.08 0.31 0.5 -0.3 -0.03 0.36 0.36 0.28 0.23 0.1 -0.0 0.02 0.17
0.0 0.01 0.01 -0.01 -0.12 -0.17 -0.17 -0.26 -0.24 -0.18 -0.14 -0.14 -0.21 0.0 0.01 0.02 -0.01 -0.22 -0.63 -0.63 -0.45 -0.41 -0.32 -0.24 -0.25 -0.36
-0.0 -0.02 -0.03 0.02 -0.66 -0.2 -0.2 -0.29 -0.36 -0.5 -0.63 -0.61 -0.43 0.0 0.01 0.01 -0.01 -0.12 -0.17 -0.17 -0.26 -0.24 -0.18 -0.14 -0.14 -0.21
-0.0 -0.01 -0.02 0.01 0.21 -0.12 -0.12 -0.17 -0.28 -0.53 0.18 0.15 -0.03 -0.0 -0.02 -0.03 0.02 -0.66 -0.2 -0.2 -0.29 -0.36 -0.5 -0.63 -0.61 -0.43
-0.0 -0.0 -0.0 0.0 0.03 -0.02 -0.02 -0.03 -0.02 0.01 -0.52 -0.17 -0.09 -0.0 -0.01 -0.02 0.01 0.21 -0.12 -0.12 -0.17 -0.28 -0.53 0.18 0.15 -0.03
-0.0 -0.01 -0.01 0.01 0.11 -0.06 -0.06 -0.09 -0.05 0.02 -0.28 -0.59 -0.31 -0.0 -0.0 -0.0 0.0 0.03 -0.02 -0.02 -0.03 -0.02 0.01 -0.52 -0.17 -0.09
-0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 -0.0 -0.01 -0.01 0.01 0.11 -0.06 -0.06 -0.09 -0.05 0.02 -0.28 -0.59 -0.31
0.0 0.01 0.02 -0.01 -0.22 0.37 0.37 -0.45 -0.41 -0.32 -0.24 -0.25 -0.36 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0
0.0 0.01 0.02 -0.01 -0.21 0.12 0.12 0.17 -0.72 -0.47 -0.18 -0.15 0.03 0.0 0.01 0.02 -0.01 -0.22 0.37 0.37 -0.45 -0.41 -0.32 -0.24 -0.25 -0.36
0.0 0.01 0.01 -0.01 -0.14 0.08 0.08 0.12 0.07 -0.03 -0.2 -0.24 -0.6 0.0 0.01 0.02 -0.01 -0.21 0.12 0.12 0.17 -0.72 -0.47 -0.18 -0.15 0.03
0.0 0.01 0.02 -0.01 -0.21 0.12 0.12 0.17 0.28 -0.47 -0.18 -0.15 0.03 0.0 0.01 0.01 -0.01 -0.14 0.08 0.08 0.12 0.07 -0.03 -0.2 -0.24 -0.6
-0.0 -0.0 -0.0 0.0 0.03 -0.02 -0.02 -0.03 -0.02 0.01 0.48 -0.17 -0.09 0.0 0.01 0.02 -0.01 -0.21 0.12 0.12 0.17 0.28 -0.47 -0.18 -0.15 0.03
-0.0 -0.01 -0.01 0.01 0.14 -0.08 -0.08 -0.12 -0.07 0.03 0.2 0.24 -0.4 -0.0 -0.0 -0.0 0.0 0.03 -0.02 -0.02 -0.03 -0.02 0.01 0.48 -0.17 -0.09
] -0.0 -0.01 -0.01 0.01 0.14 -0.08 -0.08 -0.12 -0.07 0.03 0.2 0.24 -0.4
]
end
end end
@testset "_line_outage_factors" begin @testset "_line_outage_factors" begin
instance = UnitCommitment.read("$FIXTURES/case14.json.gz") instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
isf_before = UnitCommitment._injection_shift_factors( for sc in instance.scenarios
lines = instance.lines, isf_before = UnitCommitment._injection_shift_factors(
buses = instance.buses, lines = sc.lines,
) buses = sc.buses,
lodf = UnitCommitment._line_outage_factors( )
lines = instance.lines, lodf = UnitCommitment._line_outage_factors(
buses = instance.buses, lines = sc.lines,
isf = isf_before, buses = sc.buses,
) isf = isf_before,
for contingency in instance.contingencies )
for lc in contingency.lines for contingency in sc.contingencies
prev_susceptance = lc.susceptance for lc in contingency.lines
lc.susceptance = 0.0 prev_susceptance = lc.susceptance
isf_after = UnitCommitment._injection_shift_factors( lc.susceptance = 0.0
lines = instance.lines, isf_after = UnitCommitment._injection_shift_factors(
buses = instance.buses, lines = sc.lines,
) buses = sc.buses,
lc.susceptance = prev_susceptance )
for lm in instance.lines lc.susceptance = prev_susceptance
expected = isf_after[lm.offset, :] for lm in sc.lines
actual = expected = isf_after[lm.offset, :]
isf_before[lm.offset, :] + actual =
lodf[lm.offset, lc.offset] * isf_before[lc.offset, :] isf_before[lm.offset, :] +
@test norm(expected - actual) < 1e-6 lodf[lm.offset, lc.offset] * isf_before[lc.offset, :]
@test norm(expected - actual) < 1e-6
end
end end
end end
end end

25
test/transform/initcond_test.jl
Unescape View File

@ -8,20 +8,21 @@ using UnitCommitment, Cbc, JuMP
# Load instance # Load instance
instance = UnitCommitment.read("$FIXTURES/case118-initcond.json.gz") instance = UnitCommitment.read("$FIXTURES/case118-initcond.json.gz")
optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0) optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
for sc in instance.scenarios
# All units should have unknown initial conditions
for g in sc.units
@test g.initial_power === nothing
@test g.initial_status === nothing
end
# All units should have unknown initial conditions # Generate initial conditions
for g in instance.units UnitCommitment.generate_initial_conditions!(sc, optimizer)
@test g.initial_power === nothing
@test g.initial_status === nothing
end
# Generate initial conditions
UnitCommitment.generate_initial_conditions!(instance, optimizer)
# All units should now have known initial conditions # All units should now have known initial conditions
for g in instance.units for g in sc.units
@test g.initial_power !== nothing @test g.initial_power !== nothing
@test g.initial_status !== nothing @test g.initial_status !== nothing
end
end end
# TODO: Check that initial conditions are feasible # TODO: Check that initial conditions are feasible

96
test/transform/randomize/XavQiuAhm2021_test.jl
Unescape View File

@ -10,60 +10,78 @@ using Random
using UnitCommitment, Cbc, JuMP using UnitCommitment, Cbc, JuMP
get_instance() = UnitCommitment.read_benchmark("matpower/case118/2017-02-01") get_instance() = UnitCommitment.read_benchmark("matpower/case118/2017-02-01")
system_load(instance) = sum(b.load for b in instance.buses) system_load(sc) = sum(b.load for b in sc.buses)
test_approx(x, y) = @test isapprox(x, y, atol = 1e-3) test_approx(x, y) = @test isapprox(x, y, atol = 1e-3)
@testset "XavQiuAhm2021" begin @testset "XavQiuAhm2021" begin
@testset "cost and load share" begin @testset "cost and load share" begin
instance = get_instance() instance = get_instance()
for sc in instance.scenarios
# Check original costs
unit = sc.units[10]
test_approx(unit.min_power_cost[1], 825.023)
test_approx(unit.cost_segments[1].cost[1], 36.659)
test_approx(unit.startup_categories[1].cost[1], 7570.42)
# Check original costs # Check original load share
unit = instance.units[10] bus = sc.buses[1]
test_approx(unit.min_power_cost[1], 825.023) prev_system_load = system_load(sc)
test_approx(unit.cost_segments[1].cost[1], 36.659) test_approx(bus.load[1] / prev_system_load[1], 0.012)
test_approx(unit.startup_categories[1].cost[1], 7570.42)
# Check original load share randomize!(
bus = instance.buses[1] sc,
prev_system_load = system_load(instance) method = XavQiuAhm2021.Randomization(
test_approx(bus.load[1] / prev_system_load[1], 0.012) randomize_load_profile = false,
),
rng = MersenneTwister(42),
)
randomize!( # Check randomized costs
instance, test_approx(unit.min_power_cost[1], 831.977)
method = XavQiuAhm2021.Randomization( test_approx(unit.cost_segments[1].cost[1], 36.968)
randomize_load_profile = false, test_approx(unit.startup_categories[1].cost[1], 7634.226)
),
rng = MersenneTwister(42),
)
# Check randomized costs # Check randomized load share
test_approx(unit.min_power_cost[1], 831.977) curr_system_load = system_load(sc)
test_approx(unit.cost_segments[1].cost[1], 36.968) test_approx(bus.load[1] / curr_system_load[1], 0.013)
test_approx(unit.startup_categories[1].cost[1], 7634.226)
# Check randomized load share # System load should not change
curr_system_load = system_load(instance) @test prev_system_load ≈ curr_system_load
test_approx(bus.load[1] / curr_system_load[1], 0.013) end
# System load should not change
@test prev_system_load ≈ curr_system_load
end end
@testset "load profile" begin @testset "load profile" begin
instance = get_instance() instance = get_instance()
for sc in instance.scenarios
# Check original load profile
@test round.(system_load(sc), digits = 1)[1:8] ≈ [
3059.5,
2983.2,
2937.5,
2953.9,
3073.1,
3356.4,
4068.5,
4018.8,
]
# Check original load profile randomize!(
@test round.(system_load(instance), digits = 1)[1:8] ≈ sc;
[3059.5, 2983.2, 2937.5, 2953.9, 3073.1, 3356.4, 4068.5, 4018.8] method = XavQiuAhm2021.Randomization(),
rng = MersenneTwister(42),
randomize!( )
instance,
XavQiuAhm2021.Randomization(),
rng = MersenneTwister(42),
)
# Check randomized load profile # Check randomized load profile
@test round.(system_load(instance), digits = 1)[1:8] ≈ @test round.(system_load(sc), digits = 1)[1:8] ≈ [
[4854.7, 4849.2, 4732.7, 4848.2, 4948.4, 5231.1, 5874.8, 5934.8] 4854.7,
4849.2,
4732.7,
4848.2,
4948.4,
5231.1,
5874.8,
5934.8,
]
end
end end
end end

46
test/transform/slice_test.jl
Unescape View File

@ -10,29 +10,31 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
# Should update all time-dependent fields # Should update all time-dependent fields
@test modified.time == 2 @test modified.time == 2
@test length(modified.power_balance_penalty) == 2 for sc in modified.scenarios
@test length(modified.reserves_by_name["r1"].amount) == 2 @test length(sc.power_balance_penalty) == 2
for u in modified.units @test length(sc.reserves_by_name["r1"].amount) == 2
@test length(u.max_power) == 2 for u in sc.units
@test length(u.min_power) == 2 @test length(u.max_power) == 2
@test length(u.must_run) == 2 @test length(u.min_power) == 2
@test length(u.min_power_cost) == 2 @test length(u.must_run) == 2
for s in u.cost_segments @test length(u.min_power_cost) == 2
@test length(s.mw) == 2 for s in u.cost_segments
@test length(s.cost) == 2 @test length(s.mw) == 2
@test length(s.cost) == 2
end
end
for b in sc.buses
@test length(b.load) == 2
end
for l in sc.lines
@test length(l.normal_flow_limit) == 2
@test length(l.emergency_flow_limit) == 2
@test length(l.flow_limit_penalty) == 2
end
for ps in sc.price_sensitive_loads
@test length(ps.demand) == 2
@test length(ps.revenue) == 2
end end
end
for b in modified.buses
@test length(b.load) == 2
end
for l in modified.lines
@test length(l.normal_flow_limit) == 2
@test length(l.emergency_flow_limit) == 2
@test length(l.flow_limit_penalty) == 2
end
for ps in modified.price_sensitive_loads
@test length(ps.demand) == 2
@test length(ps.revenue) == 2
end end
# Should be able to build model without errors # Should be able to build model without errors
optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0) optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)

9
test/usage.jl
Unescape View File

@ -6,8 +6,13 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON
@testset "usage" begin @testset "usage" begin
instance = UnitCommitment.read("$FIXTURES/case14.json.gz") instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
for line in instance.lines, t in 1:4 for sc in instance.scenarios
line.normal_flow_limit[t] = 10.0 sc.power_balance_penalty = [100000 for _ in 1:instance.time]
end
for sc in instance.scenarios
for line in sc.lines, t in 1:4
line.normal_flow_limit[t] = 10.0
end
end end
optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0) optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
model = UnitCommitment.build_model( model = UnitCommitment.build_model(

14
test/validation/repair_test.jl
Unescape View File

@ -16,8 +16,10 @@ end
json = parse_case14() json = parse_case14()
json["Generators"]["g1"]["Production cost curve (MW)"] = [100, 150, 200] json["Generators"]["g1"]["Production cost curve (MW)"] = [100, 150, 200]
json["Generators"]["g1"]["Production cost curve (\$)"] = [10, 25, 30] json["Generators"]["g1"]["Production cost curve (\$)"] = [10, 25, 30]
instance = UnitCommitment._from_json(json, repair = false) sc = UnitCommitment._from_json(json, repair = false)
@test UnitCommitment.repair!(instance) == 4 sc.name = "s1"
sc.probability = 1.0
@test UnitCommitment.repair!(sc) == 4
end end
@testset "Startup limit must be greater than Pmin" begin @testset "Startup limit must be greater than Pmin" begin
@ -25,15 +27,15 @@ end
json["Generators"]["g1"]["Production cost curve (MW)"] = [100, 150] json["Generators"]["g1"]["Production cost curve (MW)"] = [100, 150]
json["Generators"]["g1"]["Production cost curve (\$)"] = [100, 150] json["Generators"]["g1"]["Production cost curve (\$)"] = [100, 150]
json["Generators"]["g1"]["Startup limit (MW)"] = 80 json["Generators"]["g1"]["Startup limit (MW)"] = 80
instance = UnitCommitment._from_json(json, repair = false) sc = UnitCommitment._from_json(json, repair = false)
@test UnitCommitment.repair!(instance) == 1 @test UnitCommitment.repair!(sc) == 1
end end
@testset "Startup costs and delays must be increasing" begin @testset "Startup costs and delays must be increasing" begin
json = parse_case14() json = parse_case14()
json["Generators"]["g1"]["Startup costs (\$)"] = [300, 200, 100] json["Generators"]["g1"]["Startup costs (\$)"] = [300, 200, 100]
json["Generators"]["g1"]["Startup delays (h)"] = [8, 4, 2] json["Generators"]["g1"]["Startup delays (h)"] = [8, 4, 2]
instance = UnitCommitment._from_json(json, repair = false) sc = UnitCommitment._from_json(json, repair = false)
@test UnitCommitment.repair!(instance) == 4 @test UnitCommitment.repair!(sc) == 4
end end
end end

Write Preview
Loading…
Cancel
Save