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

Break down model.jl

Browse Source
This commit is contained in:
Alinson S. Xavier
2021-05-29 18:33:16 -05:00
parent 4e8426beba
commit 483c793d49
17 changed files with 545 additions and 465 deletions
Download Patch File Download Diff File Expand all files Collapse all files

32
benchmark/Makefile
View File

@@ -8,21 +8,21 @@ TIMESTAMP := $(shell date "+%Y-%m-%d %H:%M")
SRC_FILES := $(wildcard ../src/*.jl) SRC_FILES := $(wildcard ../src/*.jl)
INSTANCES_PGLIB := \ INSTANCES_PGLIB := \
pglib-uc/ca/2014-09-01_reserves_0 \ pglib-uc/ca/2014-09-01_reserves_0 \
pglib-uc/ca/2014-09-01_reserves_1 \ pglib-uc/ca/2014-09-01_reserves_1 \
pglib-uc/ca/2015-03-01_reserves_0 \ pglib-uc/ca/2015-03-01_reserves_0 \
pglib-uc/ca/2015-06-01_reserves_0 \ pglib-uc/ca/2015-06-01_reserves_0 \
pglib-uc/ca/Scenario400_reserves_1 \ pglib-uc/ca/Scenario400_reserves_1 \
pglib-uc/ferc/2015-01-01_lw \ pglib-uc/ferc/2015-01-01_lw \
pglib-uc/ferc/2015-05-01_lw \ pglib-uc/ferc/2015-05-01_lw \
pglib-uc/ferc/2015-07-01_hw \ pglib-uc/ferc/2015-07-01_hw \
pglib-uc/ferc/2015-10-01_lw \ pglib-uc/ferc/2015-10-01_lw \
pglib-uc/ferc/2015-12-01_lw \ pglib-uc/ferc/2015-12-01_lw \
pglib-uc/rts_gmlc/2020-04-03 \ pglib-uc/rts_gmlc/2020-04-03 \
pglib-uc/rts_gmlc/2020-09-20 \ pglib-uc/rts_gmlc/2020-09-20 \
pglib-uc/rts_gmlc/2020-10-27 \ pglib-uc/rts_gmlc/2020-10-27 \
pglib-uc/rts_gmlc/2020-11-25 \ pglib-uc/rts_gmlc/2020-11-25 \
pglib-uc/rts_gmlc/2020-12-23 pglib-uc/rts_gmlc/2020-12-23
INSTANCES_MATPOWER := \ INSTANCES_MATPOWER := \
matpower/case118/2017-02-01 \ matpower/case118/2017-02-01 \
@@ -68,7 +68,7 @@ SOLUTIONS_PGLIB := $(foreach s,$(SAMPLES),$(addprefix results/,$(addsuffix .$(s)
SOLUTIONS_ORLIB := $(foreach s,$(SAMPLES),$(addprefix results/,$(addsuffix .$(s).sol.json,$(INSTANCES_ORLIB)))) SOLUTIONS_ORLIB := $(foreach s,$(SAMPLES),$(addprefix results/,$(addsuffix .$(s).sol.json,$(INSTANCES_ORLIB))))
SOLUTIONS_TEJADA19 := $(foreach s,$(SAMPLES),$(addprefix results/,$(addsuffix .$(s).sol.json,$(INSTANCES_TEJADA19)))) SOLUTIONS_TEJADA19 := $(foreach s,$(SAMPLES),$(addprefix results/,$(addsuffix .$(s).sol.json,$(INSTANCES_TEJADA19))))
.PHONY: tables save small large clean-mps matpower pglib orlib .PHONY: matpower pglib orlib tejada19 clean clean-mps clean-sol save tables
all: matpower pglib orlib tejada19 all: matpower pglib orlib tejada19

32
benchmark/run.jl
View File

@@ -16,49 +16,37 @@ function main()
basename, suffix = split(ARGS[1], ".") basename, suffix = split(ARGS[1], ".")
solution_filename = "results/$basename.$suffix.sol.json" solution_filename = "results/$basename.$suffix.sol.json"
model_filename = "results/$basename.$suffix.mps.gz" model_filename = "results/$basename.$suffix.mps.gz"
time_limit = 60 * 20
BLAS.set_num_threads(4) BLAS.set_num_threads(4)
total_time = @elapsed begin total_time = @elapsed begin
@info "Reading: $basename" @info "Reading: $basename"
time_read = @elapsed begin time_read = @elapsed begin
instance = UnitCommitment.read_benchmark(basename) instance = UnitCommitment.read_benchmark(basename)
end end
@info @sprintf("Read problem in %.2f seconds", time_read) @info @sprintf("Read problem in %.2f seconds", time_read)
model = UnitCommitment.build_model(
time_model = @elapsed begin instance = instance,
model = build_model( optimizer = optimizer_with_attributes(
instance = instance, Gurobi.Optimizer,
optimizer = optimizer_with_attributes( "Threads" => 4,
Gurobi.Optimizer, "Seed" => rand(1:1000),
"Threads" => 4, ),
"Seed" => rand(1:1000), variable_names = true,
), )
variable_names = true,
)
end
@info "Optimizing..." @info "Optimizing..."
BLAS.set_num_threads(1) BLAS.set_num_threads(1)
UnitCommitment.optimize!( UnitCommitment.optimize!(
model, model,
time_limit = time_limit, UnitCommitment._XaQiWaTh19(time_limit = 3600.0),
gap_limit = 1e-3,
) )
end end
@info @sprintf("Total time was %.2f seconds", total_time) @info @sprintf("Total time was %.2f seconds", total_time)
@info "Writing: $solution_filename" @info "Writing: $solution_filename"
solution = UnitCommitment.solution(model) solution = UnitCommitment.solution(model)
open(solution_filename, "w") do file open(solution_filename, "w") do file
return JSON.print(file, solution, 2) return JSON.print(file, solution, 2)
end end
@info "Verifying solution..." @info "Verifying solution..."
UnitCommitment.validate(instance, solution) UnitCommitment.validate(instance, solution)
@info "Exporting model..." @info "Exporting model..."
return JuMP.write_to_file(model, model_filename) return JuMP.write_to_file(model, model_filename)
end end

10
src/UnitCommitment.jl
View File

@@ -5,13 +5,17 @@
module UnitCommitment module UnitCommitment
include("instance/structs.jl") include("instance/structs.jl")
include("transmission/structs.jl")
include("solution/structs.jl") include("solution/structs.jl")
include("solution/methods/XaQiWaTh19/structs.jl") include("solution/methods/XaQiWaTh19/structs.jl")
include("import/egret.jl") include("import/egret.jl")
include("instance/read.jl") include("instance/read.jl")
include("model/build.jl") include("model/build.jl")
include("model/formulations/base/bus.jl")
include("model/formulations/base/line.jl")
include("model/formulations/base/psload.jl")
include("model/formulations/base/system.jl")
include("model/formulations/base/unit.jl")
include("model/jumpext.jl") include("model/jumpext.jl")
include("solution/fix.jl") include("solution/fix.jl")
include("solution/methods/XaQiWaTh19/enforce.jl") include("solution/methods/XaQiWaTh19/enforce.jl")
@@ -22,8 +26,8 @@ include("solution/optimize.jl")
include("solution/solution.jl") include("solution/solution.jl")
include("solution/warmstart.jl") include("solution/warmstart.jl")
include("solution/write.jl") include("solution/write.jl")
include("transforms/initcond.jl") include("transform/initcond.jl")
include("transforms/slice.jl") include("transform/slice.jl")
include("transmission/sensitivity.jl") include("transmission/sensitivity.jl")
include("utils/log.jl") include("utils/log.jl")
include("validation/repair.jl") include("validation/repair.jl")

409
src/model/build.jl
View File

@@ -75,7 +75,6 @@ function build_model(;
) )
end end
@info @sprintf("Computed ISF in %.2f seconds", time_isf) @info @sprintf("Computed ISF in %.2f seconds", time_isf)
@info "Computing line outage factors..." @info "Computing line outage factors..."
time_lodf = @elapsed begin time_lodf = @elapsed begin
lodf = UnitCommitment._line_outage_factors( lodf = UnitCommitment._line_outage_factors(
@@ -95,7 +94,6 @@ function build_model(;
lodf[abs.(lodf).<lodf_cutoff] .= 0 lodf[abs.(lodf).<lodf_cutoff] .= 0
end end
end end
@info "Building model..." @info "Building model..."
time_model = @elapsed begin time_model = @elapsed begin
model = Model() model = Model()
@@ -106,411 +104,16 @@ function build_model(;
model[:instance] = instance model[:instance] = instance
model[:isf] = isf model[:isf] = isf
model[:lodf] = lodf model[:lodf] = lodf
for field in [ _add_transmission_line!.(model, instance.lines)
:prod_above, _add_bus!.(model, instance.buses)
:segprod, _add_unit!.(model, instance.units)
:reserve, _add_price_sensitive_load!.(model, instance.price_sensitive_loads)
:is_on, _add_system_wide_eqs!(model)
:switch_on, @objective(model, Min, model[:obj])
:switch_off,
:net_injection,
:curtail,
:overflow,
:loads,
:startup,
:eq_startup_choose,
:eq_startup_restrict,
:eq_segprod_limit,
:eq_prod_above_def,
:eq_prod_limit,
:eq_binary_link,
:eq_switch_on_off,
:eq_ramp_up,
:eq_ramp_down,
:eq_startup_limit,
:eq_shutdown_limit,
:eq_min_uptime,
:eq_min_downtime,
:eq_power_balance,
:eq_net_injection_def,
:eq_min_reserve,
:expr_inj,
:expr_reserve,
:expr_net_injection,
]
model[field] = OrderedDict()
end
for lm in instance.lines
_add_transmission_line!(model, lm)
end
for b in instance.buses
_add_bus!(model, b)
end
for g in instance.units
_add_unit!(model, g)
end
for ps in instance.price_sensitive_loads
_add_price_sensitive_load!(model, ps)
end
_build_net_injection_eqs!(model)
_build_reserve_eqs!(model)
_build_obj_function!(model)
end end
@info @sprintf("Built model in %.2f seconds", time_model) @info @sprintf("Built model in %.2f seconds", time_model)
if variable_names if variable_names
_set_names!(model) _set_names!(model)
end end
return model return model
end end
function _add_transmission_line!(model, lm)
obj, T = model[:obj], model[:instance].time
overflow = model[:overflow]
for t in 1:T
v = overflow[lm.name, t] = @variable(model, lower_bound = 0)
add_to_expression!(obj, v, lm.flow_limit_penalty[t])
end
end
function _add_bus!(model::JuMP.Model, b::Bus)
mip = model
net_injection = model[:expr_net_injection]
reserve = model[:expr_reserve]
curtail = model[:curtail]
for t in 1:model[:instance].time
# Fixed load
net_injection[b.name, t] = AffExpr(-b.load[t])
# Reserves
reserve[b.name, t] = AffExpr()
# Load curtailment
curtail[b.name, t] =
@variable(mip, 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!(
model[:obj],
curtail[b.name, t],
model[:instance].power_balance_penalty[t],
)
end
end
function _add_price_sensitive_load!(model::JuMP.Model, ps::PriceSensitiveLoad)
mip = model
loads = model[:loads]
net_injection = model[:expr_net_injection]
for t in 1:model[:instance].time
# Decision variable
loads[ps.name, t] =
@variable(mip, lower_bound = 0, upper_bound = ps.demand[t])
# Objective function terms
add_to_expression!(model[:obj], loads[ps.name, t], -ps.revenue[t])
# Net injection
add_to_expression!(
net_injection[ps.bus.name, t],
loads[ps.name, t],
-1.0,
)
end
end
function _add_unit!(model::JuMP.Model, g::Unit)
mip, T = model, model[:instance].time
gi, K, S = g.name, length(g.cost_segments), length(g.startup_categories)
segprod = model[:segprod]
prod_above = model[:prod_above]
reserve = model[:reserve]
startup = model[:startup]
is_on = model[:is_on]
switch_on = model[:switch_on]
switch_off = model[:switch_off]
expr_net_injection = model[:expr_net_injection]
expr_reserve = model[:expr_reserve]
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
is_initially_on = (g.initial_status > 0 ? 1.0 : 0.0)
# Decision variables
for t in 1:T
for k in 1:K
segprod[gi, t, k] = @variable(model, lower_bound = 0)
end
prod_above[gi, t] = @variable(model, lower_bound = 0)
if g.provides_spinning_reserves[t]
reserve[gi, t] = @variable(model, lower_bound = 0)
else
reserve[gi, t] = 0.0
end
for s in 1:S
startup[gi, t, s] = @variable(model, binary = true)
end
if g.must_run[t]
is_on[gi, t] = 1.0
switch_on[gi, t] = (t == 1 ? 1.0 - is_initially_on : 0.0)
switch_off[gi, t] = 0.0
else
is_on[gi, t] = @variable(model, binary = true)
switch_on[gi, t] = @variable(model, binary = true)
switch_off[gi, t] = @variable(model, binary = true)
end
end
for t in 1:T
# Time-dependent start-up costs
for s in 1:S
# If unit is switching on, we must choose a startup category
model[:eq_startup_choose][gi, t, s] = @constraint(
mip,
switch_on[gi, t] == sum(startup[gi, t, s] for s in 1:S)
)
# If unit has not switched off in the last `delay` time periods, startup category is forbidden.
# The last startup category is always allowed.
if s < S
range_start = t - g.startup_categories[s+1].delay + 1
range_end = t - g.startup_categories[s].delay
range = (range_start:range_end)
initial_sum = (
g.initial_status < 0 && (g.initial_status + 1 in range) ? 1.0 : 0.0
)
model[:eq_startup_restrict][gi, t, s] = @constraint(
mip,
startup[gi, t, s] <=
initial_sum +
sum(switch_off[gi, i] for i in range if i >= 1)
)
end
# Objective function terms for start-up costs
add_to_expression!(
model[:obj],
startup[gi, t, s],
g.startup_categories[s].cost,
)
end
# Objective function terms for production costs
add_to_expression!(model[:obj], is_on[gi, t], g.min_power_cost[t])
for k in 1:K
add_to_expression!(
model[:obj],
segprod[gi, t, k],
g.cost_segments[k].cost[t],
)
end
# Production limits (piecewise-linear segments)
for k in 1:K
model[:eq_segprod_limit][gi, t, k] = @constraint(
mip,
segprod[gi, t, k] <= g.cost_segments[k].mw[t] * is_on[gi, t]
)
end
# Definition of production
model[:eq_prod_above_def][gi, t] = @constraint(
mip,
prod_above[gi, t] == sum(segprod[gi, t, k] for k in 1:K)
)
# Production limit
model[:eq_prod_limit][gi, t] = @constraint(
mip,
prod_above[gi, t] + reserve[gi, t] <=
(g.max_power[t] - g.min_power[t]) * is_on[gi, t]
)
# Binary variable equations for economic units
if !g.must_run[t]
# Link binary variables
if t == 1
model[:eq_binary_link][gi, t] = @constraint(
mip,
is_on[gi, t] - is_initially_on ==
switch_on[gi, t] - switch_off[gi, t]
)
else
model[:eq_binary_link][gi, t] = @constraint(
mip,
is_on[gi, t] - is_on[gi, t-1] ==
switch_on[gi, t] - switch_off[gi, t]
)
end
# Cannot switch on and off at the same time
model[:eq_switch_on_off][gi, t] =
@constraint(mip, switch_on[gi, t] + switch_off[gi, t] <= 1)
end
# Ramp up limit
if t == 1
if is_initially_on == 1
model[:eq_ramp_up][gi, t] = @constraint(
mip,
prod_above[gi, t] + reserve[gi, t] <=
(g.initial_power - g.min_power[t]) + g.ramp_up_limit
)
end
else
model[:eq_ramp_up][gi, t] = @constraint(
mip,
prod_above[gi, t] + reserve[gi, t] <=
prod_above[gi, t-1] + g.ramp_up_limit
)
end
# Ramp down limit
if t == 1
if is_initially_on == 1
model[:eq_ramp_down][gi, t] = @constraint(
mip,
prod_above[gi, t] >=
(g.initial_power - g.min_power[t]) - g.ramp_down_limit
)
end
else
model[:eq_ramp_down][gi, t] = @constraint(
mip,
prod_above[gi, t] >= prod_above[gi, t-1] - g.ramp_down_limit
)
end
# Startup limit
model[:eq_startup_limit][gi, t] = @constraint(
mip,
prod_above[gi, t] + reserve[gi, t] <=
(g.max_power[t] - g.min_power[t]) * is_on[gi, t] -
max(0, g.max_power[t] - g.startup_limit) * switch_on[gi, t]
)
# Shutdown limit
if g.initial_power > g.shutdown_limit
model[:eq_shutdown_limit][gi, 0] =
@constraint(mip, switch_off[gi, 1] <= 0)
end
if t < T
model[:eq_shutdown_limit][gi, t] = @constraint(
mip,
prod_above[gi, t] <=
(g.max_power[t] - g.min_power[t]) * is_on[gi, t] -
max(0, g.max_power[t] - g.shutdown_limit) * switch_off[gi, t+1]
)
end
# Minimum up-time
model[:eq_min_uptime][gi, t] = @constraint(
mip,
sum(switch_on[gi, i] for i in (t-g.min_uptime+1):t if i >= 1) <=
is_on[gi, t]
)
# # Minimum down-time
model[:eq_min_downtime][gi, t] = @constraint(
mip,
sum(switch_off[gi, i] for i in (t-g.min_downtime+1):t if i >= 1) <= 1 - is_on[gi, t]
)
# Minimum up/down-time for initial periods
if t == 1
if g.initial_status > 0
model[:eq_min_uptime][gi, 0] = @constraint(
mip,
sum(
switch_off[gi, i] for
i in 1:(g.min_uptime-g.initial_status) if i <= T
) == 0
)
else
model[:eq_min_downtime][gi, 0] = @constraint(
mip,
sum(
switch_on[gi, i] for
i in 1:(g.min_downtime+g.initial_status) if i <= T
) == 0
)
end
end
# Add to net injection expression
add_to_expression!(
expr_net_injection[g.bus.name, t],
prod_above[g.name, t],
1.0,
)
add_to_expression!(
expr_net_injection[g.bus.name, t],
is_on[g.name, t],
g.min_power[t],
)
# Add to reserves expression
add_to_expression!(expr_reserve[g.bus.name, t], reserve[gi, t], 1.0)
end
end
function _build_obj_function!(model::JuMP.Model)
@objective(model, Min, model[:obj])
end
function _build_net_injection_eqs!(model::JuMP.Model)
T = model[:instance].time
net_injection = model[:net_injection]
for t in 1:T, b in model[:instance].buses
n = net_injection[b.name, t] = @variable(model)
model[:eq_net_injection_def][t, b.name] =
@constraint(model, n == model[:expr_net_injection][b.name, t])
end
for t in 1:T
model[:eq_power_balance][t] = @constraint(
model,
sum(net_injection[b.name, t] for b in model[:instance].buses) == 0
)
end
end
function _build_reserve_eqs!(model::JuMP.Model)
reserves = model[:instance].reserves
for t in 1:model[:instance].time
model[:eq_min_reserve][t] = @constraint(
model,
sum(
model[:expr_reserve][b.name, t] for b in model[:instance].buses
) >= reserves.spinning[t]
)
end
end
function _set_names!(model::JuMP.Model)
@info "Setting variable and constraint names..."
time_varnames = @elapsed begin
_set_names!(object_dictionary(model))
end
@info @sprintf("Set names in %.2f seconds", time_varnames)
end
function _set_names!(dict::Dict)
for name in keys(dict)
dict[name] isa AbstractDict || continue
for idx in keys(dict[name])
if dict[name][idx] isa AffExpr
continue
end
idx_str = join(map(string, idx), ",")
set_name(dict[name][idx], "$name[$idx_str]")
end
end
end

28
src/model/formulations/base/bus.jl Normal file
View File

@@ -0,0 +1,28 @@
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# 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
net_injection = _get(model, :expr_net_injection)
reserve = _get(model, :expr_reserve)
curtail = _get(model, :curtail)
for t in 1:model[:instance].time
# Fixed load
net_injection[b.name, t] = AffExpr(-b.load[t])
# Reserves
reserve[b.name, t] = AffExpr()
# Load curtailment
curtail[b.name, t] =
@variable(model, lower_bound = 0, upper_bound = b.load[t])
add_to_expression!(net_injection[b.name, t], curtail[b.name, t], 1.0)
add_to_expression!(
model[:obj],
curtail[b.name, t],
model[:instance].power_balance_penalty[t],
)
end
return
end

12
src/model/formulations/base/line.jl Normal file
View File

@@ -0,0 +1,12 @@
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
function _add_transmission_line!(model, lm)::Nothing
overflow = _get(model, :overflow)
for t in 1:model[:instance].time
v = overflow[lm.name, t] = @variable(model, lower_bound = 0)
add_to_expression!(model[:obj], v, lm.flow_limit_penalty[t])
end
return
end

27
src/model/formulations/base/psload.jl Normal file
View File

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

41
src/model/formulations/base/system.jl Normal file
View File

@@ -0,0 +1,41 @@
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# 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_reserve_eqs!(model)
return
end
function _add_net_injection_eqs!(model::JuMP.Model)::Nothing
T = model[:instance].time
net_injection = _get(model, :net_injection)
eq_net_injection_def = _get(model, :eq_net_injection_def)
eq_power_balance = _get(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_def[t, b.name] =
@constraint(model, n == model[:expr_net_injection][b.name, t])
end
for t in 1:T
eq_power_balance[t] = @constraint(
model,
sum(net_injection[b.name, t] for b in model[:instance].buses) == 0
)
end
return
end
function _add_reserve_eqs!(model::JuMP.Model)::Nothing
eq_min_reserve = _get(model, :eq_min_reserve)
for t in 1:model[:instance].time
eq_min_reserve[t] = @constraint(
model,
sum(
model[:expr_reserve][b.name, t] for b in model[:instance].buses
) >= model[:instance].reserves.spinning[t]
)
end
return
end

349
src/model/formulations/base/unit.jl Normal file
View File

@@ -0,0 +1,349 @@
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# 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)
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)
_add_reserve_vars!(model, g)
_add_startup_shutdown_vars!(model, g)
_add_status_vars!(model, g)
# Constraints and objective function
_add_min_uptime_downtime_eqs!(model, g)
_add_net_injection_eqs!(model, g)
_add_production_eqs!(model, g)
_add_ramp_eqs!(model, g)
_add_startup_shutdown_costs_eqs!(model, g)
_add_startup_shutdown_limit_eqs!(model, g)
return _add_status_eqs!(model, g)
end
_is_initially_on(g::Unit)::Float64 = (g.initial_status > 0 ? 1.0 : 0.0)
function _add_production_vars!(model::JuMP.Model, g::Unit)::Nothing
prod_above = _get(model, :prod_above)
segprod = _get(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)
end
prod_above[g.name, t] = @variable(model, lower_bound = 0)
end
return
end
function _add_production_eqs!(model::JuMP.Model, g::Unit)::Nothing
eq_prod_above_def = _get(model, :eq_prod_above_def)
eq_prod_limit = _get(model, :eq_prod_limit)
eq_segprod_limit = _get(model, :eq_segprod_limit)
is_on = model[:is_on]
K = length(g.cost_segments)
prod_above = model[:prod_above]
reserve = model[:reserve]
segprod = model[:segprod]
for t in 1:model[:instance].time
# Objective function terms for production costs
add_to_expression!(model[:obj], is_on[g.name, t], g.min_power_cost[t])
for k in 1:K
add_to_expression!(
model[:obj],
segprod[g.name, t, k],
g.cost_segments[k].cost[t],
)
end
# Production limits (piecewise-linear segments)
for k in 1:K
eq_segprod_limit[g.name, t, k] = @constraint(
model,
segprod[g.name, t, k] <=
g.cost_segments[k].mw[t] * is_on[g.name, t]
)
end
# Definition of production
eq_prod_above_def[g.name, t] = @constraint(
model,
prod_above[g.name, t] == sum(segprod[g.name, t, k] for k in 1:K)
)
# Production limit
eq_prod_limit[g.name, t] = @constraint(
model,
prod_above[g.name, t] + reserve[g.name, t] <=
(g.max_power[t] - g.min_power[t]) * is_on[g.name, t]
)
end
return
end
function _add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
reserve = _get(model, :reserve)
for t in 1:model[:instance].time
if g.provides_spinning_reserves[t]
reserve[g.name, t] = @variable(model, lower_bound = 0)
else
reserve[g.name, t] = 0.0
end
end
return
end
function _add_reserve_eqs!(model::JuMP.Model, g::Unit)::Nothing
reserve = model[:reserve]
for t in 1:model[:instance].time
add_to_expression!(expr_reserve[g.bus.name, t], reserve[g.name, t], 1.0)
end
return
end
function _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing
startup = _get(model, :startup)
for t in 1:model[:instance].time
for s in 1:length(g.startup_categories)
startup[g.name, t, s] = @variable(model, binary = true)
end
end
return
end
function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
eq_shutdown_limit = _get(model, :eq_shutdown_limit)
eq_startup_limit = _get(model, :eq_startup_limit)
is_on = model[:is_on]
prod_above = model[:prod_above]
reserve = model[:reserve]
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(
model,
prod_above[g.name, t] + reserve[g.name, 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] =
@constraint(model, switch_off[g.name, 1] <= 0)
end
if t < T
eq_shutdown_limit[g.name, t] = @constraint(
model,
prod_above[g.name, t] <=
(g.max_power[t] - g.min_power[t]) * is_on[g.name, t] -
max(0, g.max_power[t] - g.shutdown_limit) *
switch_off[g.name, t+1]
)
end
end
return
end
function _add_startup_shutdown_costs_eqs!(model::JuMP.Model, g::Unit)::Nothing
eq_startup_choose = _get(model, :eq_startup_choose)
eq_startup_restrict = _get(model, :eq_startup_restrict)
S = length(g.startup_categories)
startup = model[:startup]
for t in 1:model[:instance].time
for s in 1:S
# If unit is switching on, we must choose a startup category
eq_startup_choose[g.name, t, s] = @constraint(
model,
model[:switch_on][g.name, t] ==
sum(startup[g.name, t, s] for s in 1:S)
)
# If unit has not switched off in the last `delay` time periods, startup category is forbidden.
# The last startup category is always allowed.
if s < S
range_start = t - g.startup_categories[s+1].delay + 1
range_end = t - g.startup_categories[s].delay
range = (range_start:range_end)
initial_sum = (
g.initial_status < 0 && (g.initial_status + 1 in range) ? 1.0 : 0.0
)
eq_startup_restrict[g.name, t, s] = @constraint(
model,
startup[g.name, t, s] <=
initial_sum + sum(
model[:switch_off][g.name, i] for i in range if i >= 1
)
)
end
# Objective function terms for start-up costs
add_to_expression!(
model[:obj],
startup[g.name, t, s],
g.startup_categories[s].cost,
)
end
end
return
end
function _add_status_vars!(model::JuMP.Model, g::Unit)::Nothing
is_on = _get(model, :is_on)
switch_on = _get(model, :switch_on)
switch_off = _get(model, :switch_off)
for t in 1:model[:instance].time
if g.must_run[t]
is_on[g.name, t] = 1.0
switch_on[g.name, t] = (t == 1 ? 1.0 - _is_initially_on(g) : 0.0)
switch_off[g.name, t] = 0.0
else
is_on[g.name, t] = @variable(model, binary = true)
switch_on[g.name, t] = @variable(model, binary = true)
switch_off[g.name, t] = @variable(model, binary = true)
end
end
return
end
function _add_status_eqs!(model::JuMP.Model, g::Unit)::Nothing
eq_binary_link = _get(model, :eq_binary_link)
eq_switch_on_off = _get(model, :eq_switch_on_off)
is_on = model[:is_on]
switch_off = model[:switch_off]
switch_on = model[:switch_on]
for t in 1:model[:instance].time
if !g.must_run[t]
# Link binary variables
if t == 1
eq_binary_link[g.name, t] = @constraint(
model,
is_on[g.name, t] - _is_initially_on(g) ==
switch_on[g.name, t] - switch_off[g.name, t]
)
else
eq_binary_link[g.name, t] = @constraint(
model,
is_on[g.name, t] - is_on[g.name, t-1] ==
switch_on[g.name, t] - switch_off[g.name, t]
)
end
# Cannot switch on and off at the same time
eq_switch_on_off[g.name, t] = @constraint(
model,
switch_on[g.name, t] + switch_off[g.name, t] <= 1
)
end
end
return
end
function _add_ramp_eqs!(model::JuMP.Model, g::Unit)::Nothing
prod_above = model[:prod_above]
reserve = model[:reserve]
eq_ramp_up = _get(model, :eq_ramp_up)
eq_ramp_down = _get(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(
model,
prod_above[g.name, t] + reserve[g.name, t] <=
(g.initial_power - g.min_power[t]) + g.ramp_up_limit
)
end
else
eq_ramp_up[g.name, t] = @constraint(
model,
prod_above[g.name, t] + reserve[g.name, t] <=
prod_above[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(
model,
prod_above[g.name, t] >=
(g.initial_power - g.min_power[t]) - g.ramp_down_limit
)
end
else
eq_ramp_down[g.name, t] = @constraint(
model,
prod_above[g.name, t] >=
prod_above[g.name, t-1] - g.ramp_down_limit
)
end
end
end
function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
is_on = model[:is_on]
switch_off = model[:switch_off]
switch_on = model[:switch_on]
eq_min_uptime = _get(model, :eq_min_uptime)
eq_min_downtime = _get(model, :eq_min_downtime)
T = model[:instance].time
for t in 1:T
# Minimum up-time
eq_min_uptime[g.name, t] = @constraint(
model,
sum(switch_on[g.name, i] for i in (t-g.min_uptime+1):t if i >= 1) <= is_on[g.name, t]
)
# Minimum down-time
eq_min_downtime[g.name, t] = @constraint(
model,
sum(
switch_off[g.name, i] for i in (t-g.min_downtime+1):t if i >= 1
) <= 1 - is_on[g.name, t]
)
# Minimum up/down-time for initial periods
if t == 1
if g.initial_status > 0
eq_min_uptime[g.name, 0] = @constraint(
model,
sum(
switch_off[g.name, i] for
i in 1:(g.min_uptime-g.initial_status) if i <= T
) == 0
)
else
eq_min_downtime[g.name, 0] = @constraint(
model,
sum(
switch_on[g.name, i] for
i in 1:(g.min_downtime+g.initial_status) if i <= T
) == 0
)
end
end
end
end
function _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing
expr_net_injection = model[:expr_net_injection]
expr_reserve = model[:expr_reserve]
is_on = model[:is_on]
prod_above = model[:prod_above]
reserve = model[:reserve]
for t in 1:model[:instance].time
# Add to net injection expression
add_to_expression!(
expr_net_injection[g.bus.name, t],
prod_above[g.name, t],
1.0,
)
add_to_expression!(
expr_net_injection[g.bus.name, t],
is_on[g.name, t],
g.min_power[t],
)
# Add to reserves expression
add_to_expression!(expr_reserve[g.bus.name, t], reserve[g.name, t], 1.0)
end
end

28
src/model/jumpext.jl
View File

@@ -18,3 +18,31 @@ end
function set_name(x::Float64, n::String) function set_name(x::Float64, n::String)
# nop # nop
end end
function _get(model::JuMP.Model, key::Symbol)::OrderedDict
if !(key in keys(object_dictionary(model)))
model[key] = OrderedDict()
end
return model[key]
end
function _set_names!(model::JuMP.Model)
@info "Setting variable and constraint names..."
time_varnames = @elapsed begin
_set_names!(object_dictionary(model))
end
@info @sprintf("Set names in %.2f seconds", time_varnames)
end
function _set_names!(dict::Dict)
for name in keys(dict)
dict[name] isa AbstractDict || continue
for idx in keys(dict[name])
if dict[name][idx] isa AffExpr
continue
end
idx_str = join(map(string, idx), ",")
set_name(dict[name][idx], "$name[$idx_str]")
end
end
end

2
src/transmission/structs.jl → src/model/structs.jl
View File

@@ -1,3 +1,5 @@
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment # UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details. # Released under the modified BSD license. See COPYING.md for more details.
abstract type Formulation end

4
src/solution/fix.jl
View File

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

23
src/solution/methods/XaQiWaTh19/structs.jl
View File

@@ -9,6 +9,8 @@ import DataStructures: PriorityQueue
time_limit::Float64 time_limit::Float64
gap_limit::Float64 gap_limit::Float64
two_phase_gap::Bool two_phase_gap::Bool
max_violations_per_line::Int
max_violations_per_period::Int
end end
Lazy constraint solution method described in: Lazy constraint solution method described in:
@@ -17,8 +19,8 @@ Lazy constraint solution method described in:
constraint filtering in large-scale security-constrained unit commitment. constraint filtering in large-scale security-constrained unit commitment.
IEEE Transactions on Power Systems, 34(3), 2457-2460. IEEE Transactions on Power Systems, 34(3), 2457-2460.
Fields ## Fields
=========
- `time_limit`: - `time_limit`:
the time limit over the entire optimization procedure. the time limit over the entire optimization procedure.
- `gap_limit`: - `gap_limit`:
@@ -26,6 +28,13 @@ Fields
- `two_phase_gap`: - `two_phase_gap`:
if true, solve the problem with large gap tolerance first, then reduce if true, solve the problem with large gap tolerance first, then reduce
the gap tolerance when no further violated constraints are found. the gap tolerance when no further violated constraints are found.
- `max_violations_per_line`:
maximum number of violated transmission constraints to add to the
formulation per transmission line.
- `max_violations_per_period`:
maximum number of violated transmission constraints to add to the
formulation per time period.
""" """
struct _XaQiWaTh19 struct _XaQiWaTh19
time_limit::Float64 time_limit::Float64
@@ -35,11 +44,11 @@ struct _XaQiWaTh19
max_violations_per_period::Int max_violations_per_period::Int
function _XaQiWaTh19(; function _XaQiWaTh19(;
time_limit::Float64, time_limit::Float64 = 86400.0,
gap_limit::Float64, gap_limit::Float64 = 1e-3,
two_phase_gap::Bool, two_phase_gap::Bool = true,
max_violations_per_line::Int, max_violations_per_line::Int = 1,
max_violations_per_period::Int, max_violations_per_period::Int = 5,
) )
return new( return new(
time_limit, time_limit,

11
src/solution/optimize.jl
View File

@@ -10,14 +10,5 @@ advanced methods to accelerate the solution process and to enforce transmission
and N-1 security constraints. and N-1 security constraints.
""" """
function optimize!(model::JuMP.Model)::Nothing function optimize!(model::JuMP.Model)::Nothing
return UnitCommitment.optimize!( return UnitCommitment.optimize!(model, _XaQiWaTh19())
model,
_XaQiWaTh19(
time_limit = 3600.0,
gap_limit = 1e-4,
two_phase_gap = true,
max_violations_per_line = 1,
max_violations_per_period = 5,
),
)
end end

2
src/solution/warmstart.jl
View File

@@ -5,8 +5,6 @@
function set_warm_start!(model::JuMP.Model, solution::AbstractDict)::Nothing function set_warm_start!(model::JuMP.Model, solution::AbstractDict)::Nothing
instance, T = model[:instance], model[:instance].time instance, T = model[:instance], model[:instance].time
is_on = model[:is_on] is_on = model[:is_on]
prod_above = model[:prod_above]
reserve = model[:reserve]
for g in instance.units for g in instance.units
for t in 1:T for t in 1:T
JuMP.set_start_value(is_on[g.name, t], solution["Is on"][g.name][t]) JuMP.set_start_value(is_on[g.name, t], solution["Is on"][g.name][t])

0
src/transforms/initcond.jl → src/transform/initcond.jl
View File

0
src/transforms/slice.jl → src/transform/slice.jl
View File