ANL-CEEESA/UnitCommitment.jl
1
0
Fork 0
mirror of https://github.com/ANL-CEEESA/UnitCommitment.jl.git synced 2025-12-06 00:08:52 -06:00
Code Issues Packages Projects Releases Wiki Activity

stochastic extension w/ scenarios

Browse Source
This commit is contained in:
oyurdakul
2023-02-08 23:46:10 -06:00
parent 8fc84412eb
commit c95b01dadf
27 changed files with 566 additions and 386 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
Project.toml
View File

@@ -9,6 +9,7 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63"
Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

1
src/UnitCommitment.jl
View File

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

72
src/instance/read.jl
View File

@@ -7,6 +7,7 @@ using JSON
using DataStructures
using GZip
import Base: getindex, time
using Glob
const INSTANCES_URL = "https://axavier.org/UnitCommitment.jl/0.3/instances"
@@ -43,6 +44,25 @@ function read_benchmark(
return UnitCommitment.read(filename)
end
function read_scenarios(path::AbstractString)::UnitCommitmentInstance
scenario_paths = glob("*.json", path)
scenarios = Vector{UnitCommitmentScenario}()
total_number_of_scenarios = length(scenario_paths)
for (scenario_index, scenario_path) in enumerate(scenario_paths)
scenario = read_scenario(scenario_path,
total_number_of_scenarios = total_number_of_scenarios,
scenario_index = scenario_index)
push!(scenarios, scenario)
end
instance = UnitCommitmentInstance(
time = scenarios[1].time,
scenarios = scenarios
)
abs(sum(scenario.probability for scenario in instance.scenarios) - 1.0) <= 0.01 ||
error("scenario probabilities do not add up to one")
return instance
end
"""
read(path::AbstractString)::UnitCommitmentInstance
@@ -54,18 +74,33 @@ Read instance from a file. The file may be gzipped.
instance = UnitCommitment.read("/path/to/input.json.gz")
```
"""
function read(path::AbstractString)::UnitCommitmentInstance
function read_scenario(path::AbstractString; total_number_of_scenarios = 1, scenario_index = 1)::UnitCommitmentScenario
if endswith(path, ".gz")
return _read(gzopen(path))
return _read(gzopen(path),total_number_of_scenarios, scenario_index)
else
return _read(open(path))
return _read(open(path), total_number_of_scenarios, scenario_index)
end
end
function _read(file::IO)::UnitCommitmentInstance
function read(path::AbstractString; total_number_of_scenarios = 1, scenario_index = 1)::UnitCommitmentInstance
if endswith(path, ".gz")
scenario = _read(gzopen(path),total_number_of_scenarios, scenario_index)
else
scenario = _read(open(path), total_number_of_scenarios, scenario_index)
end
instance = UnitCommitmentInstance(
time = scenario.time,
scenarios = [scenario]
)
return instance
end
function _read(file::IO, total_number_of_scenarios::Int, scenario_index::Int)::UnitCommitmentScenario
return _from_json(
JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)),
)
total_number_of_scenarios,
scenario_index
)
end
function _read_json(path::String)::OrderedDict
@@ -77,7 +112,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, total_number_of_scenarios::Int, scenario_index::Int; repair = true)::UnitCommitmentScenario
_migrate(json)
units = Unit[]
buses = Bus[]
@@ -101,6 +136,17 @@ function _from_json(json; repair = true)
error("Time step $time_step is not a divisor of 60")
time_multiplier = 60 ÷ time_step
T = time_horizon * time_multiplier
#####
probability = json["Parameters"]["Scenario probability"]
if probability === nothing
probability = (1 / total_number_of_scenarios)
end
scenario_name = json["Parameters"]["Scenario name"]
if scenario_name === nothing
scenario_name = "s$(scenario_index)"
end
######
name_to_bus = Dict{String,Bus}()
name_to_line = Dict{String,TransmissionLine}()
@@ -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,16 @@ 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,
# 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

37
src/instance/structs.jl
View File

@@ -73,7 +73,9 @@ mutable struct PriceSensitiveLoad
revenue::Vector{Float64}
end
Base.@kwdef mutable struct UnitCommitmentInstance
Base.@kwdef mutable struct UnitCommitmentScenario
name::String
probability::Float64
buses_by_name::Dict{AbstractString,Bus}
buses::Vector{Bus}
contingencies_by_name::Dict{AbstractString,Contingency}
@@ -85,23 +87,34 @@ Base.@kwdef mutable struct UnitCommitmentInstance
price_sensitive_loads::Vector{PriceSensitiveLoad}
reserves::Vector{Reserve}
reserves_by_name::Dict{AbstractString,Reserve}
shortfall_penalty::Vector{Float64}
flexiramp_shortfall_penalty::Vector{Float64}
time::Int
# shortfall_penalty::Vector{Float64}
# flexiramp_shortfall_penalty::Vector{Float64}
units_by_name::Dict{AbstractString,Unit}
units::Vector{Unit}
time::Int
isf::Array{Float64,2}
lodf::Array{Float64,2}
end
Base.@kwdef mutable struct UnitCommitmentInstance
time::Int
scenarios::Vector{UnitCommitmentScenario}
end
function Base.show(io::IO, instance::UnitCommitmentInstance)
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: ")
for scenario in instance.scenarios
print(io, "Scenario $(scenario.name): ")
print(io, "$(length(scenario.units)) units, ")
print(io, "$(length(scenario.buses)) buses, ")
print(io, "$(length(scenario.lines)) lines, ")
print(io, "$(length(scenario.contingencies)) contingencies, ")
print(
io,
"$(length(scenario.price_sensitive_loads)) price sensitive loads, ",
)
end
print(io, "$(instance.time) time steps")
print(io, ")")
return

32
src/model/build.jl
View File

@@ -77,20 +77,28 @@ 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_first_stage!(model, g, formulation)
end
for b in instance.buses
_add_bus!(model, b)
for scenario in instance.scenarios
@info "Building scenario $(scenario.name) with" *
"probability $(scenario.probability)"
_setup_transmission(model, formulation.transmission, scenario)
for l in scenario.lines
_add_transmission_line!(model, l, formulation.transmission,
scenario)
end
for b in scenario.buses
_add_bus!(model, b, scenario)
end
for ps in scenario.price_sensitive_loads
_add_price_sensitive_load!(model, ps, scenario)
end
for g in scenario.units
_add_unit_second_stage!(model, g, formulation, scenario)
end
_add_system_wide_eqs!(model, scenario)
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)

23
src/model/formulations/ArrCon2000/ramp.jl
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,24 @@ 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]

15
src/model/formulations/CarArr2006/pwlcosts.jl
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,28 @@ 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

19
src/model/formulations/DamKucRajAta2016/ramp.jl
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,7 @@ 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 +77,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 +89,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 +99,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 +111,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

17
src/model/formulations/Gar1962/prod.jl
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,20 @@ 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])
### Moving this term to another function
# 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 +47,9 @@ 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

15
src/model/formulations/Gar1962/pwlcosts.jl
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,9 @@ 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 +38,21 @@ 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
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

21
src/model/formulations/KnuOstWat2018/pwlcosts.jl
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,22 @@ 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

31
src/model/formulations/MorLatRam2013/ramp.jl
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,13 @@ 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 +75,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 +91,9 @@ 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 +103,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 +119,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

15
src/model/formulations/PanGua2016/ramp.jl
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,9 +68,9 @@ 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(
eq_prod_limit_ramp_up_extra_period[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] - sum(
@@ -84,9 +85,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(

61
src/model/formulations/WanHob2016/ramp.jl
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
@@ -29,7 +30,7 @@ function _add_ramp_eqs!(
error("Each generator may only provide one flexiramp reserve")
end
for r in g.reserves
if r.type !== "flexiramp"
if r.type !== "up-frp" && r.type !== "down-frp"
error(
"This formulation only supports flexiramp reserves, not $(r.type)",
)
@@ -38,41 +39,41 @@ 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, rn, 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, rn, 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] +
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, rn, 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, rn, 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, rn, 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 +81,8 @@ 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 +90,7 @@ function _add_ramp_eqs!(
else
@constraint(
model,
mfg[rn, gn, t] <=
mfg[sc.name, rn, gn, t] <=
initial_power +
(RU * is_initially_on) +
(SU * (is_on[gn, t] - is_initially_on)) +
@@ -98,7 +99,7 @@ function _add_ramp_eqs!(
@constraint(
model,
initial_power -
(prod_above[gn, t] + (is_on[gn, t] * minp[t])) <=
(prod_above[sc.name, gn, t] + (is_on[gn, t] * minp[t])) <=
RD * is_on[gn, t] +
SD * (is_initially_on - is_on[gn, t]) +
maxp[t] * (1 - is_initially_on)
@@ -106,7 +107,7 @@ function _add_ramp_eqs!(
end
@constraint(
model,
mfg[rn, gn, t] <=
mfg[sc.name, rn, 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 +115,11 @@ 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 +127,11 @@ 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 +139,26 @@ 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, rn, 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 +166,8 @@ 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])

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

@@ -2,22 +2,22 @@
# 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

22
src/model/formulations/base/line.jl
View File

@@ -6,14 +6,15 @@ 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
@@ -22,27 +23,28 @@ end
function _setup_transmission(
model::JuMP.Model,
formulation::ShiftFactorsFormulation,
scenario::UnitCommitmentScenario
)::Nothing
instance = model[:instance]
isf = formulation.precomputed_isf
lodf = formulation.precomputed_lodf
if length(instance.buses) == 1
if length(scenario.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,
lines = scenario.lines,
buses = scenario.buses,
)
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,
lines = scenario.lines,
buses = scenario.buses,
isf = isf,
)
end
@@ -55,7 +57,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
scenario.isf = isf
scenario.lodf = lodf
return
end

10
src/model/formulations/base/psload.jl
View File

@@ -5,21 +5,23 @@
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[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
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

100
src/model/formulations/base/system.jl
View File

@@ -2,54 +2,54 @@
# 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 +57,7 @@ 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,31 +65,41 @@ 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
r.type == "flexiramp" || continue
for t in 1:instance.time
# Eq. (17) in Wang & Hobbs (2016)
eq_min_upflexiramp[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]
)
# Eq. (18) in Wang & Hobbs (2016)
eq_min_dwflexiramp[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]
)
# Account for flexiramp shortfall contribution to objective
if r.shortfall_penalty >= 0
add_to_expression!(
model[:obj],
r.shortfall_penalty,
(
model[:upflexiramp_shortfall][r.name, t] +
model[:dwflexiramp_shortfall][r.name, t]
),
T = model[:instance].time
for r in sc.reserves
if r.type == "up-frp"
for t in 1:T
# Eq. (17) in Wang & Hobbs (2016)
eq_min_upflexiramp[sc.name, r.name, t] = @constraint(
model,
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]
)
# Account for flexiramp shortfall contribution to objective
if r.shortfall_penalty >= 0
add_to_expression!(
model[:obj],
r.shortfall_penalty * sc.probability,
model[:upflexiramp_shortfall][sc.name, r.name, t]
)
end
end
elseif r.type == "down-frp"
for t in 1:T
# Eq. (18) in Wang & Hobbs (2016)
eq_min_dwflexiramp[sc.name, r.name, t] = @constraint(
model,
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 * sc.probability,
model[:dwflexiramp_shortfall][sc.name, r.name, t]
)
end
end
end
end

165
src/model/formulations/base/unit.jl
View File

@@ -2,7 +2,7 @@
# 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 _add_unit_first_stage!(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 +11,34 @@ 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 _add_unit_second_stage!(model::JuMP.Model, g::Unit, formulation::Formulation,
scenario::UnitCommitmentScenario)
# Variables
_add_production_vars!(model, g, formulation.prod_vars, scenario)
_add_spinning_reserve_vars!(model, g, scenario)
_add_flexiramp_reserve_vars!(model, g, scenario)
# Constraints and objective function
_add_net_injection_eqs!(model, g, scenario)
_add_production_limit_eqs!(model, g, formulation.prod_vars, scenario)
_add_production_piecewise_linear_eqs!(
model,
g,
formulation.prod_vars,
formulation.pwl_costs,
formulation.status_vars,
scenario
)
_add_ramp_eqs!(
model,
@@ -34,26 +46,64 @@ function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
formulation.prod_vars,
formulation.ramping,
formulation.status_vars,
scenario
)
_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, scenario)
return
end
# function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
# if !all(g.must_run) && any(g.must_run)
# error("Partially must-run units are not currently supported")
# end
# if g.initial_power === nothing || g.initial_status === nothing
# error("Initial conditions for $(g.name) must be provided")
# end
# # Variables
# _add_production_vars!(model, g, formulation.prod_vars)
# _add_spinning_reserve_vars!(model, g)
# _add_flexiramp_reserve_vars!(model, g)
# _add_startup_shutdown_vars!(model, g)
# _add_status_vars!(model, g, formulation.status_vars)
# # Constraints and objective function
# _add_min_uptime_downtime_eqs!(model, g)
# _add_net_injection_eqs!(model, g)
# _add_production_limit_eqs!(model, g, formulation.prod_vars)
# _add_production_piecewise_linear_eqs!(
# model,
# g,
# formulation.prod_vars,
# formulation.pwl_costs,
# formulation.status_vars,
# )
# _add_ramp_eqs!(
# model,
# g,
# formulation.prod_vars,
# formulation.ramping,
# formulation.status_vars,
# )
# _add_startup_cost_eqs!(model, g, formulation.startup_costs)
# _add_startup_shutdown_limit_eqs!(model, g)
# _add_status_eqs!(model, g, formulation.status_vars)
# return
# end
_is_initially_on(g::Unit)::Float64 = (g.initial_status > 0 ? 1.0 : 0.0)
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 +111,35 @@ 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] =
@variable(model, lower_bound = 0)
dwflexiramp_shortfall[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)
if r.type == "up-frp"
for t in 1:model[:instance].time
# maximum feasible generation, \bar{g_{its}} in Wang & Hobbs (2016)
mfg[sc.name, r.name, g.name, t] = @variable(model, lower_bound = 0)
upflexiramp[sc.name, r.name, g.name, t] = @variable(model) # up-flexiramp, ur_{it} in Wang & Hobbs (2016)
if (sc.name, r.name, t) ∉ keys(upflexiramp_shortfall)
upflexiramp_shortfall[sc.name, r.name, t] =
@variable(model, lower_bound = 0)
if r.shortfall_penalty < 0
set_upper_bound(upflexiramp_shortfall[sc.name, r.name, t], 0.0)
end
end
end
elseif r.type == "down-frp"
for t in 1:model[:instance].time
dwflexiramp[sc.name, r.name, g.name, t] = @variable(model) # down-flexiramp, dr_{it} in Wang & Hobbs (2016)
if (sc.name, r.name, t) ∉ keys(dwflexiramp_shortfall)
dwflexiramp_shortfall[sc.name, r.name, t] =
@variable(model, lower_bound = 0)
if r.shortfall_penalty < 0
set_upper_bound(dwflexiramp_shortfall[sc.name, r.name, t], 0.0)
end
end
end
end
@@ -99,32 +157,32 @@ 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 +196,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 +282,30 @@ 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

17
src/solution/methods/XavQiuWanThi2019/enforce.jl
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},
@@ -51,7 +54,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 +62,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

37
src/solution/methods/XavQiuWanThi2019/find.jl
View File

@@ -5,32 +5,34 @@
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]
non_slack_buses = [b for b in sc.buses if b.offset > 0]
net_injection_values = [
value(net_injection[b.name, t]) for b in non_slack_buses,
value(net_injection[sc.name, 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,
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 = model[:isf],
lodf = model[:lodf],
isf = sc.isf,
lodf = sc.lodf,
max_per_line = max_per_line,
max_per_period = max_per_period,
)
@@ -64,15 +66,16 @@ 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},
lodf::Array{Float64,2},
max_per_line::Int,
max_per_period::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()
@@ -94,16 +97,16 @@ function _find_violations(;
normal_limits::Array{Float64,2} = [
l.normal_flow_limit[t] + overflow[l.offset, t] for
l in instance.lines, t in 1:T
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 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
@@ -111,7 +114,7 @@ function _find_violations(;
k = threadid()
# Pre-contingency flows
pre_flow[:, k] = isf * net_injections[:, t]
pre_flow[:, k] = isf * net_injections[ :, t]
# Post-contingency flows
for lc in 1:L, lm in 1:L
@@ -144,7 +147,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 +162,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],
),
)

70
src/solution/methods/XavQiuWanThi2019/optimize.jl
View File

@@ -10,43 +10,47 @@ function optimize!(model::JuMP.Model, method::XavQiuWanThi2019.Method)::Nothing
JuMP.set_optimizer_attribute(model, "MIPGap", gap)
@info @sprintf("MIP gap tolerance set to %f", gap)
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
end
while true
time_elapsed = time() - initial_time
time_remaining = method.time_limit - time_elapsed
if time_remaining < 0
@info "Time limit exceeded"
break
for scenario in model[:instance].scenarios
large_gap = false
has_transmission = (length(scenario.isf) > 0)
if has_transmission && method.two_phase_gap
set_gap(1e-2)
large_gap = true
end
@info @sprintf(
"Setting MILP time limit to %.2f seconds",
time_remaining
)
JuMP.set_time_limit_sec(model, time_remaining)
@info "Solving MILP..."
JuMP.optimize!(model)
has_transmission || break
violations = _find_violations(
model,
max_per_line = method.max_violations_per_line,
max_per_period = method.max_violations_per_period,
)
if isempty(violations)
@info "No violations found"
if large_gap
large_gap = false
set_gap(method.gap_limit)
else
while true
initial_time = time()
time_elapsed = time() - initial_time
time_remaining = method.time_limit - time_elapsed
if time_remaining < 0
@info "Time limit exceeded"
break
end
else
_enforce_transmission(model, violations)
@info @sprintf(
"Setting MILP time limit to %.2f seconds",
time_remaining
)
JuMP.set_time_limit_sec(model, time_remaining)
@info "Solving MILP..."
JuMP.optimize!(model)
has_transmission || break
violations = _find_violations(
model,
scenario,
max_per_line = method.max_violations_per_line,
max_per_period = method.max_violations_per_period,
)
if isempty(violations)
@info "No violations found"
if large_gap
large_gap = false
set_gap(method.gap_limit)
else
break
end
else
_enforce_transmission(model, violations, scenario)
end
end
end
return

142
src/solution/solution.jl
View File

@@ -16,34 +16,40 @@ solution = UnitCommitment.solution(model)
"""
function solution(model::JuMP.Model)::OrderedDict
instance, T = model[:instance], model[:instance].time
function timeseries(vars, collection)
function timeseries_first_stage(vars, collection)
return OrderedDict(
b.name => [round(value(vars[b.name, t]), digits = 5) for t in 1:T]
for b in collection
)
end
function production_cost(g)
function timeseries_second_stage(vars, collection, sc)
return OrderedDict(
b.name => [round(value(vars[sc.name, b.name, t]), digits = 5) for t in 1:T]
for b in collection
)
end
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 +59,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_first_stage(model[:is_on], sc.units)
sol[sc.name]["Switch on"] = timeseries_first_stage(model[:switch_on], sc.units)
sol[sc.name]["Switch off"] = timeseries_first_stage(model[:switch_off], sc.units)
sol[sc.name]["Net injection (MW)"] =
timeseries_second_stage(model[:net_injection], sc.buses, sc)
sol[sc.name]["Load curtail (MW)"] = timeseries_second_stage(model[:curtail], sc.buses, sc)
if !isempty(sc.lines)
sol[sc.name]["Line overflow (MW)"] = timeseries_second_stage(model[:overflow], sc.lines, sc)
end
if !isempty(sc.price_sensitive_loads)
sol[sc.name]["Price-sensitive loads (MW)"] =
timeseries_second_stage(model[:loads], sc.price_sensitive_loads, 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 == "up-frp"
)
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 == "up-frp"
)
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 == "down-frp"
)
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 == "down-frp"
)
end
if !isempty(instance.price_sensitive_loads)
sol["Price-sensitive loads (MW)"] =
timeseries(model[:loads], instance.price_sensitive_loads)
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"
)
length(keys(sol)) > 1 ? nothing : sol = Dict(sol_key => sol_val for scen_key in keys(sol) for (sol_key, sol_val) in sol[scen_key])
return sol
end

10
src/validation/repair.jl
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]

4
src/validation/validate.jl
View File

@@ -356,7 +356,7 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
required,
)
end
elseif r.type == "flexiramp"
elseif r.type == "up-frp"
upflexiramp = sum(
solution["Up-flexiramp (MW)"][r.name][g.name][t] for
g in r.units
@@ -374,7 +374,7 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
)
err_count += 1
end
elseif r.type == "down-frp"
dwflexiramp = sum(
solution["Down-flexiramp (MW)"][r.name][g.name][t] for
g in r.units