diff --git a/benchmark/Makefile b/benchmark/Makefile index 12418b8..e71a55c 100644 --- a/benchmark/Makefile +++ b/benchmark/Makefile @@ -8,21 +8,21 @@ TIMESTAMP := $(shell date "+%Y-%m-%d %H:%M") SRC_FILES := $(wildcard ../src/*.jl) INSTANCES_PGLIB := \ - pglib-uc/ca/2014-09-01_reserves_0 \ - pglib-uc/ca/2014-09-01_reserves_1 \ - pglib-uc/ca/2015-03-01_reserves_0 \ - pglib-uc/ca/2015-06-01_reserves_0 \ - pglib-uc/ca/Scenario400_reserves_1 \ - pglib-uc/ferc/2015-01-01_lw \ - pglib-uc/ferc/2015-05-01_lw \ - pglib-uc/ferc/2015-07-01_hw \ - pglib-uc/ferc/2015-10-01_lw \ - pglib-uc/ferc/2015-12-01_lw \ - pglib-uc/rts_gmlc/2020-04-03 \ - pglib-uc/rts_gmlc/2020-09-20 \ - pglib-uc/rts_gmlc/2020-10-27 \ - pglib-uc/rts_gmlc/2020-11-25 \ - pglib-uc/rts_gmlc/2020-12-23 + pglib-uc/ca/2014-09-01_reserves_0 \ + pglib-uc/ca/2014-09-01_reserves_1 \ + pglib-uc/ca/2015-03-01_reserves_0 \ + pglib-uc/ca/2015-06-01_reserves_0 \ + pglib-uc/ca/Scenario400_reserves_1 \ + pglib-uc/ferc/2015-01-01_lw \ + pglib-uc/ferc/2015-05-01_lw \ + pglib-uc/ferc/2015-07-01_hw \ + pglib-uc/ferc/2015-10-01_lw \ + pglib-uc/ferc/2015-12-01_lw \ + pglib-uc/rts_gmlc/2020-04-03 \ + pglib-uc/rts_gmlc/2020-09-20 \ + pglib-uc/rts_gmlc/2020-10-27 \ + pglib-uc/rts_gmlc/2020-11-25 \ + pglib-uc/rts_gmlc/2020-12-23 INSTANCES_MATPOWER := \ 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_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 diff --git a/benchmark/run.jl b/benchmark/run.jl index 60c35c0..ac950a6 100644 --- a/benchmark/run.jl +++ b/benchmark/run.jl @@ -16,49 +16,37 @@ function main() basename, suffix = split(ARGS[1], ".") solution_filename = "results/$basename.$suffix.sol.json" model_filename = "results/$basename.$suffix.mps.gz" - - time_limit = 60 * 20 - BLAS.set_num_threads(4) - total_time = @elapsed begin @info "Reading: $basename" time_read = @elapsed begin instance = UnitCommitment.read_benchmark(basename) end @info @sprintf("Read problem in %.2f seconds", time_read) - - time_model = @elapsed begin - model = build_model( - instance = instance, - optimizer = optimizer_with_attributes( - Gurobi.Optimizer, - "Threads" => 4, - "Seed" => rand(1:1000), - ), - variable_names = true, - ) - end - + model = UnitCommitment.build_model( + instance = instance, + optimizer = optimizer_with_attributes( + Gurobi.Optimizer, + "Threads" => 4, + "Seed" => rand(1:1000), + ), + variable_names = true, + ) @info "Optimizing..." BLAS.set_num_threads(1) UnitCommitment.optimize!( model, - time_limit = time_limit, - gap_limit = 1e-3, + UnitCommitment._XaQiWaTh19(time_limit = 3600.0), ) end @info @sprintf("Total time was %.2f seconds", total_time) - @info "Writing: $solution_filename" solution = UnitCommitment.solution(model) open(solution_filename, "w") do file return JSON.print(file, solution, 2) end - @info "Verifying solution..." UnitCommitment.validate(instance, solution) - @info "Exporting model..." return JuMP.write_to_file(model, model_filename) end diff --git a/src/UnitCommitment.jl b/src/UnitCommitment.jl index 51858c5..e7692ae 100644 --- a/src/UnitCommitment.jl +++ b/src/UnitCommitment.jl @@ -5,13 +5,17 @@ module UnitCommitment include("instance/structs.jl") -include("transmission/structs.jl") include("solution/structs.jl") include("solution/methods/XaQiWaTh19/structs.jl") include("import/egret.jl") include("instance/read.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("solution/fix.jl") include("solution/methods/XaQiWaTh19/enforce.jl") @@ -22,8 +26,8 @@ include("solution/optimize.jl") include("solution/solution.jl") include("solution/warmstart.jl") include("solution/write.jl") -include("transforms/initcond.jl") -include("transforms/slice.jl") +include("transform/initcond.jl") +include("transform/slice.jl") include("transmission/sensitivity.jl") include("utils/log.jl") include("validation/repair.jl") diff --git a/src/model/build.jl b/src/model/build.jl index 68be9cb..9ff8fe4 100644 --- a/src/model/build.jl +++ b/src/model/build.jl @@ -75,7 +75,6 @@ function build_model(; ) end @info @sprintf("Computed ISF in %.2f seconds", time_isf) - @info "Computing line outage factors..." time_lodf = @elapsed begin lodf = UnitCommitment._line_outage_factors( @@ -95,7 +94,6 @@ function build_model(; lodf[abs.(lodf). 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 diff --git a/src/model/formulations/base/bus.jl b/src/model/formulations/base/bus.jl new file mode 100644 index 0000000..cc176ca --- /dev/null +++ b/src/model/formulations/base/bus.jl @@ -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 diff --git a/src/model/formulations/base/line.jl b/src/model/formulations/base/line.jl new file mode 100644 index 0000000..8370336 --- /dev/null +++ b/src/model/formulations/base/line.jl @@ -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 diff --git a/src/model/formulations/base/psload.jl b/src/model/formulations/base/psload.jl new file mode 100644 index 0000000..6303f4b --- /dev/null +++ b/src/model/formulations/base/psload.jl @@ -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 diff --git a/src/model/formulations/base/system.jl b/src/model/formulations/base/system.jl new file mode 100644 index 0000000..8929680 --- /dev/null +++ b/src/model/formulations/base/system.jl @@ -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 diff --git a/src/model/formulations/base/unit.jl b/src/model/formulations/base/unit.jl new file mode 100644 index 0000000..dd33257 --- /dev/null +++ b/src/model/formulations/base/unit.jl @@ -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 diff --git a/src/model/jumpext.jl b/src/model/jumpext.jl index b493e84..1586e92 100644 --- a/src/model/jumpext.jl +++ b/src/model/jumpext.jl @@ -18,3 +18,31 @@ end function set_name(x::Float64, n::String) # nop 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 diff --git a/src/transmission/structs.jl b/src/model/structs.jl similarity index 87% rename from src/transmission/structs.jl rename to src/model/structs.jl index 260195c..567bec2 100644 --- a/src/transmission/structs.jl +++ b/src/model/structs.jl @@ -1,3 +1,5 @@ # 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. + +abstract type Formulation end diff --git a/src/solution/fix.jl b/src/solution/fix.jl index a8b4877..f9deacb 100644 --- a/src/solution/fix.jl +++ b/src/solution/fix.jl @@ -16,14 +16,14 @@ function fix!(model::JuMP.Model, solution::AbstractDict)::Nothing for g in instance.units for t in 1:T is_on_value = round(solution["Is on"][g.name][t]) - production_value = + prod_value = round(solution["Production (MW)"][g.name][t], digits = 5) reserve_value = round(solution["Reserve (MW)"][g.name][t], digits = 5) JuMP.fix(is_on[g.name, t], is_on_value, force = true) JuMP.fix( prod_above[g.name, t], - production_value - is_on_value * g.min_power[t], + prod_value - is_on_value * g.min_power[t], force = true, ) JuMP.fix(reserve[g.name, t], reserve_value, force = true) diff --git a/src/solution/methods/XaQiWaTh19/structs.jl b/src/solution/methods/XaQiWaTh19/structs.jl index feeb6d4..15f688b 100644 --- a/src/solution/methods/XaQiWaTh19/structs.jl +++ b/src/solution/methods/XaQiWaTh19/structs.jl @@ -9,6 +9,8 @@ import DataStructures: PriorityQueue time_limit::Float64 gap_limit::Float64 two_phase_gap::Bool + max_violations_per_line::Int + max_violations_per_period::Int end 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. IEEE Transactions on Power Systems, 34(3), 2457-2460. -Fields -========= +## Fields + - `time_limit`: the time limit over the entire optimization procedure. - `gap_limit`: @@ -26,6 +28,13 @@ Fields - `two_phase_gap`: if true, solve the problem with large gap tolerance first, then reduce 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 time_limit::Float64 @@ -35,11 +44,11 @@ struct _XaQiWaTh19 max_violations_per_period::Int function _XaQiWaTh19(; - time_limit::Float64, - gap_limit::Float64, - two_phase_gap::Bool, - max_violations_per_line::Int, - max_violations_per_period::Int, + time_limit::Float64 = 86400.0, + gap_limit::Float64 = 1e-3, + two_phase_gap::Bool = true, + max_violations_per_line::Int = 1, + max_violations_per_period::Int = 5, ) return new( time_limit, diff --git a/src/solution/optimize.jl b/src/solution/optimize.jl index 6f65fc7..ef0e2f7 100644 --- a/src/solution/optimize.jl +++ b/src/solution/optimize.jl @@ -10,14 +10,5 @@ advanced methods to accelerate the solution process and to enforce transmission and N-1 security constraints. """ function optimize!(model::JuMP.Model)::Nothing - return UnitCommitment.optimize!( - model, - _XaQiWaTh19( - time_limit = 3600.0, - gap_limit = 1e-4, - two_phase_gap = true, - max_violations_per_line = 1, - max_violations_per_period = 5, - ), - ) + return UnitCommitment.optimize!(model, _XaQiWaTh19()) end diff --git a/src/solution/warmstart.jl b/src/solution/warmstart.jl index 5f09352..678f250 100644 --- a/src/solution/warmstart.jl +++ b/src/solution/warmstart.jl @@ -5,8 +5,6 @@ function set_warm_start!(model::JuMP.Model, solution::AbstractDict)::Nothing instance, T = model[:instance], model[:instance].time is_on = model[:is_on] - prod_above = model[:prod_above] - reserve = model[:reserve] for g in instance.units for t in 1:T JuMP.set_start_value(is_on[g.name, t], solution["Is on"][g.name][t]) diff --git a/src/transforms/initcond.jl b/src/transform/initcond.jl similarity index 100% rename from src/transforms/initcond.jl rename to src/transform/initcond.jl diff --git a/src/transforms/slice.jl b/src/transform/slice.jl similarity index 100% rename from src/transforms/slice.jl rename to src/transform/slice.jl