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

Merge pull request #25 from oyurdakul/stochastic-extension

Browse Source 
stochastic extension w/ scenarios
pull/29/head
Alinson S. Xavier 3 years ago committed by GitHub
parent 8fc84412eb 414128cc0b
commit 52f1ff9a27
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
40 changed files with 1218 additions and 904 deletions
Show Stats Download Patch File Download Diff File

2
src/UnitCommitment.jl
Unescape View File

@ -4,6 +4,8 @@
module UnitCommitment
using Base: String
include("instance/structs.jl")
include("model/formulations/base/structs.jl")
include("solution/structs.jl")

90
src/instance/read.jl
Unescape View File

@ -43,26 +43,76 @@ function read_benchmark(
return UnitCommitment.read(filename)
end
function _repair_scenario_names_and_probabilities!(
scenarios::Vector{UnitCommitmentScenario},
path::Vector{String},
)::Nothing
total_weight = sum([sc.probability for sc in scenarios])
for (sc_path, sc) in zip(path, scenarios)
sc.name !== "" ||
(sc.name = first(split(last(split(sc_path, "/")), ".")))
sc.probability = (sc.probability / total_weight)
end
return
end
"""
read(path::AbstractString)::UnitCommitmentInstance
Read instance from a file. The file may be gzipped.
Read a deterministic test case from the given file. The file may be gzipped.
# Example
```julia
instance = UnitCommitment.read("s1.json.gz")
```
"""
function read(path::String)::UnitCommitmentInstance
scenarios = Vector{UnitCommitmentScenario}()
scenario = _read_scenario(path)
scenario.name = "s1"
scenario.probability = 1.0
scenarios = [scenario]
instance =
UnitCommitmentInstance(time = scenario.time, scenarios = scenarios)
return instance
end
"""
read(path::Vector{String})::UnitCommitmentInstance
Read a stochastic unit commitment instance from the given files. Each file
describes a scenario. The files may be gzipped.
# Example
```julia
instance = UnitCommitment.read("/path/to/input.json.gz")
instance = UnitCommitment.read(["s1.json.gz", "s2.json.gz"])
```
"""
function read(path::AbstractString)::UnitCommitmentInstance
function read(paths::Vector{String})::UnitCommitmentInstance
scenarios = UnitCommitmentScenario[]
for p in paths
push!(scenarios, _read_scenario(p))
end
_repair_scenario_names_and_probabilities!(scenarios, paths)
instance =
UnitCommitmentInstance(time = scenarios[1].time, scenarios = scenarios)
return instance
end
function _read_scenario(path::String)::UnitCommitmentScenario
if endswith(path, ".gz")
return _read(gzopen(path))
scenario = _read(gzopen(path))
elseif endswith(path, ".json")
scenario = _read(open(path))
else
return _read(open(path))
error("Unsupported input format")
end
return scenario
end
function _read(file::IO)::UnitCommitmentInstance
function _read(file::IO)::UnitCommitmentScenario
return _from_json(
JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)),
)
@ -77,7 +127,7 @@ function _read_json(path::String)::OrderedDict
return JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing))
end
function _from_json(json; repair = true)
function _from_json(json; repair = true)::UnitCommitmentScenario
_migrate(json)
units = Unit[]
buses = Bus[]
@ -102,6 +152,11 @@ function _from_json(json; repair = true)
time_multiplier = 60 ÷ time_step
T = time_horizon * time_multiplier
probability = json["Parameters"]["Scenario weight"]
probability !== nothing || (probability = 1)
scenario_name = json["Parameters"]["Scenario name"]
scenario_name !== nothing || (scenario_name = "")
name_to_bus = Dict{String,Bus}()
name_to_line = Dict{String,TransmissionLine}()
name_to_unit = Dict{String,Unit}()
@ -118,15 +173,6 @@ function _from_json(json; repair = true)
json["Parameters"]["Power balance penalty (\$/MW)"],
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
for (bus_name, dict) in json["Buses"]
@ -308,7 +354,9 @@ function _from_json(json; repair = true)
end
end
instance = UnitCommitmentInstance(
scenario = UnitCommitmentScenario(
name = scenario_name,
probability = probability,
buses_by_name = Dict(b.name => b for b in buses),
buses = buses,
contingencies_by_name = Dict(c.name => c for c in contingencies),
@ -320,14 +368,14 @@ function _from_json(json; repair = true)
price_sensitive_loads = loads,
reserves = reserves,
reserves_by_name = name_to_reserve,
shortfall_penalty = shortfall_penalty,
flexiramp_shortfall_penalty = flexiramp_shortfall_penalty,
time = T,
units_by_name = Dict(g.name => g for g in units),
units = units,
isf = spzeros(Float64, length(lines), length(buses) - 1),
lodf = spzeros(Float64, length(lines), length(lines)),
)
if repair
UnitCommitment.repair!(instance)
UnitCommitment.repair!(scenario)
end
return instance
return scenario
end

30
src/instance/structs.jl
Unescape View File

@ -73,35 +73,41 @@ mutable struct PriceSensitiveLoad
revenue::Vector{Float64}
end
Base.@kwdef mutable struct UnitCommitmentInstance
Base.@kwdef mutable struct UnitCommitmentScenario
buses_by_name::Dict{AbstractString,Bus}
buses::Vector{Bus}
contingencies_by_name::Dict{AbstractString,Contingency}
contingencies::Vector{Contingency}
isf::Array{Float64,2}
lines_by_name::Dict{AbstractString,TransmissionLine}
lines::Vector{TransmissionLine}
lodf::Array{Float64,2}
name::String
power_balance_penalty::Vector{Float64}
price_sensitive_loads_by_name::Dict{AbstractString,PriceSensitiveLoad}
price_sensitive_loads::Vector{PriceSensitiveLoad}
reserves::Vector{Reserve}
probability::Float64
reserves_by_name::Dict{AbstractString,Reserve}
shortfall_penalty::Vector{Float64}
flexiramp_shortfall_penalty::Vector{Float64}
reserves::Vector{Reserve}
time::Int
units_by_name::Dict{AbstractString,Unit}
units::Vector{Unit}
end
Base.@kwdef mutable struct UnitCommitmentInstance
time::Int
scenarios::Vector{UnitCommitmentScenario}
end
function Base.show(io::IO, instance::UnitCommitmentInstance)
sc = instance.scenarios[1]
print(io, "UnitCommitmentInstance(")
print(io, "$(length(instance.units)) units, ")
print(io, "$(length(instance.buses)) buses, ")
print(io, "$(length(instance.lines)) lines, ")
print(io, "$(length(instance.contingencies)) contingencies, ")
print(
io,
"$(length(instance.price_sensitive_loads)) price sensitive loads, ",
)
print(io, "$(length(instance.scenarios)) scenarios, ")
print(io, "$(length(sc.units)) units, ")
print(io, "$(length(sc.buses)) buses, ")
print(io, "$(length(sc.lines)) lines, ")
print(io, "$(length(sc.contingencies)) contingencies, ")
print(io, "$(length(sc.price_sensitive_loads)) price sensitive loads, ")
print(io, "$(instance.time) time steps")
print(io, ")")
return

31
src/model/build.jl
Unescape View File

@ -77,20 +77,27 @@ function build_model(;
end
model[:obj] = AffExpr()
model[:instance] = instance
_setup_transmission(model, formulation.transmission)
for l in instance.lines
_add_transmission_line!(model, l, formulation.transmission)
for g in instance.scenarios[1].units
_add_unit_commitment!(model, g, formulation)
end
for b in instance.buses
_add_bus!(model, b)
for sc in instance.scenarios
@info "Building scenario $(sc.name) with " *
"probability $(sc.probability)"
_setup_transmission(formulation.transmission, sc)
for l in sc.lines
_add_transmission_line!(model, l, formulation.transmission, sc)
end
for b in sc.buses
_add_bus!(model, b, sc)
end
for ps in sc.price_sensitive_loads
_add_price_sensitive_load!(model, ps, sc)
end
for g in sc.units
_add_unit_dispatch!(model, g, formulation, sc)
end
_add_system_wide_eqs!(model, sc)
end
for g in instance.units
_add_unit!(model, g, formulation)
end
for ps in instance.price_sensitive_loads
_add_price_sensitive_load!(model, ps)
end
_add_system_wide_eqs!(model)
@objective(model, Min, model[:obj])
end
@info @sprintf("Built model in %.2f seconds", time_model)

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

@ -8,6 +8,7 @@ function _add_ramp_eqs!(
formulation_prod_vars::Gar1962.ProdVars,
formulation_ramping::ArrCon2000.Ramping,
formulation_status_vars::Gar1962.StatusVars,
sc::UnitCommitmentScenario,
)::Nothing
# TODO: Move upper case constants to model[:instance]
RESERVES_WHEN_START_UP = true
@ -22,7 +23,7 @@ function _add_ramp_eqs!(
eq_ramp_down = _init(model, :eq_ramp_down)
eq_ramp_up = _init(model, :eq_ramp_up)
is_initially_on = (g.initial_status > 0)
reserve = _total_reserves(model, g)
reserve = _total_reserves(model, g, sc)
# Gar1962.ProdVars
prod_above = model[:prod_above]
@ -37,10 +38,10 @@ function _add_ramp_eqs!(
if t == 1
if is_initially_on
# min power is _not_ multiplied by is_on because if !is_on, then ramp up is irrelevant
eq_ramp_up[gn, t] = @constraint(
eq_ramp_up[sc.name, gn, t] = @constraint(
model,
g.min_power[t] +
prod_above[gn, t] +
prod_above[sc.name, gn, t] +
(RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) <=
g.initial_power + RU
)
@ -48,16 +49,16 @@ function _add_ramp_eqs!(
else
max_prod_this_period =
g.min_power[t] * is_on[gn, t] +
prod_above[gn, t] +
prod_above[sc.name, gn, t] +
(
RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ?
reserve[t] : 0.0
)
min_prod_last_period =
g.min_power[t-1] * is_on[gn, t-1] + prod_above[gn, t-1]
g.min_power[t-1] * is_on[gn, t-1] + prod_above[sc.name, gn, t-1]
# Equation (24) in Kneuven et al. (2020)
eq_ramp_up[gn, t] = @constraint(
eq_ramp_up[sc.name, gn, t] = @constraint(
model,
max_prod_this_period - min_prod_last_period <=
RU * is_on[gn, t-1] + SU * switch_on[gn, t]
@ -71,24 +72,25 @@ function _add_ramp_eqs!(
# min_power + RD < initial_power < SD
# then the generator should be able to shut down at time t = 1,
# but the constraint below will force the unit to produce power
eq_ramp_down[gn, t] = @constraint(
eq_ramp_down[sc.name, gn, t] = @constraint(
model,
g.initial_power - (g.min_power[t] + prod_above[gn, t]) <= RD
g.initial_power -
(g.min_power[t] + prod_above[sc.name, gn, t]) <= RD
)
end
else
max_prod_last_period =
g.min_power[t-1] * is_on[gn, t-1] +
prod_above[gn, t-1] +
prod_above[sc.name, gn, t-1] +
(
RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ?
reserve[t-1] : 0.0
)
min_prod_this_period =
g.min_power[t] * is_on[gn, t] + prod_above[gn, t]
g.min_power[t] * is_on[gn, t] + prod_above[sc.name, gn, t]
# Equation (25) in Kneuven et al. (2020)
eq_ramp_down[gn, t] = @constraint(
eq_ramp_down[sc.name, gn, t] = @constraint(
model,
max_prod_last_period - min_prod_this_period <=
RD * is_on[gn, t] + SD * switch_off[gn, t]

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

@ -8,6 +8,7 @@ function _add_production_piecewise_linear_eqs!(
formulation_prod_vars::Gar1962.ProdVars,
formulation_pwl_costs::CarArr2006.PwlCosts,
formulation_status_vars::StatusVarsFormulation,
sc::UnitCommitmentScenario,
)::Nothing
eq_prod_above_def = _init(model, :eq_prod_above_def)
eq_segprod_limit = _init(model, :eq_segprod_limit)
@ -26,28 +27,32 @@ function _add_production_piecewise_linear_eqs!(
# difference between max power for segments k and k-1 so the
# value of cost_segments[k].mw[t] is the max production *for
# that segment*
eq_segprod_limit[gn, t, k] = @constraint(
eq_segprod_limit[sc.name, gn, t, k] = @constraint(
model,
segprod[gn, t, k] <= g.cost_segments[k].mw[t]
segprod[sc.name, gn, t, k] <= g.cost_segments[k].mw[t]
)
# Also add this as an explicit upper bound on segprod to make the
# solver's work a bit easier
set_upper_bound(segprod[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
# Equation (43) in Kneuven et al. (2020)
eq_prod_above_def[gn, t] = @constraint(
eq_prod_above_def[sc.name, gn, t] = @constraint(
model,
prod_above[gn, t] == sum(segprod[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
# Equation (44) in Kneuven et al. (2020)
add_to_expression!(
model[:obj],
segprod[gn, t, k],
g.cost_segments[k].cost[t],
segprod[sc.name, gn, t, k],
sc.probability * g.cost_segments[k].cost[t],
)
end
end

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

@ -8,6 +8,7 @@ function _add_ramp_eqs!(
formulation_prod_vars::Gar1962.ProdVars,
formulation_ramping::DamKucRajAta2016.Ramping,
formulation_status_vars::Gar1962.StatusVars,
sc::UnitCommitmentScenario,
)::Nothing
# TODO: Move upper case constants to model[:instance]
RESERVES_WHEN_START_UP = true
@ -23,7 +24,7 @@ function _add_ramp_eqs!(
gn = g.name
eq_str_ramp_down = _init(model, :eq_str_ramp_down)
eq_str_ramp_up = _init(model, :eq_str_ramp_up)
reserve = _total_reserves(model, g)
reserve = _total_reserves(model, g, sc)
# Gar1962.ProdVars
prod_above = model[:prod_above]
@ -48,15 +49,15 @@ function _add_ramp_eqs!(
# end
max_prod_this_period =
prod_above[gn, t] +
prod_above[sc.name, gn, t] +
(RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0)
min_prod_last_period = 0.0
if t > 1 && time_invariant
min_prod_last_period = prod_above[gn, t-1]
min_prod_last_period = prod_above[sc.name, gn, t-1]
# Equation (35) in Kneuven et al. (2020)
# Sparser version of (24)
eq_str_ramp_up[gn, t] = @constraint(
eq_str_ramp_up[sc.name, gn, t] = @constraint(
model,
max_prod_this_period - min_prod_last_period <=
(SU - g.min_power[t] - RU) * switch_on[gn, t] +
@ -65,7 +66,8 @@ function _add_ramp_eqs!(
elseif (t == 1 && is_initially_on) || (t > 1 && !time_invariant)
if t > 1
min_prod_last_period =
prod_above[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
min_prod_last_period = max(g.initial_power, 0.0)
end
@ -76,7 +78,7 @@ function _add_ramp_eqs!(
# Modified version of equation (35) in Kneuven et al. (2020)
# Equivalent to (24)
eq_str_ramp_up[gn, t] = @constraint(
eq_str_ramp_up[sc.name, gn, t] = @constraint(
model,
max_prod_this_period - min_prod_last_period <=
(SU - RU) * switch_on[gn, t] + RU * is_on[gn, t]
@ -88,7 +90,7 @@ function _add_ramp_eqs!(
t > 1 && (RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN) ?
reserve[t-1] : 0.0
)
min_prod_this_period = prod_above[gn, t]
min_prod_this_period = prod_above[sc.name, gn, t]
on_last_period = 0.0
if t > 1
on_last_period = is_on[gn, t-1]
@ -98,7 +100,7 @@ function _add_ramp_eqs!(
if t > 1 && time_invariant
# Equation (36) in Kneuven et al. (2020)
eq_str_ramp_down[gn, t] = @constraint(
eq_str_ramp_down[sc.name, gn, t] = @constraint(
model,
max_prod_last_period - min_prod_this_period <=
(SD - g.min_power[t] - RD) * switch_off[gn, t] +
@ -110,7 +112,7 @@ function _add_ramp_eqs!(
# Modified version of equation (36) in Kneuven et al. (2020)
# Equivalent to (25)
eq_str_ramp_down[gn, t] = @constraint(
eq_str_ramp_down[sc.name, gn, t] = @constraint(
model,
max_prod_last_period - min_prod_this_period <=
(SD - RD) * switch_off[gn, t] + RD * on_last_period

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

@ -6,14 +6,15 @@ function _add_production_vars!(
model::JuMP.Model,
g::Unit,
formulation_prod_vars::Gar1962.ProdVars,
sc::UnitCommitmentScenario,
)::Nothing
prod_above = _init(model, :prod_above)
segprod = _init(model, :segprod)
for t in 1:model[:instance].time
for k in 1:length(g.cost_segments)
segprod[g.name, t, k] = @variable(model, lower_bound = 0)
segprod[sc.name, g.name, t, k] = @variable(model, lower_bound = 0)
end
prod_above[g.name, t] = @variable(model, lower_bound = 0)
prod_above[sc.name, g.name, t] = @variable(model, lower_bound = 0)
end
return
end
@ -22,16 +23,16 @@ function _add_production_limit_eqs!(
model::JuMP.Model,
g::Unit,
formulation_prod_vars::Gar1962.ProdVars,
sc::UnitCommitmentScenario,
)::Nothing
eq_prod_limit = _init(model, :eq_prod_limit)
is_on = model[:is_on]
prod_above = model[:prod_above]
reserve = _total_reserves(model, g)
reserve = _total_reserves(model, g, sc)
gn = g.name
for t in 1:model[:instance].time
# Objective function terms for production costs
# Part of (69) of Kneuven et al. (2020) as C^R_g * u_g(t) term
add_to_expression!(model[:obj], is_on[gn, t], g.min_power_cost[t])
# Production limit
# Equation (18) in Kneuven et al. (2020)
@ -42,9 +43,10 @@ function _add_production_limit_eqs!(
if power_diff < 1e-7
power_diff = 0.0
end
eq_prod_limit[gn, t] = @constraint(
eq_prod_limit[sc.name, gn, t] = @constraint(
model,
prod_above[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

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

@ -8,6 +8,7 @@ function _add_production_piecewise_linear_eqs!(
formulation_prod_vars::Gar1962.ProdVars,
formulation_pwl_costs::Gar1962.PwlCosts,
formulation_status_vars::Gar1962.StatusVars,
sc::UnitCommitmentScenario,
)::Nothing
eq_prod_above_def = _init(model, :eq_prod_above_def)
eq_segprod_limit = _init(model, :eq_segprod_limit)
@ -24,9 +25,10 @@ function _add_production_piecewise_linear_eqs!(
for t in 1:model[:instance].time
# Definition of production
# Equation (43) in Kneuven et al. (2020)
eq_prod_above_def[gn, t] = @constraint(
eq_prod_above_def[sc.name, gn, t] = @constraint(
model,
prod_above[gn, t] == sum(segprod[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
@ -37,21 +39,25 @@ function _add_production_piecewise_linear_eqs!(
# difference between max power for segments k and k-1 so the
# value of cost_segments[k].mw[t] is the max production *for
# that segment*
eq_segprod_limit[gn, t, k] = @constraint(
eq_segprod_limit[sc.name, gn, t, k] = @constraint(
model,
segprod[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
# solver's work a bit easier
set_upper_bound(segprod[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
# Equation (44) in Kneuven et al. (2020)
add_to_expression!(
model[:obj],
segprod[gn, t, k],
g.cost_segments[k].cost[t],
segprod[sc.name, gn, t, k],
sc.probability * g.cost_segments[k].cost[t],
)
end
end

1
src/model/formulations/Gar1962/status.jl
Unescape View File

@ -20,6 +20,7 @@ function _add_status_vars!(
switch_on[g.name, t] = @variable(model, binary = true)
switch_off[g.name, t] = @variable(model, binary = true)
end
add_to_expression!(model[:obj], is_on[g.name, t], g.min_power_cost[t])
end
return
end

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

@ -8,6 +8,7 @@ function _add_production_piecewise_linear_eqs!(
formulation_prod_vars::Gar1962.ProdVars,
formulation_pwl_costs::KnuOstWat2018.PwlCosts,
formulation_status_vars::Gar1962.StatusVars,
sc::UnitCommitmentScenario,
)::Nothing
eq_prod_above_def = _init(model, :eq_prod_above_def)
eq_segprod_limit_a = _init(model, :eq_segprod_limit_a)
@ -58,27 +59,27 @@ function _add_production_piecewise_linear_eqs!(
if g.min_uptime > 1
# Equation (46) in Kneuven et al. (2020)
eq_segprod_limit_a[gn, t, k] = @constraint(
eq_segprod_limit_a[sc.name, gn, t, k] = @constraint(
model,
segprod[gn, t, k] <=
segprod[sc.name, gn, t, k] <=
g.cost_segments[k].mw[t] * is_on[gn, t] -
Cv * switch_on[gn, t] -
(t < T ? Cw * switch_off[gn, t+1] : 0.0)
)
else
# Equation (47a)/(48a) in Kneuven et al. (2020)
eq_segprod_limit_b[gn, t, k] = @constraint(
eq_segprod_limit_b[sc.name, gn, t, k] = @constraint(
model,
segprod[gn, t, k] <=
segprod[sc.name, gn, t, k] <=
g.cost_segments[k].mw[t] * is_on[gn, t] -
Cv * switch_on[gn, t] -
(t < T ? max(0, Cv - Cw) * switch_off[gn, t+1] : 0.0)
)
# Equation (47b)/(48b) in Kneuven et al. (2020)
eq_segprod_limit_c[gn, t, k] = @constraint(
eq_segprod_limit_c[sc.name, gn, t, k] = @constraint(
model,
segprod[gn, t, k] <=
segprod[sc.name, gn, t, k] <=
g.cost_segments[k].mw[t] * is_on[gn, t] -
max(0, Cw - Cv) * switch_on[gn, t] -
(t < T ? Cw * switch_off[gn, t+1] : 0.0)
@ -87,22 +88,26 @@ function _add_production_piecewise_linear_eqs!(
# Definition of production
# Equation (43) in Kneuven et al. (2020)
eq_prod_above_def[gn, t] = @constraint(
eq_prod_above_def[sc.name, gn, t] = @constraint(
model,
prod_above[gn, t] == sum(segprod[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
# Equation (44) in Kneuven et al. (2020)
add_to_expression!(
model[:obj],
segprod[gn, t, k],
segprod[sc.name, gn, t, k],
g.cost_segments[k].cost[t],
)
# Also add an explicit upper bound on segprod to make the solver's
# work a bit easier
set_upper_bound(segprod[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

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

@ -8,6 +8,7 @@ function _add_ramp_eqs!(
formulation_prod_vars::Gar1962.ProdVars,
formulation_ramping::MorLatRam2013.Ramping,
formulation_status_vars::Gar1962.StatusVars,
sc::UnitCommitmentScenario,
)::Nothing
# TODO: Move upper case constants to model[:instance]
RESERVES_WHEN_START_UP = true
@ -22,7 +23,7 @@ function _add_ramp_eqs!(
gn = g.name
eq_ramp_down = _init(model, :eq_ramp_down)
eq_ramp_up = _init(model, :eq_str_ramp_up)
reserve = _total_reserves(model, g)
reserve = _total_reserves(model, g, sc)
# Gar1962.ProdVars
prod_above = model[:prod_above]
@ -39,10 +40,10 @@ function _add_ramp_eqs!(
# Ramp up limit
if t == 1
if is_initially_on
eq_ramp_up[gn, t] = @constraint(
eq_ramp_up[sc.name, gn, t] = @constraint(
model,
g.min_power[t] +
prod_above[gn, t] +
prod_above[sc.name, gn, t] +
(RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) <=
g.initial_power + RU
)
@ -58,13 +59,14 @@ function _add_ramp_eqs!(
SU = g.startup_limit
max_prod_this_period =
g.min_power[t] * is_on[gn, t] +
prod_above[gn, t] +
prod_above[sc.name, gn, t] +
(
RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ?
reserve[t] : 0.0
)
min_prod_last_period =
g.min_power[t-1] * is_on[gn, t-1] + prod_above[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(
model,
max_prod_this_period - min_prod_last_period <=
@ -74,11 +76,11 @@ function _add_ramp_eqs!(
# Equation (26) in Kneuven et al. (2020)
# TODO: what if RU < SU? places too stringent upper bound
# prod_above[gn, t] when starting up, and creates diff with (24).
eq_ramp_up[gn, t] = @constraint(
eq_ramp_up[sc.name, gn, t] = @constraint(
model,
prod_above[gn, t] +
prod_above[sc.name, gn, t] +
(RESERVES_WHEN_RAMP_UP ? reserve[t] : 0.0) -
prod_above[gn, t-1] <= RU
prod_above[sc.name, gn, t-1] <= RU
)
end
end
@ -90,9 +92,10 @@ function _add_ramp_eqs!(
# min_power + RD < initial_power < SD
# then the generator should be able to shut down at time t = 1,
# but the constraint below will force the unit to produce power
eq_ramp_down[gn, t] = @constraint(
eq_ramp_down[sc.name, gn, t] = @constraint(
model,
g.initial_power - (g.min_power[t] + prod_above[gn, t]) <= RD
g.initial_power -
(g.min_power[t] + prod_above[sc.name, gn, t]) <= RD
)
end
else
@ -102,13 +105,13 @@ function _add_ramp_eqs!(
SD = g.shutdown_limit
max_prod_last_period =
g.min_power[t-1] * is_on[gn, t-1] +
prod_above[gn, t-1] +
prod_above[sc.name, gn, t-1] +
(
RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ?
reserve[t-1] : 0.0
)
min_prod_this_period =
g.min_power[t] * is_on[gn, t] + prod_above[gn, t]
g.min_power[t] * is_on[gn, t] + prod_above[sc.name, gn, t]
eq_ramp_down[gn, t] = @constraint(
model,
max_prod_last_period - min_prod_this_period <=
@ -118,11 +121,11 @@ function _add_ramp_eqs!(
# Equation (27) in Kneuven et al. (2020)
# TODO: Similar to above, what to do if shutting down in time t
# and RD < SD? There is a difference with (25).
eq_ramp_down[gn, t] = @constraint(
eq_ramp_down[sc.name, gn, t] = @constraint(
model,
prod_above[gn, t-1] +
prod_above[sc.name, gn, t-1] +
(RESERVES_WHEN_RAMP_DOWN ? reserve[t-1] : 0.0) -
prod_above[gn, t] <= RD
prod_above[sc.name, gn, t] <= RD
)
end
end

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

@ -8,11 +8,12 @@ function _add_ramp_eqs!(
formulation_prod_vars::Gar1962.ProdVars,
formulation_ramping::PanGua2016.Ramping,
formulation_status_vars::Gar1962.StatusVars,
sc::UnitCommitmentScenario,
)::Nothing
# TODO: Move upper case constants to model[:instance]
RESERVES_WHEN_SHUT_DOWN = true
gn = g.name
reserve = _total_reserves(model, g)
reserve = _total_reserves(model, g, sc)
eq_str_prod_limit = _init(model, :eq_str_prod_limit)
eq_prod_limit_ramp_up_extra_period =
_init(model, :eq_prod_limit_ramp_up_extra_period)
@ -52,9 +53,9 @@ function _add_ramp_eqs!(
# Generalization of (20)
# Necessary that if any of the switch_on = 1 in the sum,
# then switch_off[gn, t+1] = 0
eq_str_prod_limit[gn, t] = @constraint(
eq_str_prod_limit[sc.name, gn, t] = @constraint(
model,
prod_above[gn, t] +
prod_above[sc.name, gn, t] +
g.min_power[t] * is_on[gn, t] +
reserve[t] <=
Pbar * is_on[gn, t] -
@ -67,16 +68,17 @@ function _add_ramp_eqs!(
if UT - 2 < TRU
# Equation (40) in Kneuven et al. (2020)
# Covers an additional time period of the ramp-up trajectory, compared to (38)
eq_prod_limit_ramp_up_extra_period[gn, t] = @constraint(
model,
prod_above[gn, t] +
g.min_power[t] * is_on[gn, t] +
reserve[t] <=
Pbar * is_on[gn, t] - sum(
(Pbar - (SU + i * RU)) * switch_on[gn, t-i] for
i in 0:min(UT - 1, TRU, t - 1)
eq_prod_limit_ramp_up_extra_period[sc.name, gn, t] =
@constraint(
model,
prod_above[sc.name, gn, t] +
g.min_power[t] * is_on[gn, t] +
reserve[t] <=
Pbar * is_on[gn, t] - sum(
(Pbar - (SU + i * RU)) * switch_on[gn, t-i] for
i in 0:min(UT - 1, TRU, t - 1)
)
)
)
end
# Add in shutdown trajectory if KSD >= 0 (else this is dominated by (38))
@ -84,9 +86,9 @@ function _add_ramp_eqs!(
if KSD > 0
KSU = min(TRU, UT - 2 - KSD, t - 1)
# Equation (41) in Kneuven et al. (2020)
eq_prod_limit_shutdown_trajectory[gn, t] = @constraint(
eq_prod_limit_shutdown_trajectory[sc.name, gn, t] = @constraint(
model,
prod_above[gn, t] +
prod_above[sc.name, gn, t] +
g.min_power[t] * is_on[gn, t] +
(RESERVES_WHEN_SHUT_DOWN ? reserve[t] : 0.0) <=
Pbar * is_on[gn, t] - sum(

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

@ -8,6 +8,7 @@ function _add_ramp_eqs!(
::Gar1962.ProdVars,
::WanHob2016.Ramping,
::Gar1962.StatusVars,
sc::UnitCommitmentScenario,
)::Nothing
is_initially_on = (g.initial_status > 0)
SU = g.startup_limit
@ -38,41 +39,43 @@ function _add_ramp_eqs!(
for t in 1:model[:instance].time
@constraint(
model,
prod_above[gn, t] + (is_on[gn, t] * minp[t]) <= mfg[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)
@constraint(model, mfg[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
@constraint(
model,
minp[t] * (is_on[gn, t+1] + is_on[gn, t] - 1) <=
prod_above[gn, t] - dwflexiramp[rn, gn, t] +
(is_on[gn, t] * minp[t])
prod_above[sc.name, gn, t] -
dwflexiramp[sc.name, rn, gn, t] + (is_on[gn, t] * minp[t])
) # first inequality of Eq. (20) in Wang & Hobbs (2016)
@constraint(
model,
prod_above[gn, t] - dwflexiramp[rn, gn, t] +
prod_above[sc.name, gn, t] -
dwflexiramp[sc.name, rn, gn, t] +
(is_on[gn, t] * minp[t]) <=
mfg[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)
@constraint(
model,
minp[t] * (is_on[gn, t+1] + is_on[gn, t] - 1) <=
prod_above[gn, t] +
upflexiramp[rn, gn, t] +
prod_above[sc.name, gn, t] +
upflexiramp[sc.name, rn, gn, t] +
(is_on[gn, t] * minp[t])
) # first inequality of Eq. (21) in Wang & Hobbs (2016)
@constraint(
model,
prod_above[gn, t] +
upflexiramp[rn, gn, t] +
prod_above[sc.name, gn, t] +
upflexiramp[sc.name, rn, gn, t] +
(is_on[gn, t] * minp[t]) <=
mfg[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)
if t != 1
@constraint(
model,
mfg[rn, gn, t] <=
prod_above[gn, t-1] +
mfg[sc.name, gn, t] <=
prod_above[sc.name, gn, t-1] +
(is_on[gn, t-1] * minp[t]) +
(RU * is_on[gn, t-1]) +
(SU * (is_on[gn, t] - is_on[gn, t-1])) +
@ -80,8 +83,13 @@ function _add_ramp_eqs!(
) # Eq. (23) in Wang & Hobbs (2016)
@constraint(
model,
(prod_above[gn, t-1] + (is_on[gn, t-1] * minp[t])) -
(prod_above[gn, t] + (is_on[gn, t] * minp[t])) <=
(
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] +
SD * (is_on[gn, t-1] - is_on[gn, t]) +
maxp[t] * (1 - is_on[gn, t-1])
@ -89,7 +97,7 @@ function _add_ramp_eqs!(
else
@constraint(
model,
mfg[rn, gn, t] <=
mfg[sc.name, gn, t] <=
initial_power +
(RU * is_initially_on) +
(SU * (is_on[gn, t] - is_initially_on)) +
@ -97,8 +105,10 @@ function _add_ramp_eqs!(
) # Eq. (23) in Wang & Hobbs (2016) for the first time period
@constraint(
model,
initial_power -
(prod_above[gn, t] + (is_on[gn, t] * minp[t])) <=
initial_power - (
prod_above[sc.name, gn, t] +
(is_on[gn, t] * minp[t])
) <=
RD * is_on[gn, t] +
SD * (is_initially_on - is_on[gn, t]) +
maxp[t] * (1 - is_initially_on)
@ -106,7 +116,7 @@ function _add_ramp_eqs!(
end
@constraint(
model,
mfg[rn, gn, t] <=
mfg[sc.name, gn, t] <=
(SD * (is_on[gn, t] - is_on[gn, t+1])) +
(maxp[t] * is_on[gn, t+1])
) # Eq. (24) in Wang & Hobbs (2016)
@ -114,11 +124,12 @@ function _add_ramp_eqs!(
model,
-RD * is_on[gn, t+1] -
SD * (is_on[gn, t] - is_on[gn, t+1]) -
maxp[t] * (1 - is_on[gn, t]) <= upflexiramp[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)
@constraint(
model,
upflexiramp[rn, gn, t] <=
upflexiramp[sc.name, rn, 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])
@ -126,11 +137,12 @@ function _add_ramp_eqs!(
@constraint(
model,
-RU * is_on[gn, t] - SU * (is_on[gn, t+1] - is_on[gn, t]) -
maxp[t] * (1 - is_on[gn, t+1]) <= dwflexiramp[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)
@constraint(
model,
dwflexiramp[rn, gn, t] <=
dwflexiramp[sc.name, rn, gn, t] <=
RD * is_on[gn, t+1] +
SD * (is_on[gn, t] - is_on[gn, t+1]) +
maxp[t] * (1 - is_on[gn, t])
@ -138,26 +150,27 @@ function _add_ramp_eqs!(
@constraint(
model,
-maxp[t] * is_on[gn, t] + minp[t] * is_on[gn, t+1] <=
upflexiramp[rn, gn, t]
upflexiramp[sc.name, rn, gn, t]
) # first inequality of Eq. (28) in Wang & Hobbs (2016)
@constraint(
model,
upflexiramp[rn, gn, t] <= maxp[t] * is_on[gn, t+1]
upflexiramp[sc.name, rn, gn, t] <= maxp[t] * is_on[gn, t+1]
) # second inequality of Eq. (28) in Wang & Hobbs (2016)
@constraint(
model,
-maxp[t] * is_on[gn, t+1] <= dwflexiramp[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)
@constraint(
model,
dwflexiramp[rn, gn, t] <=
dwflexiramp[sc.name, rn, gn, t] <=
(maxp[t] * is_on[gn, t]) - (minp[t] * is_on[gn, t+1])
) # second inequality of Eq. (29) in Wang & Hobbs (2016)
else
@constraint(
model,
mfg[rn, gn, t] <=
prod_above[gn, t-1] +
mfg[sc.name, gn, t] <=
prod_above[sc.name, gn, t-1] +
(is_on[gn, t-1] * minp[t]) +
(RU * is_on[gn, t-1]) +
(SU * (is_on[gn, t] - is_on[gn, t-1])) +
@ -165,8 +178,11 @@ function _add_ramp_eqs!(
) # Eq. (23) in Wang & Hobbs (2016) for the last time period
@constraint(
model,
(prod_above[gn, t-1] + (is_on[gn, t-1] * minp[t])) -
(prod_above[gn, t] + (is_on[gn, t] * minp[t])) <=
(
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] +
SD * (is_on[gn, t-1] - is_on[gn, t]) +
maxp[t] * (1 - is_on[gn, t-1])

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

@ -2,22 +2,30 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
function _add_bus!(model::JuMP.Model, b::Bus)::Nothing
function _add_bus!(
model::JuMP.Model,
b::Bus,
sc::UnitCommitmentScenario,
)::Nothing
net_injection = _init(model, :expr_net_injection)
curtail = _init(model, :curtail)
for t in 1:model[:instance].time
# Fixed load
net_injection[b.name, t] = AffExpr(-b.load[t])
net_injection[sc.name, b.name, t] = AffExpr(-b.load[t])
# Load curtailment
curtail[b.name, t] =
curtail[sc.name, b.name, t] =
@variable(model, lower_bound = 0, upper_bound = b.load[t])
add_to_expression!(net_injection[b.name, t], curtail[b.name, t], 1.0)
add_to_expression!(
net_injection[sc.name, b.name, t],
curtail[sc.name, b.name, t],
1.0,
)
add_to_expression!(
model[:obj],
curtail[b.name, t],
model[:instance].power_balance_penalty[t],
curtail[sc.name, b.name, t],
sc.power_balance_penalty[t] * sc.probability,
)
end
return

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

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

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

@ -5,21 +5,26 @@
function _add_price_sensitive_load!(
model::JuMP.Model,
ps::PriceSensitiveLoad,
sc::UnitCommitmentScenario,
)::Nothing
loads = _init(model, :loads)
net_injection = _init(model, :expr_net_injection)
for t in 1:model[:instance].time
# Decision variable
loads[ps.name, t] =
loads[sc.name, ps.name, t] =
@variable(model, lower_bound = 0, upper_bound = ps.demand[t])
# Objective function terms
add_to_expression!(model[:obj], loads[ps.name, t], -ps.revenue[t])
add_to_expression!(
model[:obj],
loads[sc.name, ps.name, t],
-ps.revenue[t] * sc.probability,
)
# Net injection
add_to_expression!(
net_injection[ps.bus.name, t],
loads[ps.name, t],
net_injection[sc.name, ps.bus.name, t],
loads[sc.name, ps.name, t],
-1.0,
)
end

2
src/model/formulations/base/sensitivity.jl
Unescape View File

@ -18,7 +18,7 @@ function _injection_shift_factors(;
lines::Array{TransmissionLine},
)
susceptance = _susceptance_matrix(lines)
incidence = _reduced_incidence_matrix(lines = lines, buses = buses)
incidence = _reduced_incidence_matrix(buses = buses, lines = lines)
laplacian = transpose(incidence) * susceptance * incidence
isf = susceptance * incidence * inv(Array(laplacian))
return isf

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

@ -2,54 +2,68 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
function _add_system_wide_eqs!(model::JuMP.Model)::Nothing
_add_net_injection_eqs!(model)
_add_spinning_reserve_eqs!(model)
_add_flexiramp_reserve_eqs!(model)
function _add_system_wide_eqs!(
model::JuMP.Model,
sc::UnitCommitmentScenario,
)::Nothing
_add_net_injection_eqs!(model, sc)
_add_spinning_reserve_eqs!(model, sc)
_add_flexiramp_reserve_eqs!(model, sc)
return
end
function _add_net_injection_eqs!(model::JuMP.Model)::Nothing
function _add_net_injection_eqs!(
model::JuMP.Model,
sc::UnitCommitmentScenario,
)::Nothing
T = model[:instance].time
net_injection = _init(model, :net_injection)
eq_net_injection = _init(model, :eq_net_injection)
eq_power_balance = _init(model, :eq_power_balance)
for t in 1:T, b in model[:instance].buses
n = net_injection[b.name, t] = @variable(model)
eq_net_injection[b.name, t] =
@constraint(model, -n + model[:expr_net_injection][b.name, t] == 0)
for t in 1:T, b in sc.buses
n = net_injection[sc.name, b.name, t] = @variable(model)
eq_net_injection[sc.name, b.name, t] = @constraint(
model,
-n + model[:expr_net_injection][sc.name, b.name, t] == 0
)
end
for t in 1:T
eq_power_balance[t] = @constraint(
eq_power_balance[sc.name, t] = @constraint(
model,
sum(net_injection[b.name, t] for b in model[:instance].buses) == 0
sum(net_injection[sc.name, b.name, t] for b in sc.buses) == 0
)
end
return
end
function _add_spinning_reserve_eqs!(model::JuMP.Model)::Nothing
instance = model[:instance]
function _add_spinning_reserve_eqs!(
model::JuMP.Model,
sc::UnitCommitmentScenario,
)::Nothing
T = model[:instance].time
eq_min_spinning_reserve = _init(model, :eq_min_spinning_reserve)
for r in instance.reserves
for r in sc.reserves
r.type == "spinning" || continue
for t in 1:instance.time
for t in 1:T
# Equation (68) in Kneuven et al. (2020)
# As in Morales-España et al. (2013a)
# Akin to the alternative formulation with max_power_avail
# from Carrión and Arroyo (2006) and Ostrowski et al. (2012)
eq_min_spinning_reserve[r.name, t] = @constraint(
eq_min_spinning_reserve[sc.name, r.name, t] = @constraint(
model,
sum(model[:reserve][r.name, g.name, t] for g in r.units) +
model[:reserve_shortfall][r.name, t] >= r.amount[t]
sum(
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
if r.shortfall_penalty >= 0
add_to_expression!(
model[:obj],
r.shortfall_penalty,
model[:reserve_shortfall][r.name, t],
r.shortfall_penalty * sc.probability,
model[:reserve_shortfall][sc.name, r.name, t],
)
end
end
@ -57,7 +71,10 @@ function _add_spinning_reserve_eqs!(model::JuMP.Model)::Nothing
return
end
function _add_flexiramp_reserve_eqs!(model::JuMP.Model)::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
# 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
@ -65,29 +82,37 @@ function _add_flexiramp_reserve_eqs!(model::JuMP.Model)::Nothing
# objective function.
eq_min_upflexiramp = _init(model, :eq_min_upflexiramp)
eq_min_dwflexiramp = _init(model, :eq_min_dwflexiramp)
instance = model[:instance]
for r in instance.reserves
T = model[:instance].time
for r in sc.reserves
r.type == "flexiramp" || continue
for t in 1:instance.time
for t in 1:T
# Eq. (17) in Wang & Hobbs (2016)
eq_min_upflexiramp[r.name, t] = @constraint(
eq_min_upflexiramp[sc.name, r.name, t] = @constraint(
model,
sum(model[:upflexiramp][r.name, g.name, t] for g in r.units) + model[:upflexiramp_shortfall][r.name, t] >= r.amount[t]
sum(
model[:upflexiramp][sc.name, r.name, g.name, t] for
g in r.units
) + model[:upflexiramp_shortfall][sc.name, r.name, t] >=
r.amount[t]
)
# Eq. (18) in Wang & Hobbs (2016)
eq_min_dwflexiramp[r.name, t] = @constraint(
eq_min_dwflexiramp[sc.name, r.name, t] = @constraint(
model,
sum(model[:dwflexiramp][r.name, g.name, t] for g in r.units) + model[:dwflexiramp_shortfall][r.name, t] >= r.amount[t]
sum(
model[:dwflexiramp][sc.name, r.name, g.name, t] for
g in r.units
) + model[:dwflexiramp_shortfall][sc.name, r.name, t] >=
r.amount[t]
)
# Account for flexiramp shortfall contribution to objective
if r.shortfall_penalty >= 0
add_to_expression!(
model[:obj],
r.shortfall_penalty,
r.shortfall_penalty * sc.probability,
(
model[:upflexiramp_shortfall][r.name, t] +
model[:dwflexiramp_shortfall][r.name, t]
model[:upflexiramp_shortfall][sc.name, r.name, t] +
model[:dwflexiramp_shortfall][sc.name, r.name, t]
),
)
end

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

@ -2,7 +2,13 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
function _add_unit!(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)
error("Partially must-run units are not currently supported")
end
@ -11,22 +17,40 @@ function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
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_startup_cost_eqs!(model, g, formulation.startup_costs)
_add_status_eqs!(model, g, formulation.status_vars)
return
end
# Function for adding variables, constraints, and objective function terms
# related to the continuous dispatch decisions of units
function _add_unit_dispatch!(
model::JuMP.Model,
g::Unit,
formulation::Formulation,
sc::UnitCommitmentScenario,
)
# Variables
_add_production_vars!(model, g, formulation.prod_vars, sc)
_add_spinning_reserve_vars!(model, g, sc)
_add_flexiramp_reserve_vars!(model, g, sc)
# Constraints and objective function
_add_net_injection_eqs!(model, g, sc)
_add_production_limit_eqs!(model, g, formulation.prod_vars, sc)
_add_production_piecewise_linear_eqs!(
model,
g,
formulation.prod_vars,
formulation.pwl_costs,
formulation.status_vars,
sc,
)
_add_ramp_eqs!(
model,
@ -34,26 +58,31 @@ function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
formulation.prod_vars,
formulation.ramping,
formulation.status_vars,
sc,
)
_add_startup_cost_eqs!(model, g, formulation.startup_costs)
_add_startup_shutdown_limit_eqs!(model, g)
_add_status_eqs!(model, g, formulation.status_vars)
_add_startup_shutdown_limit_eqs!(model, g, sc)
return
end
_is_initially_on(g::Unit)::Float64 = (g.initial_status > 0 ? 1.0 : 0.0)
function _add_spinning_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
function _add_spinning_reserve_vars!(
model::JuMP.Model,
g::Unit,
sc::UnitCommitmentScenario,
)::Nothing
reserve = _init(model, :reserve)
reserve_shortfall = _init(model, :reserve_shortfall)
for r in g.reserves
r.type == "spinning" || continue
for t in 1:model[:instance].time
reserve[r.name, g.name, t] = @variable(model, lower_bound = 0)
if (r.name, t) ∉ keys(reserve_shortfall)
reserve_shortfall[r.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)
reserve_shortfall[sc.name, r.name, t] =
@variable(model, lower_bound = 0)
if r.shortfall_penalty < 0
set_upper_bound(reserve_shortfall[r.name, t], 0.0)
set_upper_bound(reserve_shortfall[sc.name, r.name, t], 0.0)
end
end
end
@ -61,27 +90,37 @@ function _add_spinning_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
return
end
function _add_flexiramp_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
function _add_flexiramp_reserve_vars!(
model::JuMP.Model,
g::Unit,
sc::UnitCommitmentScenario,
)::Nothing
upflexiramp = _init(model, :upflexiramp)
upflexiramp_shortfall = _init(model, :upflexiramp_shortfall)
mfg = _init(model, :mfg)
dwflexiramp = _init(model, :dwflexiramp)
dwflexiramp_shortfall = _init(model, :dwflexiramp_shortfall)
for r in g.reserves
r.type == "flexiramp" || continue
for t in 1:model[:instance].time
# maximum feasible generation, \bar{g_{its}} in Wang & Hobbs (2016)
mfg[r.name, g.name, t] = @variable(model, lower_bound = 0)
upflexiramp[r.name, g.name, t] = @variable(model) # up-flexiramp, ur_{it} in Wang & Hobbs (2016)
dwflexiramp[r.name, g.name, t] = @variable(model) # down-flexiramp, dr_{it} in Wang & Hobbs (2016)
if (r.name, t) ∉ keys(upflexiramp_shortfall)
upflexiramp_shortfall[r.name, t] =
for t in 1:model[:instance].time
# maximum feasible generation, \bar{g_{its}} in Wang & Hobbs (2016)
mfg[sc.name, g.name, t] = @variable(model, lower_bound = 0)
for r in g.reserves
r.type == "flexiramp" || continue
upflexiramp[sc.name, r.name, g.name, t] = @variable(model) # up-flexiramp, ur_{it} in Wang & Hobbs (2016)
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(upflexiramp_shortfall)
upflexiramp_shortfall[sc.name, r.name, t] =
@variable(model, lower_bound = 0)
dwflexiramp_shortfall[r.name, t] =
dwflexiramp_shortfall[sc.name, r.name, t] =
@variable(model, lower_bound = 0)
if r.shortfall_penalty < 0
set_upper_bound(upflexiramp_shortfall[r.name, t], 0.0)
set_upper_bound(dwflexiramp_shortfall[r.name, t], 0.0)
set_upper_bound(
upflexiramp_shortfall[sc.name, r.name, t],
0.0,
)
set_upper_bound(
dwflexiramp_shortfall[sc.name, r.name, t],
0.0,
)
end
end
end
@ -99,32 +138,36 @@ function _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing
return
end
function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
function _add_startup_shutdown_limit_eqs!(
model::JuMP.Model,
g::Unit,
sc::UnitCommitmentScenario,
)::Nothing
eq_shutdown_limit = _init(model, :eq_shutdown_limit)
eq_startup_limit = _init(model, :eq_startup_limit)
is_on = model[:is_on]
prod_above = model[:prod_above]
reserve = _total_reserves(model, g)
reserve = _total_reserves(model, g, sc)
switch_off = model[:switch_off]
switch_on = model[:switch_on]
T = model[:instance].time
for t in 1:T
# Startup limit
eq_startup_limit[g.name, t] = @constraint(
eq_startup_limit[sc.name, g.name, t] = @constraint(
model,
prod_above[g.name, t] + reserve[t] <=
prod_above[sc.name, g.name, t] + reserve[t] <=
(g.max_power[t] - g.min_power[t]) * is_on[g.name, t] -
max(0, g.max_power[t] - g.startup_limit) * switch_on[g.name, t]
)
# Shutdown limit
if g.initial_power > g.shutdown_limit
eq_shutdown_limit[g.name, 0] =
eq_shutdown_limit[sc.name, g.name, 0] =
@constraint(model, switch_off[g.name, 1] <= 0)
end
if t < T
eq_shutdown_limit[g.name, t] = @constraint(
eq_shutdown_limit[sc.name, g.name, t] = @constraint(
model,
prod_above[g.name, t] <=
prod_above[sc.name, g.name, t] <=
(g.max_power[t] - g.min_power[t]) * is_on[g.name, t] -
max(0, g.max_power[t] - g.shutdown_limit) *
switch_off[g.name, t+1]
@ -138,43 +181,44 @@ function _add_ramp_eqs!(
model::JuMP.Model,
g::Unit,
formulation::RampingFormulation,
sc::UnitCommitmentScenario,
)::Nothing
prod_above = model[:prod_above]
reserve = _total_reserves(model, g)
reserve = _total_reserves(model, g, sc)
eq_ramp_up = _init(model, :eq_ramp_up)
eq_ramp_down = _init(model, :eq_ramp_down)
for t in 1:model[:instance].time
# Ramp up limit
if t == 1
if _is_initially_on(g) == 1
eq_ramp_up[g.name, t] = @constraint(
eq_ramp_up[sc.name, g.name, t] = @constraint(
model,
prod_above[g.name, t] + reserve[t] <=
prod_above[sc.name, g.name, t] + reserve[t] <=
(g.initial_power - g.min_power[t]) + g.ramp_up_limit
)
end
else
eq_ramp_up[g.name, t] = @constraint(
eq_ramp_up[sc.name, g.name, t] = @constraint(
model,
prod_above[g.name, t] + reserve[t] <=
prod_above[g.name, t-1] + g.ramp_up_limit
prod_above[sc.name, g.name, t] + reserve[t] <=
prod_above[sc.name, g.name, t-1] + g.ramp_up_limit
)
end
# Ramp down limit
if t == 1
if _is_initially_on(g) == 1
eq_ramp_down[g.name, t] = @constraint(
eq_ramp_down[sc.name, g.name, t] = @constraint(
model,
prod_above[g.name, t] >=
prod_above[sc.name, g.name, t] >=
(g.initial_power - g.min_power[t]) - g.ramp_down_limit
)
end
else
eq_ramp_down[g.name, t] = @constraint(
eq_ramp_down[sc.name, g.name, t] = @constraint(
model,
prod_above[g.name, t] >=
prod_above[g.name, t-1] - g.ramp_down_limit
prod_above[sc.name, g.name, t] >=
prod_above[sc.name, g.name, t-1] - g.ramp_down_limit
)
end
end
@ -223,30 +267,37 @@ function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
end
end
function _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing
function _add_net_injection_eqs!(
model::JuMP.Model,
g::Unit,
sc::UnitCommitmentScenario,
)::Nothing
expr_net_injection = model[:expr_net_injection]
for t in 1:model[:instance].time
# Add to net injection expression
add_to_expression!(
expr_net_injection[g.bus.name, t],
model[:prod_above][g.name, t],
expr_net_injection[sc.name, g.bus.name, t],
model[:prod_above][sc.name, g.name, t],
1.0,
)
add_to_expression!(
expr_net_injection[g.bus.name, t],
expr_net_injection[sc.name, g.bus.name, t],
model[:is_on][g.name, t],
g.min_power[t],
)
end
end
function _total_reserves(model, g)::Vector
function _total_reserves(model, g, sc)::Vector
T = model[:instance].time
reserve = [0.0 for _ in 1:T]
spinning_reserves = [r for r in g.reserves if r.type == "spinning"]
if !isempty(spinning_reserves)
reserve += [
sum(model[:reserve][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
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
instance, T = model[:instance], model[:instance].time
"Production (MW)" ∈ keys(solution) ? solution = Dict("s1" => solution) :
nothing
is_on = model[:is_on]
prod_above = model[:prod_above]
reserve = model[:reserve]
for g in instance.units
for t in 1:T
is_on_value = round(solution["Is on"][g.name][t])
prod_value =
round(solution["Production (MW)"][g.name][t], digits = 5)
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 sc in instance.scenarios
for g in sc.units
for t in 1:T
reserve_value = round(
solution["Spinning reserve (MW)"][r.name][g.name][t],
is_on_value = round(solution[sc.name]["Is on"][g.name][t])
prod_value = round(
solution[sc.name]["Production (MW)"][g.name][t],
digits = 5,
)
JuMP.fix(is_on[g.name, t], is_on_value, force = true)
JuMP.fix(
reserve[r.name, g.name, t],
reserve_value,
prod_above[sc.name, g.name, t],
prod_value - is_on_value * g.min_power[t],
force = true,
)
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
return
end

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

@ -5,13 +5,15 @@
function _enforce_transmission(
model::JuMP.Model,
violations::Vector{_Violation},
sc::UnitCommitmentScenario,
)::Nothing
for v in violations
_enforce_transmission(
model = model,
sc = sc,
violation = v,
isf = model[:isf],
lodf = model[:lodf],
isf = sc.isf,
lodf = sc.lodf,
)
end
return
@ -19,6 +21,7 @@ end
function _enforce_transmission(;
model::JuMP.Model,
sc::UnitCommitmentScenario,
violation::_Violation,
isf::Matrix{Float64},
lodf::Matrix{Float64},
@ -31,19 +34,21 @@ function _enforce_transmission(;
if violation.outage_line === nothing
limit = violation.monitored_line.normal_flow_limit[violation.time]
@info @sprintf(
" %8.3f MW overflow in %-5s time %3d (pre-contingency)",
" %8.3f MW overflow in %-5s time %3d (pre-contingency, scenario %s)",
violation.amount,
violation.monitored_line.name,
violation.time,
sc.name,
)
else
limit = violation.monitored_line.emergency_flow_limit[violation.time]
@info @sprintf(
" %8.3f MW overflow in %-5s time %3d (outage: line %s)",
" %8.3f MW overflow in %-5s time %3d (outage: line %s, scenario %s)",
violation.amount,
violation.monitored_line.name,
violation.time,
violation.outage_line.name,
sc.name,
)
end
@ -51,7 +56,7 @@ function _enforce_transmission(;
t = violation.time
flow = @variable(model, base_name = "flow[$fm,$t]")
v = overflow[violation.monitored_line.name, violation.time]
v = overflow[sc.name, violation.monitored_line.name, violation.time]
@constraint(model, flow <= limit + v)
@constraint(model, -flow <= limit + v)
@ -59,23 +64,23 @@ function _enforce_transmission(;
@constraint(
model,
flow == sum(
net_injection[b.name, violation.time] *
net_injection[sc.name, b.name, violation.time] *
isf[violation.monitored_line.offset, b.offset] for
b in instance.buses if b.offset > 0
b in sc.buses if b.offset > 0
)
)
else
@constraint(
model,
flow == sum(
net_injection[b.name, violation.time] * (
net_injection[sc.name, b.name, violation.time] * (
isf[violation.monitored_line.offset, b.offset] + (
lodf[
violation.monitored_line.offset,
violation.outage_line.offset,
] * isf[violation.outage_line.offset, b.offset]
)
) for b in instance.buses if b.offset > 0
) for b in sc.buses if b.offset > 0
)
)
end

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

@ -5,39 +5,35 @@
import Base.Threads: @threads
function _find_violations(
model::JuMP.Model;
model::JuMP.Model,
sc::UnitCommitmentScenario;
max_per_line::Int,
max_per_period::Int,
)
instance = model[:instance]
net_injection = model[:net_injection]
overflow = model[:overflow]
length(instance.buses) > 1 || return []
length(sc.buses) > 1 || return []
violations = []
@info "Verifying transmission limits..."
time_screening = @elapsed begin
non_slack_buses = [b for b in instance.buses if b.offset > 0]
net_injection_values = [
value(net_injection[b.name, t]) for b in non_slack_buses,
t in 1:instance.time
]
overflow_values = [
value(overflow[lm.name, t]) for lm in instance.lines,
t in 1:instance.time
]
violations = UnitCommitment._find_violations(
instance = instance,
net_injections = net_injection_values,
overflow = overflow_values,
isf = model[:isf],
lodf = model[:lodf],
max_per_line = max_per_line,
max_per_period = max_per_period,
)
end
@info @sprintf(
"Verified transmission limits in %.2f seconds",
time_screening
non_slack_buses = [b for b in sc.buses if b.offset > 0]
net_injection_values = [
value(net_injection[sc.name, b.name, t]) for b in non_slack_buses,
t in 1:instance.time
]
overflow_values = [
value(overflow[sc.name, lm.name, t]) for lm in sc.lines,
t in 1:instance.time
]
violations = UnitCommitment._find_violations(
instance = instance,
sc = sc,
net_injections = net_injection_values,
overflow = overflow_values,
isf = sc.isf,
lodf = sc.lodf,
max_per_line = max_per_line,
max_per_period = max_per_period,
)
return violations
end
@ -64,6 +60,7 @@ matrix, where L is the number of transmission lines.
"""
function _find_violations(;
instance::UnitCommitmentInstance,
sc::UnitCommitmentScenario,
net_injections::Array{Float64,2},
overflow::Array{Float64,2},
isf::Array{Float64,2},
@ -71,8 +68,8 @@ function _find_violations(;
max_per_line::Int,
max_per_period::Int,
)::Array{_Violation,1}
B = length(instance.buses) - 1
L = length(instance.lines)
B = length(sc.buses) - 1
L = length(sc.lines)
T = instance.time
K = nthreads()
@ -93,17 +90,17 @@ function _find_violations(;
post_v::Array{Float64} = zeros(L, L, K) # post_v[lm, lc, thread]
normal_limits::Array{Float64,2} = [
l.normal_flow_limit[t] + overflow[l.offset, t] for
l in instance.lines, t in 1:T
l.normal_flow_limit[t] + overflow[l.offset, t] for l in sc.lines,
t in 1:T
]
emergency_limits::Array{Float64,2} = [
l.emergency_flow_limit[t] + overflow[l.offset, t] for
l in instance.lines, t in 1:T
l.emergency_flow_limit[t] + overflow[l.offset, t] for l in sc.lines,
t in 1:T
]
is_vulnerable::Array{Bool} = zeros(Bool, L)
for c in instance.contingencies
for c in sc.contingencies
is_vulnerable[c.lines[1].offset] = true
end
@ -144,7 +141,7 @@ function _find_violations(;
filters[t],
_Violation(
time = t,
monitored_line = instance.lines[lm],
monitored_line = sc.lines[lm],
outage_line = nothing,
amount = pre_v[lm, k],
),
@ -159,8 +156,8 @@ function _find_violations(;
filters[t],
_Violation(
time = t,
monitored_line = instance.lines[lm],
outage_line = instance.lines[lc],
monitored_line = sc.lines[lm],
outage_line = sc.lines[lc],
amount = post_v[lm, lc, k],
),
)

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

@ -12,10 +12,15 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing
end
initial_time = time()
large_gap = false
has_transmission = (length(model[:isf]) > 0)
if has_transmission && method.two_phase_gap
set_gap(1e-2)
large_gap = true
has_transmission = false
for sc in model[:instance].scenarios
if length(sc.isf) > 0
has_transmission = true
end
if has_transmission && method.two_phase_gap
set_gap(1e-2)
large_gap = true
end
end
while true
time_elapsed = time() - initial_time
@ -31,13 +36,41 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing
JuMP.set_time_limit_sec(model, time_remaining)
@info "Solving MILP..."
JuMP.optimize!(model)
has_transmission || break
violations = _find_violations(
model,
max_per_line = method.max_violations_per_line,
max_per_period = method.max_violations_per_period,
@info "Verifying transmission limits..."
time_screening = @elapsed begin
violations = []
for sc in model[:instance].scenarios
push!(
violations,
_find_violations(
model,
sc,
max_per_line = method.max_violations_per_line,
max_per_period = method.max_violations_per_period,
),
)
end
end
@info @sprintf(
"Verified transmission limits in %.2f seconds",
time_screening
)
if isempty(violations)
violations_found = false
for v in violations
if !isempty(v)
violations_found = true
end
end
if violations_found
for (i, v) in enumerate(violations)
_enforce_transmission(model, v, model[:instance].scenarios[i])
end
else
@info "No violations found"
if large_gap
large_gap = false
@ -45,8 +78,6 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing
else
break
end
else
_enforce_transmission(model, violations)
end
end
return

154
src/solution/solution.jl
Unescape View File

@ -16,34 +16,44 @@ solution = UnitCommitment.solution(model)
"""
function solution(model::JuMP.Model)::OrderedDict
instance, T = model[:instance], model[:instance].time
function timeseries(vars, collection)
return OrderedDict(
b.name => [round(value(vars[b.name, t]), digits = 5) for t in 1:T]
for b in collection
)
function timeseries(vars, collection; sc = nothing)
if sc === nothing
return OrderedDict(
b.name =>
[round(value(vars[b.name, t]), digits = 5) for t in 1:T] for
b in collection
)
else
return OrderedDict(
b.name => [
round(value(vars[sc.name, b.name, t]), digits = 5) for
t in 1:T
] for b in collection
)
end
end
function production_cost(g)
function production_cost(g, sc)
return [
value(model[:is_on][g.name, t]) * g.min_power_cost[t] + sum(
Float64[
value(model[:segprod][g.name, t, k]) *
value(model[:segprod][sc.name, g.name, t, k]) *
g.cost_segments[k].cost[t] for
k in 1:length(g.cost_segments)
],
) for t in 1:T
]
end
function production(g)
function production(g, sc)
return [
value(model[:is_on][g.name, t]) * g.min_power[t] + sum(
Float64[
value(model[:segprod][g.name, t, k]) for
value(model[:segprod][sc.name, g.name, t, k]) for
k in 1:length(g.cost_segments)
],
) for t in 1:T
]
end
function startup_cost(g)
function startup_cost(g, sc)
S = length(g.startup_categories)
return [
sum(
@ -53,66 +63,70 @@ function solution(model::JuMP.Model)::OrderedDict
]
end
sol = OrderedDict()
sol["Production (MW)"] =
OrderedDict(g.name => production(g) for g in instance.units)
sol["Production cost (\$)"] =
OrderedDict(g.name => production_cost(g) for g in instance.units)
sol["Startup cost (\$)"] =
OrderedDict(g.name => startup_cost(g) for g in instance.units)
sol["Is on"] = timeseries(model[:is_on], instance.units)
sol["Switch on"] = timeseries(model[:switch_on], instance.units)
sol["Switch off"] = timeseries(model[:switch_off], instance.units)
sol["Net injection (MW)"] =
timeseries(model[:net_injection], instance.buses)
sol["Load curtail (MW)"] = timeseries(model[:curtail], instance.buses)
if !isempty(instance.lines)
sol["Line overflow (MW)"] = timeseries(model[:overflow], instance.lines)
for sc in instance.scenarios
sol[sc.name] = OrderedDict()
sol[sc.name]["Production (MW)"] =
OrderedDict(g.name => production(g, sc) for g in sc.units)
sol[sc.name]["Production cost (\$)"] =
OrderedDict(g.name => production_cost(g, sc) for g in sc.units)
sol[sc.name]["Startup cost (\$)"] =
OrderedDict(g.name => startup_cost(g, sc) for g in sc.units)
sol[sc.name]["Is on"] = timeseries(model[:is_on], sc.units)
sol[sc.name]["Switch on"] = timeseries(model[:switch_on], sc.units)
sol[sc.name]["Switch off"] = timeseries(model[:switch_off], sc.units)
sol[sc.name]["Net injection (MW)"] =
timeseries(model[:net_injection], sc.buses, sc = sc)
sol[sc.name]["Load curtail (MW)"] =
timeseries(model[:curtail], sc.buses, sc = sc)
if !isempty(sc.lines)
sol[sc.name]["Line overflow (MW)"] =
timeseries(model[:overflow], sc.lines, sc = sc)
end
if !isempty(sc.price_sensitive_loads)
sol[sc.name]["Price-sensitive loads (MW)"] =
timeseries(model[:loads], sc.price_sensitive_loads, sc = sc)
end
sol[sc.name]["Spinning reserve (MW)"] = OrderedDict(
r.name => OrderedDict(
g.name => [
value(model[:reserve][sc.name, r.name, g.name, t]) for t in 1:instance.time
] for g in r.units
) for r in sc.reserves if r.type == "spinning"
)
sol[sc.name]["Spinning reserve shortfall (MW)"] = OrderedDict(
r.name => [
value(model[:reserve_shortfall][sc.name, r.name, t]) for
t in 1:instance.time
] for r in sc.reserves if r.type == "spinning"
)
sol[sc.name]["Up-flexiramp (MW)"] = OrderedDict(
r.name => OrderedDict(
g.name => [
value(model[:upflexiramp][sc.name, r.name, g.name, t]) for t in 1:instance.time
] for g in r.units
) for r in sc.reserves if r.type == "flexiramp"
)
sol[sc.name]["Up-flexiramp shortfall (MW)"] = OrderedDict(
r.name => [
value(model[:upflexiramp_shortfall][sc.name, r.name, t]) for t in 1:instance.time
] for r in sc.reserves if r.type == "flexiramp"
)
sol[sc.name]["Down-flexiramp (MW)"] = OrderedDict(
r.name => OrderedDict(
g.name => [
value(model[:dwflexiramp][sc.name, r.name, g.name, t]) for t in 1:instance.time
] for g in r.units
) for r in sc.reserves if r.type == "flexiramp"
)
sol[sc.name]["Down-flexiramp shortfall (MW)"] = OrderedDict(
r.name => [
value(model[:dwflexiramp_shortfall][sc.name, r.name, t]) for t in 1:instance.time
] for r in sc.reserves if r.type == "flexiramp"
)
end
if !isempty(instance.price_sensitive_loads)
sol["Price-sensitive loads (MW)"] =
timeseries(model[:loads], instance.price_sensitive_loads)
if length(instance.scenarios) == 1
return first(values(sol))
else
return sol
end
sol["Spinning reserve (MW)"] = OrderedDict(
r.name => OrderedDict(
g.name => [
value(model[:reserve][r.name, g.name, t]) for
t in 1:instance.time
] for g in r.units
) for r in instance.reserves if r.type == "spinning"
)
sol["Spinning reserve shortfall (MW)"] = OrderedDict(
r.name => [
value(model[:reserve_shortfall][r.name, t]) for
t in 1:instance.time
] for r in instance.reserves if r.type == "spinning"
)
sol["Up-flexiramp (MW)"] = OrderedDict(
r.name => OrderedDict(
g.name => [
value(model[:upflexiramp][r.name, g.name, t]) for
t in 1:instance.time
] for g in r.units
) for r in instance.reserves if r.type == "flexiramp"
)
sol["Up-flexiramp shortfall (MW)"] = OrderedDict(
r.name => [
value(model[:upflexiramp_shortfall][r.name, t]) for
t in 1:instance.time
] for r in instance.reserves if r.type == "flexiramp"
)
sol["Down-flexiramp (MW)"] = OrderedDict(
r.name => OrderedDict(
g.name => [
value(model[:dwflexiramp][r.name, g.name, t]) for
t in 1:instance.time
] for g in r.units
) for r in instance.reserves if r.type == "flexiramp"
)
sol["Down-flexiramp shortfall (MW)"] = OrderedDict(
r.name => [
value(model[:upflexiramp_shortfall][r.name, t]) for
t in 1:instance.time
] for r in instance.reserves if r.type == "flexiramp"
)
return sol
end

12
src/transform/initcond.jl
Unescape View File

@ -5,18 +5,18 @@
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
optimizer. The instance is modified in-place.
optimizer. The scenario is modified in-place.
"""
function generate_initial_conditions!(
instance::UnitCommitmentInstance,
sc::UnitCommitmentScenario,
optimizer,
)::Nothing
G = instance.units
B = instance.buses
G = sc.units
B = sc.buses
t = 1
mip = JuMP.Model(optimizer)

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

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

50
src/transform/slice.jl
Unescape View File

@ -24,31 +24,33 @@ function slice(
)::UnitCommitmentInstance
modified = deepcopy(instance)
modified.time = length(range)
modified.power_balance_penalty = modified.power_balance_penalty[range]
for r in modified.reserves
r.amount = r.amount[range]
end
for u in modified.units
u.max_power = u.max_power[range]
u.min_power = u.min_power[range]
u.must_run = u.must_run[range]
u.min_power_cost = u.min_power_cost[range]
for s in u.cost_segments
s.mw = s.mw[range]
s.cost = s.cost[range]
for sc in modified.scenarios
sc.power_balance_penalty = sc.power_balance_penalty[range]
for r in sc.reserves
r.amount = r.amount[range]
end
for u in sc.units
u.max_power = u.max_power[range]
u.min_power = u.min_power[range]
u.must_run = u.must_run[range]
u.min_power_cost = u.min_power_cost[range]
for s in u.cost_segments
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
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
return modified
end

10
src/validation/repair.jl
Unescape View File

@ -3,19 +3,19 @@
# Released under the modified BSD license. See COPYING.md for more details.
"""
repair!(instance)
repair!(sc)
Verifies that the given unit commitment instance is valid and automatically
Verifies that the given unit commitment scenario is valid and automatically
fixes some validation errors if possible, issuing a warning for each error
found. If a validation error cannot be automatically fixed, issues an
exception.
Returns the number of validation errors found.
"""
function repair!(instance::UnitCommitmentInstance)::Int
function repair!(sc::UnitCommitmentScenario)::Int
n_errors = 0
for g in instance.units
for g in sc.units
# Startup costs and delays must be increasing
for s in 2:length(g.startup_categories)
@ -38,7 +38,7 @@ function repair!(instance::UnitCommitmentInstance)::Int
end
end
for t in 1:instance.time
for t in 1:sc.time
# Production cost curve should be convex
for k in 2:length(g.cost_segments)
cost = g.cost_segments[k].cost[t]

609
src/validation/validate.jl
Unescape View File

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

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
instance = UnitCommitment.read("$FIXTURES/ucjl-0.2.json.gz")
@test length(instance.reserves_by_name["r1"].amount) == 4
@test instance.units_by_name["g2"].reserves[1].name == "r1"
@test length(instance.scenarios) == 1
sc = instance.scenarios[1]
@test length(sc.reserves_by_name["r1"].amount) == 4
@test sc.units_by_name["g2"].reserves[1].name == "r1"
end
@testset "read v0.3" begin
instance = UnitCommitment.read("$FIXTURES/ucjl-0.3.json.gz")
@test length(instance.units) == 6
@test length(instance.buses) == 14
@test length(instance.lines) == 20
@test length(instance.scenarios) == 1
sc = instance.scenarios[1]
@test length(sc.units) == 6
@test length(sc.buses) == 14
@test length(sc.lines) == 20
end

87
test/instance/read_test.jl
Unescape View File

@ -7,42 +7,49 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
@testset "read_benchmark" begin
instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
@test length(instance.lines) == 20
@test length(instance.buses) == 14
@test length(instance.units) == 6
@test length(instance.contingencies) == 19
@test length(instance.price_sensitive_loads) == 1
@test repr(instance) == (
"UnitCommitmentInstance(1 scenarios, 6 units, 14 buses, " *
"20 lines, 19 contingencies, 1 price sensitive loads, 4 time steps)"
)
@test length(instance.scenarios) == 1
sc = instance.scenarios[1]
@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 instance.lines[5].name == "l5"
@test instance.lines[5].source.name == "b2"
@test instance.lines[5].target.name == "b5"
@test instance.lines[5].reactance ≈ 0.17388
@test instance.lines[5].susceptance ≈ 10.037550333
@test instance.lines[5].normal_flow_limit == [1e8 for t in 1:4]
@test instance.lines[5].emergency_flow_limit == [1e8 for t in 1:4]
@test instance.lines[5].flow_limit_penalty == [5e3 for t in 1:4]
@test instance.lines_by_name["l5"].name == "l5"
@test sc.lines[5].name == "l5"
@test sc.lines[5].source.name == "b2"
@test sc.lines[5].target.name == "b5"
@test sc.lines[5].reactance ≈ 0.17388
@test sc.lines[5].susceptance ≈ 10.037550333
@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[1].name == "l1"
@test instance.lines[1].source.name == "b1"
@test instance.lines[1].target.name == "b2"
@test instance.lines[1].reactance ≈ 0.059170
@test instance.lines[1].susceptance ≈ 29.496860773945
@test instance.lines[1].normal_flow_limit == [300.0 for t in 1:4]
@test instance.lines[1].emergency_flow_limit == [400.0 for t in 1:4]
@test instance.lines[1].flow_limit_penalty == [1e3 for t in 1:4]
@test sc.lines[1].name == "l1"
@test sc.lines[1].source.name == "b1"
@test sc.lines[1].target.name == "b2"
@test sc.lines[1].reactance ≈ 0.059170
@test sc.lines[1].susceptance ≈ 29.496860773945
@test sc.lines[1].normal_flow_limit == [300.0 for t in 1:4]
@test sc.lines[1].emergency_flow_limit == [400.0 for t in 1:4]
@test sc.lines[1].flow_limit_penalty == [1e3 for t in 1:4]
@test instance.buses[9].name == "b9"
@test instance.buses[9].load == [35.36638, 33.25495, 31.67138, 31.14353]
@test instance.buses_by_name["b9"].name == "b9"
@test sc.buses[9].name == "b9"
@test sc.buses[9].load == [35.36638, 33.25495, 31.67138, 31.14353]
@test sc.buses_by_name["b9"].name == "b9"
@test instance.reserves[1].name == "r1"
@test instance.reserves[1].type == "spinning"
@test instance.reserves[1].amount == [100.0, 100.0, 100.0, 100.0]
@test instance.reserves_by_name["r1"].name == "r1"
@test sc.reserves[1].name == "r1"
@test sc.reserves[1].type == "spinning"
@test sc.reserves[1].amount == [100.0, 100.0, 100.0, 100.0]
@test sc.reserves_by_name["r1"].name == "r1"
unit = instance.units[1]
unit = sc.units[1]
@test unit.name == "g1"
@test unit.bus.name == "b1"
@test unit.ramp_up_limit == 1e6
@ -69,14 +76,14 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
@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"
@test sc.units_by_name["g1"].name == "g1"
unit = instance.units[2]
unit = sc.units[2]
@test unit.name == "g2"
@test unit.must_run == [false for t in 1:4]
@test length(unit.reserves) == 1
unit = instance.units[3]
unit = sc.units[3]
@test unit.name == "g3"
@test unit.bus.name == "b3"
@test unit.ramp_up_limit == 70.0
@ -98,23 +105,23 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
@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"
@test sc.contingencies[1].lines == [sc.lines[1]]
@test sc.contingencies[1].units == []
@test sc.contingencies[1].name == "c1"
@test sc.contingencies_by_name["c1"].name == "c1"
load = instance.price_sensitive_loads[1]
load = sc.price_sensitive_loads[1]
@test load.name == "ps1"
@test load.bus.name == "b3"
@test load.revenue == [100.0 for t in 1:4]
@test load.demand == [50.0 for t in 1:4]
@test instance.price_sensitive_loads_by_name["ps1"].name == "ps1"
@test sc.price_sensitive_loads_by_name["ps1"].name == "ps1"
end
@testset "read_benchmark sub-hourly" begin
instance = UnitCommitment.read("$FIXTURES/case14-sub-hourly.json.gz")
@test instance.time == 4
unit = instance.units[1]
unit = instance.scenarios[1].units[1]
@test unit.name == "g1"
@test unit.min_uptime == 2
@test unit.min_downtime == 2

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

@ -7,13 +7,13 @@ import UnitCommitment: _Violation, _offer, _query
@testset "_ViolationFilter" begin
instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
sc = instance.scenarios[1]
filter = UnitCommitment._ViolationFilter(max_per_line = 1, max_total = 2)
_offer(
filter,
_Violation(
time = 1,
monitored_line = instance.lines[1],
monitored_line = sc.lines[1],
outage_line = nothing,
amount = 100.0,
),
@ -22,8 +22,8 @@ import UnitCommitment: _Violation, _offer, _query
filter,
_Violation(
time = 1,
monitored_line = instance.lines[1],
outage_line = instance.lines[1],
monitored_line = sc.lines[1],
outage_line = sc.lines[1],
amount = 300.0,
),
)
@ -31,8 +31,8 @@ import UnitCommitment: _Violation, _offer, _query
filter,
_Violation(
time = 1,
monitored_line = instance.lines[1],
outage_line = instance.lines[5],
monitored_line = sc.lines[1],
outage_line = sc.lines[5],
amount = 500.0,
),
)
@ -40,8 +40,8 @@ import UnitCommitment: _Violation, _offer, _query
filter,
_Violation(
time = 1,
monitored_line = instance.lines[1],
outage_line = instance.lines[4],
monitored_line = sc.lines[1],
outage_line = sc.lines[4],
amount = 400.0,
),
)
@ -49,8 +49,8 @@ import UnitCommitment: _Violation, _offer, _query
filter,
_Violation(
time = 1,
monitored_line = instance.lines[2],
outage_line = instance.lines[1],
monitored_line = sc.lines[2],
outage_line = sc.lines[1],
amount = 200.0,
),
)
@ -58,8 +58,8 @@ import UnitCommitment: _Violation, _offer, _query
filter,
_Violation(
time = 1,
monitored_line = instance.lines[2],
outage_line = instance.lines[8],
monitored_line = sc.lines[2],
outage_line = sc.lines[8],
amount = 100.0,
),
)
@ -68,14 +68,14 @@ import UnitCommitment: _Violation, _offer, _query
expected = [
_Violation(
time = 1,
monitored_line = instance.lines[2],
outage_line = instance.lines[1],
monitored_line = sc.lines[2],
outage_line = sc.lines[1],
amount = 200.0,
),
_Violation(
time = 1,
monitored_line = instance.lines[1],
outage_line = instance.lines[5],
monitored_line = sc.lines[1],
outage_line = sc.lines[5],
amount = 500.0,
),
]

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

@ -7,23 +7,25 @@ import UnitCommitment: _Violation, _offer, _query
@testset "find_violations" begin
instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
for line in instance.lines, t in 1:instance.time
sc = instance.scenarios[1]
for line in sc.lines, t in 1:instance.time
line.normal_flow_limit[t] = 1.0
line.emergency_flow_limit[t] = 1.0
end
isf = UnitCommitment._injection_shift_factors(
lines = instance.lines,
buses = instance.buses,
lines = sc.lines,
buses = sc.buses,
)
lodf = UnitCommitment._line_outage_factors(
lines = instance.lines,
buses = instance.buses,
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 instance.lines, 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,

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

@ -6,7 +6,8 @@ using UnitCommitment, Test, LinearAlgebra
@testset "_susceptance_matrix" begin
instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
actual = UnitCommitment._susceptance_matrix(instance.lines)
sc = instance.scenarios[1]
actual = UnitCommitment._susceptance_matrix(sc.lines)
@test size(actual) == (20, 20)
expected = Diagonal([
29.5,
@ -35,9 +36,10 @@ end
@testset "_reduced_incidence_matrix" begin
instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
sc = instance.scenarios[1]
actual = UnitCommitment._reduced_incidence_matrix(
lines = instance.lines,
buses = instance.buses,
lines = sc.lines,
buses = sc.buses,
)
@test size(actual) == (20, 13)
@test actual[1, 1] == -1.0
@ -82,9 +84,10 @@ end
@testset "_injection_shift_factors" begin
instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
sc = instance.scenarios[1]
actual = UnitCommitment._injection_shift_factors(
lines = instance.lines,
buses = instance.buses,
lines = sc.lines,
buses = sc.buses,
)
@test size(actual) == (20, 13)
@test round.(actual, digits = 2) == [
@ -113,25 +116,26 @@ end
@testset "_line_outage_factors" begin
instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
sc = instance.scenarios[1]
isf_before = UnitCommitment._injection_shift_factors(
lines = instance.lines,
buses = instance.buses,
lines = sc.lines,
buses = sc.buses,
)
lodf = UnitCommitment._line_outage_factors(
lines = instance.lines,
buses = instance.buses,
lines = sc.lines,
buses = sc.buses,
isf = isf_before,
)
for contingency in instance.contingencies
for contingency in sc.contingencies
for lc in contingency.lines
prev_susceptance = lc.susceptance
lc.susceptance = 0.0
isf_after = UnitCommitment._injection_shift_factors(
lines = instance.lines,
buses = instance.buses,
lines = sc.lines,
buses = sc.buses,
)
lc.susceptance = prev_susceptance
for lm in instance.lines
for lm in sc.lines
expected = isf_after[lm.offset, :]
actual =
isf_before[lm.offset, :] +

8
test/transform/initcond_test.jl
Unescape View File

@ -8,18 +8,18 @@ using UnitCommitment, Cbc, JuMP
# Load instance
instance = UnitCommitment.read("$FIXTURES/case118-initcond.json.gz")
optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
sc = instance.scenarios[1]
# All units should have unknown initial conditions
for g in instance.units
for g in sc.units
@test g.initial_power === nothing
@test g.initial_status === nothing
end
# Generate initial conditions
UnitCommitment.generate_initial_conditions!(instance, optimizer)
UnitCommitment.generate_initial_conditions!(sc, optimizer)
# 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_status !== nothing
end

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

@ -9,30 +9,31 @@ using Distributions
using Random
using UnitCommitment, Cbc, JuMP
get_instance() = UnitCommitment.read_benchmark("matpower/case118/2017-02-01")
system_load(instance) = sum(b.load for b in instance.buses)
function get_scenario()
return UnitCommitment.read_benchmark(
"matpower/case118/2017-02-01",
).scenarios[1]
end
system_load(sc) = sum(b.load for b in sc.buses)
test_approx(x, y) = @test isapprox(x, y, atol = 1e-3)
@testset "XavQiuAhm2021" begin
@testset "cost and load share" begin
instance = get_instance()
sc = get_scenario()
# Check original costs
unit = instance.units[10]
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 load share
bus = instance.buses[1]
prev_system_load = system_load(instance)
bus = sc.buses[1]
prev_system_load = system_load(sc)
test_approx(bus.load[1] / prev_system_load[1], 0.012)
randomize!(
instance,
method = XavQiuAhm2021.Randomization(
randomize_load_profile = false,
),
sc,
XavQiuAhm2021.Randomization(randomize_load_profile = false),
rng = MersenneTwister(42),
)
@ -42,7 +43,7 @@ test_approx(x, y) = @test isapprox(x, y, atol = 1e-3)
test_approx(unit.startup_categories[1].cost[1], 7634.226)
# Check randomized load share
curr_system_load = system_load(instance)
curr_system_load = system_load(sc)
test_approx(bus.load[1] / curr_system_load[1], 0.013)
# System load should not change
@ -50,20 +51,15 @@ test_approx(x, y) = @test isapprox(x, y, atol = 1e-3)
end
@testset "load profile" begin
instance = get_instance()
sc = get_scenario()
# Check original load profile
@test round.(system_load(instance), digits = 1)[1:8] ≈
@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]
randomize!(
instance,
XavQiuAhm2021.Randomization(),
rng = MersenneTwister(42),
)
randomize!(sc, XavQiuAhm2021.Randomization(); rng = MersenneTwister(42))
# 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]
end
end

14
test/transform/slice_test.jl
Unescape View File

@ -7,12 +7,13 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
@testset "slice" begin
instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
modified = UnitCommitment.slice(instance, 1:2)
sc = modified.scenarios[1]
# Should update all time-dependent fields
@test modified.time == 2
@test length(modified.power_balance_penalty) == 2
@test length(modified.reserves_by_name["r1"].amount) == 2
for u in modified.units
@test length(sc.power_balance_penalty) == 2
@test length(sc.reserves_by_name["r1"].amount) == 2
for u in sc.units
@test length(u.max_power) == 2
@test length(u.min_power) == 2
@test length(u.must_run) == 2
@ -22,18 +23,19 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip
@test length(s.cost) == 2
end
end
for b in modified.buses
for b in sc.buses
@test length(b.load) == 2
end
for l in modified.lines
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 modified.price_sensitive_loads
for ps in sc.price_sensitive_loads
@test length(ps.demand) == 2
@test length(ps.revenue) == 2
end
# Should be able to build model without errors
optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
model = UnitCommitment.build_model(

81
test/usage.jl
Unescape View File

@ -4,34 +4,59 @@
using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON
function _set_flow_limits!(instance)
for sc in instance.scenarios
sc.power_balance_penalty = [100_000 for _ in 1:instance.time]
for line in sc.lines, t in 1:4
line.normal_flow_limit[t] = 10.0
end
end
end
@testset "usage" begin
instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
for line in instance.lines, t in 1:4
line.normal_flow_limit[t] = 10.0
@testset "deterministic" begin
instance = UnitCommitment.read("$FIXTURES/case14.json.gz")
_set_flow_limits!(instance)
optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
model = UnitCommitment.build_model(
instance = instance,
optimizer = optimizer,
variable_names = true,
)
@test name(model[:is_on]["g1", 1]) == "is_on[g1,1]"
# Optimize and retrieve solution
UnitCommitment.optimize!(model)
solution = UnitCommitment.solution(model)
# Write solution to a file
filename = tempname()
UnitCommitment.write(filename, solution)
loaded = JSON.parsefile(filename)
@test length(loaded["Is on"]) == 6
# Verify solution
@test UnitCommitment.validate(instance, solution)
# Reoptimize with fixed solution
UnitCommitment.fix!(model, solution)
UnitCommitment.optimize!(model)
@test UnitCommitment.validate(instance, solution)
end
@testset "stochastic" begin
instance = UnitCommitment.read([
"$FIXTURES/case14.json.gz",
"$FIXTURES/case14.json.gz",
])
_set_flow_limits!(instance)
@test length(instance.scenarios) == 2
optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
model = UnitCommitment.build_model(
instance = instance,
optimizer = optimizer,
)
UnitCommitment.optimize!(model)
solution = UnitCommitment.solution(model)
end
optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
model = UnitCommitment.build_model(
instance = instance,
optimizer = optimizer,
variable_names = true,
)
@test name(model[:is_on]["g1", 1]) == "is_on[g1,1]"
# Optimize and retrieve solution
UnitCommitment.optimize!(model)
solution = UnitCommitment.solution(model)
# Write solution to a file
filename = tempname()
UnitCommitment.write(filename, solution)
loaded = JSON.parsefile(filename)
@test length(loaded["Is on"]) == 6
# Verify solution
@test UnitCommitment.validate(instance, solution)
# Reoptimize with fixed solution
UnitCommitment.fix!(model, solution)
UnitCommitment.optimize!(model)
@test UnitCommitment.validate(instance, solution)
end

12
test/validation/repair_test.jl
Unescape View File

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

Write Preview
Loading…
Cancel
Save