From 9224cd2efbb553ecbc598486827f0e0198975731 Mon Sep 17 00:00:00 2001 From: Alinson S Xavier Date: Thu, 27 May 2021 21:37:38 -0500 Subject: [PATCH] Format source code with JuliaFormatter; set up GH Actions --- .JuliaFormatter.toml | 5 + .github/workflows/lint.yml | 28 ++ Makefile | 4 + benchmark/run.jl | 19 +- src/UnitCommitment.jl | 16 +- src/convert.jl | 12 +- src/initcond.jl | 36 ++- src/instance.jl | 142 +++++----- src/log.jl | 44 +-- src/model.jl | 538 ++++++++++++++++++++----------------- src/screening.jl | 115 ++++---- src/sensitivity.jl | 32 ++- src/sysimage.jl | 12 +- src/validate.jl | 230 ++++++++++------ test/convert_test.jl | 9 +- test/initcond_test.jl | 8 +- test/instance_test.jl | 122 ++++----- test/model_test.jl | 8 +- test/screening_test.jl | 90 +++---- test/sensitivity_test.jl | 182 +++++++------ test/validate_test.jl | 26 +- 21 files changed, 919 insertions(+), 759 deletions(-) create mode 100644 .JuliaFormatter.toml create mode 100644 .github/workflows/lint.yml diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 0000000..0f4f1be --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1,5 @@ +always_for_in = true +always_use_return = true +margin = 80 +remove_extra_newlines = true +short_to_long_function_def = true diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml new file mode 100644 index 0000000..aeff998 --- /dev/null +++ b/.github/workflows/lint.yml @@ -0,0 +1,28 @@ +name: lint +on: + push: + pull_request: +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: julia-actions/setup-julia@latest + with: + version: '1' + - uses: actions/checkout@v1 + - name: Format check + shell: julia --color=yes {0} + run: | + using Pkg + Pkg.add(PackageSpec(name="JuliaFormatter", version="0.14.4")) + using JuliaFormatter + format("src", verbose=true) + format("test", verbose=true) + format("benchmark", verbose=true) + out = String(read(Cmd(`git diff`))) + if isempty(out) + exit(0) + end + @error "Some files have not been formatted !!!" + write(stdout, out) + exit(1) diff --git a/Makefile b/Makefile index f08c844..b4e3f08 100644 --- a/Makefile +++ b/Makefile @@ -22,4 +22,8 @@ test: build/sysimage.so @echo Running tests... $(JULIA) --sysimage build/sysimage.so -e 'using Pkg; Pkg.test("UnitCommitment")' | tee build/test.log + +format: + julia -e 'using JuliaFormatter; format("src"); format("test"); format("benchmark")' + .PHONY: docs test diff --git a/benchmark/run.jl b/benchmark/run.jl index 2af77bf..60c35c0 100644 --- a/benchmark/run.jl +++ b/benchmark/run.jl @@ -30,34 +30,37 @@ function main() time_model = @elapsed begin model = build_model( - instance=instance, - optimizer=optimizer_with_attributes( + instance = instance, + optimizer = optimizer_with_attributes( Gurobi.Optimizer, "Threads" => 4, "Seed" => rand(1:1000), ), - variable_names=true, + variable_names = true, ) end @info "Optimizing..." BLAS.set_num_threads(1) - UnitCommitment.optimize!(model, time_limit=time_limit, gap_limit=1e-3) - + UnitCommitment.optimize!( + model, + time_limit = time_limit, + gap_limit = 1e-3, + ) end @info @sprintf("Total time was %.2f seconds", total_time) @info "Writing: $solution_filename" solution = UnitCommitment.solution(model) open(solution_filename, "w") do file - JSON.print(file, solution, 2) + return JSON.print(file, solution, 2) end @info "Verifying solution..." - UnitCommitment.validate(instance, solution) + UnitCommitment.validate(instance, solution) @info "Exporting model..." - JuMP.write_to_file(model, model_filename) + return JuMP.write_to_file(model, model_filename) end main() diff --git a/src/UnitCommitment.jl b/src/UnitCommitment.jl index cda6e23..4a36550 100644 --- a/src/UnitCommitment.jl +++ b/src/UnitCommitment.jl @@ -3,12 +3,12 @@ # Released under the modified BSD license. See COPYING.md for more details. module UnitCommitment - include("log.jl") - include("instance.jl") - include("screening.jl") - include("model.jl") - include("sensitivity.jl") - include("validate.jl") - include("convert.jl") - include("initcond.jl") +include("log.jl") +include("instance.jl") +include("screening.jl") +include("model.jl") +include("sensitivity.jl") +include("validate.jl") +include("convert.jl") +include("initcond.jl") end diff --git a/src/convert.jl b/src/convert.jl index b7b47fd..7599061 100644 --- a/src/convert.jl +++ b/src/convert.jl @@ -10,20 +10,20 @@ function _read_json(path::String)::OrderedDict else file = open(path) end - return JSON.parse(file, dicttype=()->DefaultOrderedDict(nothing)) + return JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)) end function _read_egret_solution(path::String)::OrderedDict egret = _read_json(path) T = length(egret["system"]["time_keys"]) - + solution = OrderedDict() is_on = solution["Is on"] = OrderedDict() production = solution["Production (MW)"] = OrderedDict() reserve = solution["Reserve (MW)"] = OrderedDict() production_cost = solution["Production cost (\$)"] = OrderedDict() startup_cost = solution["Startup cost (\$)"] = OrderedDict() - + for (gen_name, gen_dict) in egret["elements"]["generator"] if endswith(gen_name, "_T") || endswith(gen_name, "_R") gen_name = gen_name[1:end-2] @@ -39,18 +39,18 @@ function _read_egret_solution(path::String)::OrderedDict else reserve[gen_name] = zeros(T) end - startup_cost[gen_name] = zeros(T) + startup_cost[gen_name] = zeros(T) production_cost[gen_name] = zeros(T) if "commitment_cost" in keys(gen_dict) for t in 1:T x = gen_dict["commitment"]["values"][t] commitment_cost = gen_dict["commitment_cost"]["values"][t] prod_above_cost = gen_dict["production_cost"]["values"][t] - prod_base_cost = gen_dict["p_cost"]["values"][1][2] * x + prod_base_cost = gen_dict["p_cost"]["values"][1][2] * x startup_cost[gen_name][t] = commitment_cost - prod_base_cost production_cost[gen_name][t] = prod_above_cost + prod_base_cost end end end return solution -end \ No newline at end of file +end diff --git a/src/initcond.jl b/src/initcond.jl index b1be33a..6bbf0ea 100644 --- a/src/initcond.jl +++ b/src/initcond.jl @@ -19,33 +19,31 @@ function generate_initial_conditions!( B = instance.buses t = 1 mip = JuMP.Model(optimizer) - + # Decision variables @variable(mip, x[G], Bin) @variable(mip, p[G] >= 0) - + # Constraint: Minimum power - @constraint(mip, - min_power[g in G], - p[g] >= g.min_power[t] * x[g]) - + @constraint(mip, min_power[g in G], p[g] >= g.min_power[t] * x[g]) + # Constraint: Maximum power - @constraint(mip, - max_power[g in G], - p[g] <= g.max_power[t] * x[g]) - + @constraint(mip, max_power[g in G], p[g] <= g.max_power[t] * x[g]) + # Constraint: Production equals demand - @constraint(mip, - power_balance, - sum(b.load[t] for b in B) == sum(p[g] for g in G)) - + @constraint( + mip, + power_balance, + sum(b.load[t] for b in B) == sum(p[g] for g in G) + ) + # Constraint: Must run for g in G if g.must_run[t] @constraint(mip, x[g] == 1) end end - + # Objective function function cost_slope(g) mw = g.min_power[t] @@ -60,12 +58,10 @@ function generate_initial_conditions!( return c / mw end end - @objective(mip, - Min, - sum(p[g] * cost_slope(g) for g in G)) - + @objective(mip, Min, sum(p[g] * cost_slope(g) for g in G)) + JuMP.optimize!(mip) - + for g in G if JuMP.value(x[g]) > 0 g.initial_power = JuMP.value(p[g]) diff --git a/src/instance.jl b/src/instance.jl index 8cb9491..76aa049 100644 --- a/src/instance.jl +++ b/src/instance.jl @@ -8,7 +8,6 @@ using DataStructures using GZip import Base: getindex, time - mutable struct Bus name::String offset::Int @@ -17,19 +16,16 @@ mutable struct Bus price_sensitive_loads::Vector end - mutable struct CostSegment mw::Vector{Float64} cost::Vector{Float64} end - mutable struct StartupCategory delay::Int cost::Float64 end - mutable struct Unit name::String bus::Bus @@ -50,7 +46,6 @@ mutable struct Unit startup_categories::Vector{StartupCategory} end - mutable struct TransmissionLine name::String offset::Int @@ -63,19 +58,16 @@ mutable struct TransmissionLine flow_limit_penalty::Vector{Float64} end - mutable struct Reserves spinning::Vector{Float64} end - mutable struct Contingency name::String lines::Vector{TransmissionLine} units::Vector{Unit} end - mutable struct PriceSensitiveLoad name::String bus::Bus @@ -83,7 +75,6 @@ mutable struct PriceSensitiveLoad revenue::Vector{Float64} end - mutable struct UnitCommitmentInstance time::Int power_balance_penalty::Vector{Float64} @@ -95,25 +86,26 @@ mutable struct UnitCommitmentInstance price_sensitive_loads::Vector{PriceSensitiveLoad} 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.price_sensitive_loads)) price sensitive loads, ", + ) print(io, "$(instance.time) time steps") print(io, ")") + return end - -function read_benchmark(name::AbstractString) :: UnitCommitmentInstance +function read_benchmark(name::AbstractString)::UnitCommitmentInstance basedir = dirname(@__FILE__) return UnitCommitment.read("$basedir/../instances/$name.json.gz") end - function read(path::AbstractString)::UnitCommitmentInstance if endswith(path, ".gz") return _read(gzopen(path)) @@ -122,50 +114,51 @@ function read(path::AbstractString)::UnitCommitmentInstance end end - function _read(file::IO)::UnitCommitmentInstance - return _from_json(JSON.parse(file, dicttype=()->DefaultOrderedDict(nothing))) + return _from_json( + JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)), + ) end - -function _from_json(json; repair=true) +function _from_json(json; repair = true) units = Unit[] buses = Bus[] contingencies = Contingency[] lines = TransmissionLine[] loads = PriceSensitiveLoad[] - function scalar(x; default=nothing) + function scalar(x; default = nothing) x !== nothing || return default - x + return x end - + time_horizon = json["Parameters"]["Time (h)"] if time_horizon === nothing time_horizon = json["Parameters"]["Time horizon (h)"] end - time_horizon !== nothing || error("Missing required parameter: Time horizon (h)") - time_step = scalar(json["Parameters"]["Time step (min)"], default=60) - (60 % time_step == 0) || error("Time step $time_step is not a divisor of 60") + time_horizon !== nothing || error("Missing parameter: Time horizon (h)") + time_step = scalar(json["Parameters"]["Time step (min)"], default = 60) + (60 % time_step == 0) || + error("Time step $time_step is not a divisor of 60") time_multiplier = 60 ÷ time_step T = time_horizon * time_multiplier - - name_to_bus = Dict{String, Bus}() - name_to_line = Dict{String, TransmissionLine}() - name_to_unit = Dict{String, Unit}() - - function timeseries(x; default=nothing) + + name_to_bus = Dict{String,Bus}() + name_to_line = Dict{String,TransmissionLine}() + name_to_unit = Dict{String,Unit}() + + function timeseries(x; default = nothing) x !== nothing || return default x isa Array || return [x for t in 1:T] return x end - + # Read parameters power_balance_penalty = timeseries( json["Parameters"]["Power balance penalty (\$/MW)"], - default=[1000.0 for t in 1:T], + default = [1000.0 for t in 1:T], ) - + # Read buses for (bus_name, dict) in json["Buses"] bus = Bus( @@ -178,15 +171,19 @@ function _from_json(json; repair=true) name_to_bus[bus_name] = bus push!(buses, bus) end - + # Read units for (unit_name, dict) in json["Generators"] bus = name_to_bus[dict["Bus"]] - + # Read production cost curve K = length(dict["Production cost curve (MW)"]) - curve_mw = hcat([timeseries(dict["Production cost curve (MW)"][k]) for k in 1:K]...) - curve_cost = hcat([timeseries(dict["Production cost curve (\$)"][k]) for k in 1:K]...) + curve_mw = hcat( + [timeseries(dict["Production cost curve (MW)"][k]) for k in 1:K]..., + ) + curve_cost = hcat( + [timeseries(dict["Production cost curve (\$)"][k]) for k in 1:K]..., + ) min_power = curve_mw[:, 1] max_power = curve_mw[:, K] min_power_cost = curve_cost[:, 1] @@ -194,13 +191,13 @@ function _from_json(json; repair=true) for k in 2:K amount = curve_mw[:, k] - curve_mw[:, k-1] cost = (curve_cost[:, k] - curve_cost[:, k-1]) ./ amount - replace!(cost, NaN=>0.0) + replace!(cost, NaN => 0.0) push!(segments, CostSegment(amount, cost)) end - + # Read startup costs - startup_delays = scalar(dict["Startup delays (h)"], default=[1]) - startup_costs = scalar(dict["Startup costs (\$)"], default=[0.]) + startup_delays = scalar(dict["Startup delays (h)"], default = [1]) + startup_costs = scalar(dict["Startup costs (\$)"], default = [0.0]) startup_categories = StartupCategory[] for k in 1:length(startup_delays) push!( @@ -211,40 +208,43 @@ function _from_json(json; repair=true) ), ) end - + # Read and validate initial conditions - initial_power = scalar(dict["Initial power (MW)"], default=nothing) - initial_status = scalar(dict["Initial status (h)"], default=nothing) + initial_power = scalar(dict["Initial power (MW)"], default = nothing) + initial_status = scalar(dict["Initial status (h)"], default = nothing) if initial_power === nothing - initial_status === nothing || error("unit $unit_name has initial status but no initial power") + initial_status === nothing || + error("unit $unit_name has initial status but no initial power") else - initial_status !== nothing || error("unit $unit_name has initial power but no initial status") - initial_status != 0 || error("unit $unit_name has invalid initial status") + initial_status !== nothing || + error("unit $unit_name has initial power but no initial status") + initial_status != 0 || + error("unit $unit_name has invalid initial status") if initial_status < 0 && initial_power > 1e-3 error("unit $unit_name has invalid initial power") end initial_status *= time_multiplier end - + unit = Unit( unit_name, bus, max_power, min_power, - timeseries(dict["Must run?"], default=[false for t in 1:T]), + timeseries(dict["Must run?"], default = [false for t in 1:T]), min_power_cost, segments, - scalar(dict["Minimum uptime (h)"], default=1) * time_multiplier, - scalar(dict["Minimum downtime (h)"], default=1) * time_multiplier, - scalar(dict["Ramp up limit (MW)"], default=1e6), - scalar(dict["Ramp down limit (MW)"], default=1e6), - scalar(dict["Startup limit (MW)"], default=1e6), - scalar(dict["Shutdown limit (MW)"], default=1e6), + scalar(dict["Minimum uptime (h)"], default = 1) * time_multiplier, + scalar(dict["Minimum downtime (h)"], default = 1) * time_multiplier, + scalar(dict["Ramp up limit (MW)"], default = 1e6), + scalar(dict["Ramp down limit (MW)"], default = 1e6), + scalar(dict["Startup limit (MW)"], default = 1e6), + scalar(dict["Shutdown limit (MW)"], default = 1e6), initial_status, initial_power, timeseries( dict["Provides spinning reserves?"], - default=[true for t in 1:T], + default = [true for t in 1:T], ), startup_categories, ) @@ -252,16 +252,14 @@ function _from_json(json; repair=true) name_to_unit[unit_name] = unit push!(units, unit) end - + # Read reserves reserves = Reserves(zeros(T)) if "Reserves" in keys(json) - reserves.spinning = timeseries( - json["Reserves"]["Spinning (MW)"], - default=zeros(T), - ) + reserves.spinning = + timeseries(json["Reserves"]["Spinning (MW)"], default = zeros(T)) end - + # Read transmission lines if "Transmission lines" in keys(json) for (line_name, dict) in json["Transmission lines"] @@ -274,38 +272,40 @@ function _from_json(json; repair=true) scalar(dict["Susceptance (S)"]), timeseries( dict["Normal flow limit (MW)"], - default=[1e8 for t in 1:T], + default = [1e8 for t in 1:T], ), timeseries( dict["Emergency flow limit (MW)"], - default=[1e8 for t in 1:T], + default = [1e8 for t in 1:T], ), timeseries( dict["Flow limit penalty (\$/MW)"], - default=[5000.0 for t in 1:T], + default = [5000.0 for t in 1:T], ), ) name_to_line[line_name] = line push!(lines, line) end end - + # Read contingencies if "Contingencies" in keys(json) for (cont_name, dict) in json["Contingencies"] affected_units = Unit[] affected_lines = TransmissionLine[] if "Affected lines" in keys(dict) - affected_lines = [name_to_line[l] for l in dict["Affected lines"]] + affected_lines = + [name_to_line[l] for l in dict["Affected lines"]] end if "Affected units" in keys(dict) - affected_units = [name_to_unit[u] for u in dict["Affected units"]] + affected_units = + [name_to_unit[u] for u in dict["Affected units"]] end cont = Contingency(cont_name, affected_lines, affected_units) push!(contingencies, cont) end end - + # Read price-sensitive loads if "Price-sensitive loads" in keys(json) for (load_name, dict) in json["Price-sensitive loads"] @@ -320,7 +320,7 @@ function _from_json(json; repair=true) push!(loads, load) end end - + instance = UnitCommitmentInstance( T, power_balance_penalty, @@ -337,7 +337,6 @@ function _from_json(json; repair=true) return instance end - """ slice(instance, range) @@ -387,5 +386,4 @@ function slice( return modified end - export UnitCommitmentInstance diff --git a/src/log.jl b/src/log.jl index 8d986ab..e3daf87 100644 --- a/src/log.jl +++ b/src/log.jl @@ -7,35 +7,37 @@ using Base.CoreLogging, Logging, Printf struct TimeLogger <: AbstractLogger initial_time::Float64 - file::Union{Nothing, IOStream} - screen_log_level - io_log_level + file::Union{Nothing,IOStream} + screen_log_level::Any + io_log_level::Any end function TimeLogger(; - initial_time::Float64, - file::Union{Nothing, IOStream} = nothing, - screen_log_level = CoreLogging.Info, - io_log_level = CoreLogging.Info, - ) :: TimeLogger + initial_time::Float64, + file::Union{Nothing,IOStream} = nothing, + screen_log_level = CoreLogging.Info, + io_log_level = CoreLogging.Info, +)::TimeLogger return TimeLogger(initial_time, file, screen_log_level, io_log_level) end min_enabled_level(logger::TimeLogger) = logger.io_log_level shouldlog(logger::TimeLogger, level, _module, group, id) = true -function handle_message(logger::TimeLogger, - level, - message, - _module, - group, - id, - filepath, - line; - kwargs...) +function handle_message( + logger::TimeLogger, + level, + message, + _module, + group, + id, + filepath, + line; + kwargs..., +) elapsed_time = time() - logger.initial_time time_string = @sprintf("[%12.3f] ", elapsed_time) - + if level >= Logging.Error color = :light_red elseif level >= Logging.Warn @@ -43,9 +45,9 @@ function handle_message(logger::TimeLogger, else color = :light_green end - + if level >= logger.screen_log_level - printstyled(time_string, color=color) + printstyled(time_string, color = color) println(message) end if logger.file !== nothing && level >= logger.io_log_level @@ -58,5 +60,5 @@ end function _setup_logger() initial_time = time() - global_logger(TimeLogger(initial_time=initial_time)) + return global_logger(TimeLogger(initial_time = initial_time)) end diff --git a/src/model.jl b/src/model.jl index 1e546c7..acd1e1e 100644 --- a/src/model.jl +++ b/src/model.jl @@ -5,15 +5,14 @@ using JuMP, MathOptInterface, DataStructures import JuMP: value, fix, set_name - # Extend some JuMP functions so that decision variables can be safely replaced by # (constant) floating point numbers. function value(x::Float64) - x + return x end function fix(x::Float64, v::Float64; force) - abs(x - v) < 1e-6 || error("Value mismatch: $x != $v") + return abs(x - v) < 1e-6 || error("Value mismatch: $x != $v") end function set_name(x::Float64, n::String) @@ -21,20 +20,19 @@ function set_name(x::Float64, n::String) end function build_model(; - filename::Union{String, Nothing}=nothing, - instance::Union{UnitCommitmentInstance, Nothing}=nothing, - isf::Union{Matrix{Float64}, Nothing}=nothing, - lodf::Union{Matrix{Float64}, Nothing}=nothing, - isf_cutoff::Float64=0.005, - lodf_cutoff::Float64=0.001, - optimizer=nothing, - variable_names::Bool=false, + filename::Union{String,Nothing} = nothing, + instance::Union{UnitCommitmentInstance,Nothing} = nothing, + isf::Union{Matrix{Float64},Nothing} = nothing, + lodf::Union{Matrix{Float64},Nothing} = nothing, + isf_cutoff::Float64 = 0.005, + lodf_cutoff::Float64 = 0.001, + optimizer = nothing, + variable_names::Bool = false, )::JuMP.Model - if (filename === nothing) && (instance === nothing) error("Either filename or instance must be specified") end - + if filename !== nothing @info "Reading: $(filename)" time_read = @elapsed begin @@ -42,7 +40,7 @@ function build_model(; end @info @sprintf("Read problem in %.2f seconds", time_read) end - + if length(instance.buses) == 1 isf = zeros(0, 0) lodf = zeros(0, 0) @@ -51,25 +49,29 @@ function build_model(; @info "Computing injection shift factors..." time_isf = @elapsed begin isf = UnitCommitment._injection_shift_factors( - lines=instance.lines, - buses=instance.buses, + lines = instance.lines, + buses = instance.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, - isf=isf, + lines = instance.lines, + buses = instance.buses, + isf = isf, ) end @info @sprintf("Computed LODF in %.2f seconds", time_lodf) - - @info @sprintf("Applying PTDF and LODF cutoffs (%.5f, %.5f)", isf_cutoff, lodf_cutoff) - isf[abs.(isf) .< isf_cutoff] .= 0 - lodf[abs.(lodf) .< lodf_cutoff] .= 0 + + @info @sprintf( + "Applying PTDF and LODF cutoffs (%.5f, %.5f)", + isf_cutoff, + lodf_cutoff + ) + isf[abs.(isf).= 1)) + 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], @@ -267,142 +280,171 @@ function _add_unit!(model::JuMP.Model, g::Unit) 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]) + 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]) + 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)) + 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]) + 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]) + 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]) + 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) + 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) + 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) + 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) + 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]) - + 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] = + 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]) + 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]) - + 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]) - + 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) + 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) + 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_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] @@ -412,45 +454,36 @@ function _build_net_injection_eqs!(model::JuMP.Model) @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 - ) + 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] - ) + 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 _enforce_transmission( - ; +function _enforce_transmission(; model::JuMP.Model, violation::Violation, isf::Matrix{Float64}, lodf::Matrix{Float64}, )::Nothing - instance = model[:instance] + instance = model[:instance] limit::Float64 = 0.0 overflow = model[:overflow] net_injection = model[:net_injection] - + if violation.outage_line === nothing limit = violation.monitored_line.normal_flow_limit[violation.time] @info @sprintf( @@ -469,34 +502,42 @@ function _enforce_transmission( violation.outage_line.name, ) end - + fm = violation.monitored_line.name t = violation.time - flow = @variable(model, base_name="flow[$fm,$t]") - + flow = @variable(model, base_name = "flow[$fm,$t]") + v = overflow[violation.monitored_line.name, violation.time] - @constraint(model, flow <= limit + v) + @constraint(model, flow <= limit + v) @constraint(model, -flow <= limit + v) - + if violation.outage_line === nothing - @constraint(model, flow == sum(net_injection[b.name, violation.time] * - isf[violation.monitored_line.offset, b.offset] - for b in instance.buses - if b.offset > 0)) + @constraint( + model, + flow == sum( + net_injection[b.name, violation.time] * + isf[violation.monitored_line.offset, b.offset] for + b in instance.buses if b.offset > 0 + ) + ) else - @constraint(model, flow == sum(net_injection[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)) - end - nothing + @constraint( + model, + flow == sum( + net_injection[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 + ) + ) + end + return nothing end - function _set_names!(model::JuMP.Model) @info "Setting variable and constraint names..." time_varnames = @elapsed begin @@ -505,7 +546,6 @@ function _set_names!(model::JuMP.Model) @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 @@ -519,72 +559,75 @@ function _set_names!(dict::Dict) end end - function solution(model::JuMP.Model) instance, T = model[:instance], model[:instance].time function timeseries(vars, collection) return OrderedDict( - b.name => [round(value(vars[b.name, t]), digits=5) for t in 1:T] + b.name => [round(value(vars[b.name, t]), digits = 5) for t in 1:T] for b in collection ) end function production_cost(g) return [ - value(model[:is_on][g.name, t]) * g.min_power_cost[t] + - sum( + value(model[:is_on][g.name, t]) * g.min_power_cost[t] + sum( Float64[ - value(model[:segprod][g.name, t, k]) * g.cost_segments[k].cost[t] - for k in 1:length(g.cost_segments) - ] - ) - for t in 1:T + value(model[:segprod][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) return [ - value(model[:is_on][g.name, t]) * g.min_power[t] + - sum( + value(model[:is_on][g.name, t]) * g.min_power[t] + sum( Float64[ - value(model[:segprod][g.name, t, k]) - for k in 1:length(g.cost_segments) - ] - ) - for t in 1:T + value(model[:segprod][g.name, t, k]) for + k in 1:length(g.cost_segments) + ], + ) for t in 1:T ] end function startup_cost(g) S = length(g.startup_categories) - return [sum(g.startup_categories[s].cost * value(model[:startup][g.name, t, s]) - for s in 1:S) - for t in 1:T] + return [ + sum( + g.startup_categories[s].cost * + value(model[:startup][g.name, t, s]) for s in 1:S + ) for t in 1:T + ] 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["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["Reserve (MW)"] = timeseries(model[:reserve], instance.units) - sol["Net injection (MW)"] = timeseries(model[:net_injection], instance.buses) + 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) end if !isempty(instance.price_sensitive_loads) - sol["Price-sensitive loads (MW)"] = timeseries(model[:loads], instance.price_sensitive_loads) + sol["Price-sensitive loads (MW)"] = + timeseries(model[:loads], instance.price_sensitive_loads) end return sol end - function write(filename::AbstractString, solution::AbstractDict)::Nothing open(filename, "w") do file - JSON.print(file, solution, 2) - end + return JSON.print(file, solution, 2) + end + return end - function fix!(model::JuMP.Model, solution::AbstractDict)::Nothing instance, T = model[:instance], model[:instance].time is_on = model[:is_on] @@ -593,20 +636,22 @@ 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 = 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) + production_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], - force=true, + force = true, ) - JuMP.fix(reserve[g.name, t], reserve_value, force=true) + JuMP.fix(reserve[g.name, t], reserve_value, force = true) end end + return end - function set_warm_start!(model::JuMP.Model, solution::AbstractDict)::Nothing instance, T = model[:instance], model[:instance].time is_on = model[:is_on] @@ -615,20 +660,25 @@ function set_warm_start!(model::JuMP.Model, solution::AbstractDict)::Nothing 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]) - JuMP.set_start_value(switch_on[g.name, t], solution["Switch on"][g.name][t]) - JuMP.set_start_value(switch_off[g.name, t], solution["Switch off"][g.name][t]) + JuMP.set_start_value( + switch_on[g.name, t], + solution["Switch on"][g.name][t], + ) + JuMP.set_start_value( + switch_off[g.name, t], + solution["Switch off"][g.name][t], + ) end end + return end - function optimize!( model::JuMP.Model; - time_limit=3600, - gap_limit=1e-4, - two_phase_gap=true, + time_limit = 3600, + gap_limit = 1e-4, + two_phase_gap = true, )::Nothing - function set_gap(gap) try JuMP.set_optimizer_attribute(model, "MIPGap", gap) @@ -637,20 +687,20 @@ function optimize!( @warn "Could not change MIP gap tolerance" end end - + instance = model[:instance] initial_time = time() - + large_gap = false has_transmission = (length(model[:isf]) > 0) - + if has_transmission && two_phase_gap set_gap(1e-2) large_gap = true else set_gap(gap_limit) end - + while true time_elapsed = time() - initial_time time_remaining = time_limit - time_elapsed @@ -658,18 +708,21 @@ function optimize!( @info "Time limit exceeded" break end - - @info @sprintf("Setting MILP time limit to %.2f seconds", time_remaining) + + @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) if isempty(violations) - @info "No violations found" + @info "No violations found" if large_gap large_gap = false set_gap(gap_limit) @@ -680,10 +733,9 @@ function optimize!( _enforce_transmission(model, violations) end end - - nothing -end + return +end function _find_violations(model::JuMP.Model) instance = model[:instance] @@ -695,36 +747,38 @@ function _find_violations(model::JuMP.Model) time_screening = @elapsed begin non_slack_buses = [b for b in instance.buses if b.offset > 0] net_injection_values = [ - value(net_injection[b.name, t]) - for b in non_slack_buses, t in 1:instance.time + value(net_injection[b.name, t]) for b in non_slack_buses, + t in 1:instance.time ] overflow_values = [ - value(overflow[lm.name, t]) - for lm in instance.lines, t in 1:instance.time + value(overflow[lm.name, t]) for lm in instance.lines, + t in 1:instance.time ] violations = UnitCommitment._find_violations( - instance=instance, - net_injections=net_injection_values, - overflow=overflow_values, - isf=model[:isf], - lodf=model[:lodf], + instance = instance, + net_injections = net_injection_values, + overflow = overflow_values, + isf = model[:isf], + lodf = model[:lodf], ) end - @info @sprintf("Verified transmission limits in %.2f seconds", time_screening) + @info @sprintf( + "Verified transmission limits in %.2f seconds", + time_screening + ) return violations end - function _enforce_transmission( model::JuMP.Model, violations::Vector{Violation}, )::Nothing for v in violations _enforce_transmission( - model=model, - violation=v, - isf=model[:isf], - lodf=model[:lodf], + model = model, + violation = v, + isf = model[:isf], + lodf = model[:lodf], ) end return diff --git a/src/screening.jl b/src/screening.jl index 872f7b2..9a4971b 100644 --- a/src/screening.jl +++ b/src/screening.jl @@ -4,49 +4,44 @@ # Copyright (C) 2019 Argonne National Laboratory # Written by Alinson Santos Xavier - using DataStructures using Base.Threads - struct Violation time::Int monitored_line::TransmissionLine - outage_line::Union{TransmissionLine, Nothing} + outage_line::Union{TransmissionLine,Nothing} amount::Float64 # Violation amount (in MW) end - function Violation(; - time::Int, - monitored_line::TransmissionLine, - outage_line::Union{TransmissionLine, Nothing}, - amount::Float64, - )::Violation + time::Int, + monitored_line::TransmissionLine, + outage_line::Union{TransmissionLine,Nothing}, + amount::Float64, +)::Violation return Violation(time, monitored_line, outage_line, amount) end - mutable struct ViolationFilter max_per_line::Int max_total::Int - queues::Dict{Int, PriorityQueue{Violation, Float64}} + queues::Dict{Int,PriorityQueue{Violation,Float64}} end - function ViolationFilter(; - max_per_line::Int=1, - max_total::Int=5, - )::ViolationFilter + max_per_line::Int = 1, + max_total::Int = 5, +)::ViolationFilter return ViolationFilter(max_per_line, max_total, Dict()) end - function _offer(filter::ViolationFilter, v::Violation)::Nothing if v.monitored_line.offset ∉ keys(filter.queues) - filter.queues[v.monitored_line.offset] = PriorityQueue{Violation, Float64}() + filter.queues[v.monitored_line.offset] = + PriorityQueue{Violation,Float64}() end - q::PriorityQueue{Violation, Float64} = filter.queues[v.monitored_line.offset] + q::PriorityQueue{Violation,Float64} = filter.queues[v.monitored_line.offset] if length(q) < filter.max_per_line enqueue!(q, v => v.amount) else @@ -55,13 +50,12 @@ function _offer(filter::ViolationFilter, v::Violation)::Nothing enqueue!(q, v => v.amount) end end - nothing + return nothing end - -function _query(filter::ViolationFilter)::Array{Violation, 1} +function _query(filter::ViolationFilter)::Array{Violation,1} violations = Array{Violation,1}() - time_queue = PriorityQueue{Violation, Float64}() + time_queue = PriorityQueue{Violation,Float64}() for l in keys(filter.queues) line_queue = filter.queues[l] while length(line_queue) > 0 @@ -82,7 +76,6 @@ function _query(filter::ViolationFilter)::Array{Violation, 1} return violations end - """ function _find_violations( @@ -104,49 +97,46 @@ UnitCommitment.line_outage_factors. The argument `overflow` specifies how much flow above the transmission limits (in MW) is allowed. It should be an L x T matrix, where L is the number of transmission lines. """ -function _find_violations( - ; +function _find_violations(; instance::UnitCommitmentInstance, - net_injections::Array{Float64, 2}, - overflow::Array{Float64, 2}, + net_injections::Array{Float64,2}, + overflow::Array{Float64,2}, isf::Array{Float64,2}, lodf::Array{Float64,2}, max_per_line::Int = 1, max_per_period::Int = 5, -)::Array{Violation, 1} - +)::Array{Violation,1} B = length(instance.buses) - 1 L = length(instance.lines) T = instance.time K = nthreads() - + size(net_injections) == (B, T) || error("net_injections has incorrect size") size(isf) == (L, B) || error("isf has incorrect size") size(lodf) == (L, L) || error("lodf has incorrect size") - + filters = Dict( t => ViolationFilter( - max_total=max_per_period, - max_per_line=max_per_line, - ) - for t in 1:T + max_total = max_per_period, + max_per_line = max_per_line, + ) for t in 1:T ) - + pre_flow::Array{Float64} = zeros(L, K) # pre_flow[lm, thread] post_flow::Array{Float64} = zeros(L, L, K) # post_flow[lm, lc, thread] pre_v::Array{Float64} = zeros(L, K) # pre_v[lm, thread] post_v::Array{Float64} = zeros(L, L, K) # post_v[lm, lc, thread] - + normal_limits::Array{Float64,2} = [ - l.normal_flow_limit[t] + overflow[l.offset, t] - for l in instance.lines, t in 1:T + l.normal_flow_limit[t] + overflow[l.offset, t] for + l in instance.lines, t in 1:T ] - + emergency_limits::Array{Float64,2} = [ - l.emergency_flow_limit[t] + overflow[l.offset, t] - for l in instance.lines, t in 1:T + l.emergency_flow_limit[t] + overflow[l.offset, t] for + l in instance.lines, t in 1:T ] - + is_vulnerable::Array{Bool} = zeros(Bool, L) for c in instance.contingencies is_vulnerable[c.lines[1].offset] = true @@ -154,68 +144,69 @@ function _find_violations( @threads for t in 1:T k = threadid() - + # Pre-contingency flows pre_flow[:, k] = isf * net_injections[:, t] - + # Post-contingency flows for lc in 1:L, lm in 1:L - post_flow[lm, lc, k] = pre_flow[lm, k] + pre_flow[lc, k] * lodf[lm, lc] + post_flow[lm, lc, k] = + pre_flow[lm, k] + pre_flow[lc, k] * lodf[lm, lc] end - + # Pre-contingency violations for lm in 1:L pre_v[lm, k] = max( 0.0, pre_flow[lm, k] - normal_limits[lm, t], - - pre_flow[lm, k] - normal_limits[lm, t], + -pre_flow[lm, k] - normal_limits[lm, t], ) end - + # Post-contingency violations for lc in 1:L, lm in 1:L post_v[lm, lc, k] = max( 0.0, post_flow[lm, lc, k] - emergency_limits[lm, t], - - post_flow[lm, lc, k] - emergency_limits[lm, t], + -post_flow[lm, lc, k] - emergency_limits[lm, t], ) end - + # Offer pre-contingency violations for lm in 1:L if pre_v[lm, k] > 1e-5 _offer( filters[t], Violation( - time=t, - monitored_line=instance.lines[lm], - outage_line=nothing, - amount=pre_v[lm, k], + time = t, + monitored_line = instance.lines[lm], + outage_line = nothing, + amount = pre_v[lm, k], ), ) end end - + # Offer post-contingency violations for lm in 1:L, lc in 1:L if post_v[lm, lc, k] > 1e-5 && is_vulnerable[lc] _offer( filters[t], Violation( - time=t, - monitored_line=instance.lines[lm], - outage_line=instance.lines[lc], - amount=post_v[lm, lc, k], + time = t, + monitored_line = instance.lines[lm], + outage_line = instance.lines[lc], + amount = post_v[lm, lc, k], ), ) end end end - + violations = Violation[] for t in 1:instance.time append!(violations, _query(filters[t])) end - + return violations end diff --git a/src/sensitivity.jl b/src/sensitivity.jl index 48c559f..691b3a7 100644 --- a/src/sensitivity.jl +++ b/src/sensitivity.jl @@ -13,15 +13,17 @@ M[l.offset, b.offset] indicates the amount of power (in MW) that flows through transmission line l when 1 MW of power is injected at the slack bus (the bus that has offset zero) and withdrawn from b. """ -function _injection_shift_factors(; buses::Array{Bus}, lines::Array{TransmissionLine}) +function _injection_shift_factors(; + buses::Array{Bus}, + lines::Array{TransmissionLine}, +) susceptance = _susceptance_matrix(lines) - incidence = _reduced_incidence_matrix(lines=lines, buses=buses) + incidence = _reduced_incidence_matrix(lines = lines, buses = buses) laplacian = transpose(incidence) * susceptance * incidence isf = susceptance * incidence * inv(Array(laplacian)) return isf end - """ _reduced_incidence_matrix(; buses::Array{Bus}, lines::Array{TransmissionLine}) @@ -31,7 +33,10 @@ is the number of buses and L is the number of lines. For each row, there is a 1 element and a -1 element, indicating the source and target buses, respectively, for that line. """ -function _reduced_incidence_matrix(; buses::Array{Bus}, lines::Array{TransmissionLine}) +function _reduced_incidence_matrix(; + buses::Array{Bus}, + lines::Array{TransmissionLine}, +) matrix = spzeros(Float64, length(lines), length(buses) - 1) for line in lines if line.source.offset > 0 @@ -41,7 +46,7 @@ function _reduced_incidence_matrix(; buses::Array{Bus}, lines::Array{Transmissio matrix[line.offset, line.target.offset] = -1 end end - matrix + return matrix end """ @@ -54,7 +59,6 @@ function _susceptance_matrix(lines::Array{TransmissionLine}) return Diagonal([l.susceptance for l in lines]) end - """ _line_outage_factors(; buses, lines, isf) @@ -63,19 +67,13 @@ Returns a LxL matrix containing the Line Outage Distribution Factors (LODFs) for the given network. This matrix how does the pre-contingency flow change when each individual transmission line is removed. """ -function _line_outage_factors( - ; - buses::Array{Bus, 1}, - lines::Array{TransmissionLine, 1}, +function _line_outage_factors(; + buses::Array{Bus,1}, + lines::Array{TransmissionLine,1}, isf::Array{Float64,2}, -) :: Array{Float64,2} +)::Array{Float64,2} n_lines, n_buses = size(isf) - incidence = Array( - _reduced_incidence_matrix( - lines=lines, - buses=buses, - ), - ) + incidence = Array(_reduced_incidence_matrix(lines = lines, buses = buses)) lodf::Array{Float64,2} = isf * transpose(incidence) m, n = size(lodf) for i in 1:n diff --git a/src/sysimage.jl b/src/sysimage.jl index 259b6fc..804c702 100644 --- a/src/sysimage.jl +++ b/src/sysimage.jl @@ -10,17 +10,11 @@ using JuMP using MathOptInterface using SparseArrays -pkg = [ - :DataStructures, - :JSON, - :JuMP, - :MathOptInterface, - :SparseArrays, -] +pkg = [:DataStructures, :JSON, :JuMP, :MathOptInterface, :SparseArrays] @info "Building system image..." create_sysimage( pkg, - precompile_statements_file="build/precompile.jl", - sysimage_path="build/sysimage.so", + precompile_statements_file = "build/precompile.jl", + sysimage_path = "build/sysimage.so", ) diff --git a/src/validate.jl b/src/validate.jl index 989fcbf..2ac6d3e 100644 --- a/src/validate.jl +++ b/src/validate.jl @@ -18,9 +18,9 @@ Returns the number of validation errors found. """ function repair!(instance::UnitCommitmentInstance)::Int n_errors = 0 - + for g in instance.units - + # Startup costs and delays must be increasing for s in 2:length(g.startup_categories) if g.startup_categories[s].delay <= g.startup_categories[s-1].delay @@ -31,7 +31,7 @@ function repair!(instance::UnitCommitmentInstance)::Int g.startup_categories[s].delay = new_value n_errors += 1 end - + if g.startup_categories[s].cost < g.startup_categories[s-1].cost prev_value = g.startup_categories[s].cost new_value = g.startup_categories[s-1].cost @@ -40,9 +40,8 @@ function repair!(instance::UnitCommitmentInstance)::Int g.startup_categories[s].cost = new_value n_errors += 1 end - end - + for t in 1:instance.time # Production cost curve should be convex for k in 2:length(g.cost_segments) @@ -67,19 +66,16 @@ function repair!(instance::UnitCommitmentInstance)::Int end end end - - + return n_errors end - function validate(instance_filename::String, solution_filename::String) instance = UnitCommitment.read(instance_filename) solution = JSON.parse(open(solution_filename)) return validate(instance, solution) end - """ validate(instance, solution)::Bool @@ -92,38 +88,39 @@ This function is implemented independently from the optimization model in producing valid solutions. It can also be used to verify the solutions produced by other optimization packages. """ -function validate(instance::UnitCommitmentInstance, - solution::Union{Dict,OrderedDict}; - )::Bool +function validate( + instance::UnitCommitmentInstance, + solution::Union{Dict,OrderedDict}, +)::Bool err_count = 0 err_count += _validate_units(instance, solution) err_count += _validate_reserve_and_demand(instance, solution) - + if err_count > 0 @error "Found $err_count validation errors" return false end - + return true end - -function _validate_units(instance, solution; tol=0.01) +function _validate_units(instance, solution; tol = 0.01) err_count = 0 - + for unit in instance.units production = solution["Production (MW)"][unit.name] reserve = solution["Reserve (MW)"][unit.name] actual_production_cost = solution["Production cost (\$)"][unit.name] actual_startup_cost = solution["Startup cost (\$)"][unit.name] is_on = bin(solution["Is on"][unit.name]) - + for t in 1:instance.time # Auxiliary variables if t == 1 is_starting_up = (unit.initial_status < 0) && is_on[t] is_shutting_down = (unit.initial_status > 0) && !is_on[t] - ramp_up = max(0, production[t] + reserve[t] - unit.initial_power) + ramp_up = + max(0, production[t] + reserve[t] - unit.initial_power) ramp_down = max(0, unit.initial_power - production[t]) else is_starting_up = !is_on[t-1] && is_on[t] @@ -131,7 +128,7 @@ function _validate_units(instance, solution; tol=0.01) ramp_up = max(0, production[t] + reserve[t] - production[t-1]) ramp_down = max(0, production[t-1] - production[t]) end - + # Compute production costs production_cost, startup_cost = 0, 0 if is_on[t] @@ -143,84 +140,133 @@ function _validate_units(instance, solution; tol=0.01) residual = max(0, residual - s.mw[t]) end end - + # Production should be non-negative if production[t] < -tol - @error @sprintf("Unit %s produces negative amount of power at time %d (%.2f)", - unit.name, t, production[t]) + @error @sprintf( + "Unit %s produces negative amount of power at time %d (%.2f)", + unit.name, + t, + production[t] + ) err_count += 1 end - + # Verify must-run if !is_on[t] && unit.must_run[t] - @error @sprintf("Must-run unit %s is offline at time %d", - unit.name, t) + @error @sprintf( + "Must-run unit %s is offline at time %d", + unit.name, + t + ) err_count += 1 end - + # Verify reserve eligibility if !unit.provides_spinning_reserves[t] && reserve[t] > tol - @error @sprintf("Unit %s is not eligible to provide spinning reserves at time %d", - unit.name, t) + @error @sprintf( + "Unit %s is not eligible to provide spinning reserves at time %d", + unit.name, + t + ) err_count += 1 end - + # If unit is on, must produce at least its minimum power if is_on[t] && (production[t] < unit.min_power[t] - tol) - @error @sprintf("Unit %s produces below its minimum limit at time %d (%.2f < %.2f)", - unit.name, t, production[t], unit.min_power[t]) + @error @sprintf( + "Unit %s produces below its minimum limit at time %d (%.2f < %.2f)", + unit.name, + t, + production[t], + unit.min_power[t] + ) err_count += 1 end - + # If unit is on, must produce at most its maximum power - if is_on[t] && (production[t] + reserve[t] > unit.max_power[t] + tol) - @error @sprintf("Unit %s produces above its maximum limit at time %d (%.2f + %.2f> %.2f)", - unit.name, t, production[t], reserve[t], unit.max_power[t]) + if is_on[t] && + (production[t] + reserve[t] > unit.max_power[t] + tol) + @error @sprintf( + "Unit %s produces above its maximum limit at time %d (%.2f + %.2f> %.2f)", + unit.name, + t, + production[t], + reserve[t], + unit.max_power[t] + ) err_count += 1 end - + # If unit is off, must produce zero if !is_on[t] && production[t] + reserve[t] > tol - @error @sprintf("Unit %s produces power at time %d while off", - unit.name, t) + @error @sprintf( + "Unit %s produces power at time %d while off", + unit.name, + t + ) err_count += 1 end - + # Startup limit if is_starting_up && (ramp_up > unit.startup_limit + tol) - @error @sprintf("Unit %s exceeds startup limit at time %d (%.2f > %.2f)", - unit.name, t, ramp_up, unit.startup_limit) + @error @sprintf( + "Unit %s exceeds startup limit at time %d (%.2f > %.2f)", + unit.name, + t, + ramp_up, + unit.startup_limit + ) err_count += 1 end # Shutdown limit if is_shutting_down && (ramp_down > unit.shutdown_limit + tol) - @error @sprintf("Unit %s exceeds shutdown limit at time %d (%.2f > %.2f)", - unit.name, t, ramp_down, unit.shutdown_limit) + @error @sprintf( + "Unit %s exceeds shutdown limit at time %d (%.2f > %.2f)", + unit.name, + t, + ramp_down, + unit.shutdown_limit + ) err_count += 1 end # Ramp-up limit - if !is_starting_up && !is_shutting_down && (ramp_up > unit.ramp_up_limit + tol) - @error @sprintf("Unit %s exceeds ramp up limit at time %d (%.2f > %.2f)", - unit.name, t, ramp_up, unit.ramp_up_limit) + if !is_starting_up && + !is_shutting_down && + (ramp_up > unit.ramp_up_limit + tol) + @error @sprintf( + "Unit %s exceeds ramp up limit at time %d (%.2f > %.2f)", + unit.name, + t, + ramp_up, + unit.ramp_up_limit + ) err_count += 1 end # Ramp-down limit - if !is_starting_up && !is_shutting_down && (ramp_down > unit.ramp_down_limit + tol) - @error @sprintf("Unit %s exceeds ramp down limit at time %d (%.2f > %.2f)", - unit.name, t, ramp_down, unit.ramp_down_limit) + if !is_starting_up && + !is_shutting_down && + (ramp_down > unit.ramp_down_limit + tol) + @error @sprintf( + "Unit %s exceeds ramp down limit at time %d (%.2f > %.2f)", + unit.name, + t, + ramp_down, + unit.ramp_down_limit + ) err_count += 1 end - + # Verify startup costs & minimum downtime if is_starting_up - + # Calculate how much time the unit has been offline time_down = 0 for k in 1:(t-1) - if !is_on[t - k] + if !is_on[t-k] time_down += 1 else break @@ -233,29 +279,32 @@ function _validate_units(instance, solution; tol=0.01) end time_down += initial_down end - + # Calculate startup costs for c in unit.startup_categories if time_down >= c.delay startup_cost = c.cost end end - + # Check minimum downtime if time_down < unit.min_downtime - @error @sprintf("Unit %s violates minimum downtime at time %d", - unit.name, t) + @error @sprintf( + "Unit %s violates minimum downtime at time %d", + unit.name, + t + ) err_count += 1 end end - + # Verify minimum uptime if is_shutting_down - + # Calculate how much time the unit has been online time_up = 0 for k in 1:(t-1) - if is_on[t - k] + if is_on[t-k] time_up += 1 else break @@ -268,61 +317,70 @@ function _validate_units(instance, solution; tol=0.01) end time_up += initial_up end - + if (t == time_up + 1) && (unit.initial_status > 0) time_up += unit.initial_status end - + # Check minimum uptime if time_up < unit.min_uptime - @error @sprintf("Unit %s violates minimum uptime at time %d", - unit.name, t) + @error @sprintf( + "Unit %s violates minimum uptime at time %d", + unit.name, + t + ) err_count += 1 end end - + # Verify production costs if abs(actual_production_cost[t] - production_cost) > 1.00 - @error @sprintf("Unit %s has unexpected production cost at time %d (%.2f should be %.2f)", - unit.name, t, actual_production_cost[t], production_cost) + @error @sprintf( + "Unit %s has unexpected production cost at time %d (%.2f should be %.2f)", + unit.name, + t, + actual_production_cost[t], + production_cost + ) err_count += 1 end # Verify startup costs if abs(actual_startup_cost[t] - startup_cost) > 1.00 - @error @sprintf("Unit %s has unexpected startup cost at time %d (%.2f should be %.2f)", - unit.name, t, actual_startup_cost[t], startup_cost) + @error @sprintf( + "Unit %s has unexpected startup cost at time %d (%.2f should be %.2f)", + unit.name, + t, + actual_startup_cost[t], + startup_cost + ) err_count += 1 end - end end - + return err_count end - -function _validate_reserve_and_demand(instance, solution, tol=0.01) +function _validate_reserve_and_demand(instance, solution, tol = 0.01) err_count = 0 for t in 1:instance.time load_curtail = 0 fixed_load = sum(b.load[t] for b in instance.buses) ps_load = sum( - solution["Price-sensitive loads (MW)"][ps.name][t] - for ps in instance.price_sensitive_loads - ) - production = sum( - solution["Production (MW)"][g.name][t] - for g in instance.units + solution["Price-sensitive loads (MW)"][ps.name][t] for + ps in instance.price_sensitive_loads ) + production = + sum(solution["Production (MW)"][g.name][t] for g in instance.units) if "Load curtail (MW)" in keys(solution) load_curtail = sum( - solution["Load curtail (MW)"][b.name][t] - for b in instance.buses + solution["Load curtail (MW)"][b.name][t] for + b in instance.buses ) end balance = fixed_load - load_curtail - production + ps_load - + # Verify that production equals demand if abs(balance) > tol @error @sprintf( @@ -335,9 +393,10 @@ function _validate_reserve_and_demand(instance, solution, tol=0.01) ) err_count += 1 end - + # Verify spinning reserves - reserve = sum(solution["Reserve (MW)"][g.name][t] for g in instance.units) + reserve = + sum(solution["Reserve (MW)"][g.name][t] for g in instance.units) if reserve < instance.reserves.spinning[t] - tol @error @sprintf( "Insufficient spinning reserves at time %d (%.2f should be %.2f)", @@ -348,7 +407,6 @@ function _validate_reserve_and_demand(instance, solution, tol=0.01) err_count += 1 end end - + return err_count end - \ No newline at end of file diff --git a/test/convert_test.jl b/test/convert_test.jl index f6ee042..6c8fb0f 100644 --- a/test/convert_test.jl +++ b/test/convert_test.jl @@ -6,14 +6,17 @@ using UnitCommitment @testset "convert" begin @testset "EGRET solution" begin - solution = UnitCommitment._read_egret_solution("fixtures/egret_output.json.gz") + solution = + UnitCommitment._read_egret_solution("fixtures/egret_output.json.gz") for attr in ["Is on", "Production (MW)", "Production cost (\$)"] @test attr in keys(solution) @test "115_STEAM_1" in keys(solution[attr]) @test length(solution[attr]["115_STEAM_1"]) == 48 end - @test solution["Production cost (\$)"]["315_CT_6"][15:20] == [0., 0., 884.44, 1470.71, 1470.71, 884.44] - @test solution["Startup cost (\$)"]["315_CT_6"][15:20] == [0., 0., 5665.23, 0., 0., 0.] + @test solution["Production cost (\$)"]["315_CT_6"][15:20] == + [0.0, 0.0, 884.44, 1470.71, 1470.71, 884.44] + @test solution["Startup cost (\$)"]["315_CT_6"][15:20] == + [0.0, 0.0, 5665.23, 0.0, 0.0, 0.0] @test length(keys(solution["Is on"])) == 154 end end diff --git a/test/initcond_test.jl b/test/initcond_test.jl index ba773f1..4cb127c 100644 --- a/test/initcond_test.jl +++ b/test/initcond_test.jl @@ -8,21 +8,21 @@ using UnitCommitment, Cbc, JuMP # Load instance instance = UnitCommitment.read("$(pwd())/fixtures/case118-initcond.json.gz") optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0) - + # All units should have unknown initial conditions for g in instance.units @test g.initial_power === nothing @test g.initial_status === nothing end - + # Generate initial conditions UnitCommitment.generate_initial_conditions!(instance, optimizer) - + # All units should now have known initial conditions for g in instance.units @test g.initial_power !== nothing @test g.initial_status !== nothing end - + # TODO: Check that initial conditions are feasible end diff --git a/test/instance_test.jl b/test/instance_test.jl index d4658ee..4fe9b9f 100644 --- a/test/instance_test.jl +++ b/test/instance_test.jl @@ -15,46 +15,46 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip @test length(instance.price_sensitive_loads) == 1 @test instance.time == 4 - @test instance.lines[5].name == "l5" - @test instance.lines[5].source.name == "b2" - @test instance.lines[5].target.name == "b5" - @test instance.lines[5].reactance ≈ 0.17388 - @test instance.lines[5].susceptance ≈ 10.037550333 - @test instance.lines[5].normal_flow_limit == [1e8 for t in 1:4] + @test instance.lines[5].name == "l5" + @test instance.lines[5].source.name == "b2" + @test instance.lines[5].target.name == "b5" + @test instance.lines[5].reactance ≈ 0.17388 + @test instance.lines[5].susceptance ≈ 10.037550333 + @test instance.lines[5].normal_flow_limit == [1e8 for t in 1:4] @test instance.lines[5].emergency_flow_limit == [1e8 for t in 1:4] - @test instance.lines[5].flow_limit_penalty == [5e3 for t in 1:4] - - @test instance.lines[1].name == "l1" - @test instance.lines[1].source.name == "b1" - @test instance.lines[1].target.name == "b2" - @test instance.lines[1].reactance ≈ 0.059170 - @test instance.lines[1].susceptance ≈ 29.496860773945 - @test instance.lines[1].normal_flow_limit == [300.0 for t in 1:4] + @test instance.lines[5].flow_limit_penalty == [5e3 for t in 1:4] + + @test instance.lines[1].name == "l1" + @test instance.lines[1].source.name == "b1" + @test instance.lines[1].target.name == "b2" + @test instance.lines[1].reactance ≈ 0.059170 + @test instance.lines[1].susceptance ≈ 29.496860773945 + @test instance.lines[1].normal_flow_limit == [300.0 for t in 1:4] @test instance.lines[1].emergency_flow_limit == [400.0 for t in 1:4] - @test instance.lines[1].flow_limit_penalty == [1e3 for t in 1:4] + @test instance.lines[1].flow_limit_penalty == [1e3 for t in 1:4] @test instance.buses[9].name == "b9" @test instance.buses[9].load == [35.36638, 33.25495, 31.67138, 31.14353] unit = instance.units[1] - @test unit.name == "g1" - @test unit.bus.name == "b1" - @test unit.ramp_up_limit == 1e6 - @test unit.ramp_down_limit == 1e6 - @test unit.startup_limit == 1e6 - @test unit.shutdown_limit == 1e6 - @test unit.must_run == [false for t in 1:4] - @test unit.min_power_cost == [1400. for t in 1:4] - @test unit.min_uptime == 1 - @test unit.min_downtime == 1 - @test unit.provides_spinning_reserves == [true for t in 1:4] + @test unit.name == "g1" + @test unit.bus.name == "b1" + @test unit.ramp_up_limit == 1e6 + @test unit.ramp_down_limit == 1e6 + @test unit.startup_limit == 1e6 + @test unit.shutdown_limit == 1e6 + @test unit.must_run == [false for t in 1:4] + @test unit.min_power_cost == [1400.0 for t in 1:4] + @test unit.min_uptime == 1 + @test unit.min_downtime == 1 + @test unit.provides_spinning_reserves == [true for t in 1:4] for t in 1:1 - @test unit.cost_segments[1].mw[t] == 10.0 - @test unit.cost_segments[2].mw[t] == 20.0 - @test unit.cost_segments[3].mw[t] == 5.0 - @test unit.cost_segments[1].cost[t] ≈ 20.0 - @test unit.cost_segments[2].cost[t] ≈ 30.0 - @test unit.cost_segments[3].cost[t] ≈ 40.0 + @test unit.cost_segments[1].mw[t] == 10.0 + @test unit.cost_segments[2].mw[t] == 20.0 + @test unit.cost_segments[3].mw[t] == 5.0 + @test unit.cost_segments[1].cost[t] ≈ 20.0 + @test unit.cost_segments[2].cost[t] ≈ 30.0 + @test unit.cost_segments[3].cost[t] ≈ 40.0 end @test length(unit.startup_categories) == 3 @test unit.startup_categories[1].delay == 1 @@ -63,42 +63,42 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip @test unit.startup_categories[1].cost == 1000.0 @test unit.startup_categories[2].cost == 1500.0 @test unit.startup_categories[3].cost == 2000.0 - + unit = instance.units[2] - @test unit.name == "g2" + @test unit.name == "g2" @test unit.must_run == [false for t in 1:4] unit = instance.units[3] - @test unit.name == "g3" - @test unit.bus.name == "b3" - @test unit.ramp_up_limit == 70.0 - @test unit.ramp_down_limit == 70.0 - @test unit.startup_limit == 70.0 - @test unit.shutdown_limit == 70.0 - @test unit.must_run == [true for t in 1:4] - @test unit.min_power_cost == [0. for t in 1:4] - @test unit.min_uptime == 1 - @test unit.min_downtime == 1 - @test unit.provides_spinning_reserves == [true for t in 1:4] + @test unit.name == "g3" + @test unit.bus.name == "b3" + @test unit.ramp_up_limit == 70.0 + @test unit.ramp_down_limit == 70.0 + @test unit.startup_limit == 70.0 + @test unit.shutdown_limit == 70.0 + @test unit.must_run == [true for t in 1:4] + @test unit.min_power_cost == [0.0 for t in 1:4] + @test unit.min_uptime == 1 + @test unit.min_downtime == 1 + @test unit.provides_spinning_reserves == [true for t in 1:4] for t in 1:4 - @test unit.cost_segments[1].mw[t] ≈ 33 - @test unit.cost_segments[2].mw[t] ≈ 33 - @test unit.cost_segments[3].mw[t] ≈ 34 + @test unit.cost_segments[1].mw[t] ≈ 33 + @test unit.cost_segments[2].mw[t] ≈ 33 + @test unit.cost_segments[3].mw[t] ≈ 34 @test unit.cost_segments[1].cost[t] ≈ 33.75 @test unit.cost_segments[2].cost[t] ≈ 38.04 @test unit.cost_segments[3].cost[t] ≈ 44.77853 end - + @test instance.reserves.spinning == zeros(4) - + @test instance.contingencies[1].lines == [instance.lines[1]] @test instance.contingencies[1].units == [] - + load = instance.price_sensitive_loads[1] - @test load.name == "ps1" + @test load.name == "ps1" @test load.bus.name == "b3" - @test load.revenue == [100. for t in 1:4] - @test load.demand == [50. for t in 1:4] + @test load.revenue == [100.0 for t in 1:4] + @test load.demand == [50.0 for t in 1:4] end @testset "read sub-hourly" begin @@ -114,11 +114,11 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip @test unit.startup_categories[3].delay == 6 @test unit.initial_status == -200 end - + @testset "slice" begin instance = UnitCommitment.read_benchmark("test/case14") modified = UnitCommitment.slice(instance, 1:2) - + # Should update all time-dependent fields @test modified.time == 2 @test length(modified.power_balance_penalty) == 2 @@ -146,11 +146,13 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON, GZip @test length(ps.demand) == 2 @test length(ps.revenue) == 2 end - + # Should be able to build model without errors optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0) - model = build_model(instance=modified, - optimizer=optimizer, - variable_names=true) + model = build_model( + instance = modified, + optimizer = optimizer, + variable_names = true, + ) end end diff --git a/test/model_test.jl b/test/model_test.jl index c736489..dae7515 100644 --- a/test/model_test.jl +++ b/test/model_test.jl @@ -12,9 +12,9 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP end optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0) model = build_model( - instance=instance, - optimizer=optimizer, - variable_names=true, + instance = instance, + optimizer = optimizer, + variable_names = true, ) @test name(model[:is_on]["g1", 1]) == "is_on[g1,1]" @@ -27,7 +27,7 @@ using UnitCommitment, LinearAlgebra, Cbc, JuMP UnitCommitment.write(filename, solution) loaded = JSON.parsefile(filename) @test length(loaded["Is on"]) == 6 - + # Verify solution @test UnitCommitment.validate(instance, solution) diff --git a/test/screening_test.jl b/test/screening_test.jl index 372f780..7d2201f 100644 --- a/test/screening_test.jl +++ b/test/screening_test.jl @@ -8,81 +8,81 @@ import UnitCommitment: Violation, _offer, _query @testset "Screening" begin @testset "Violation filter" begin instance = UnitCommitment.read_benchmark("test/case14") - filter = UnitCommitment.ViolationFilter(max_per_line=1, max_total=2) + filter = UnitCommitment.ViolationFilter(max_per_line = 1, max_total = 2) _offer( filter, Violation( - time=1, - monitored_line=instance.lines[1], - outage_line=nothing, - amount=100., + time = 1, + monitored_line = instance.lines[1], + outage_line = nothing, + amount = 100.0, ), ) _offer( filter, Violation( - time=1, - monitored_line=instance.lines[1], - outage_line=instance.lines[1], - amount=300., + time = 1, + monitored_line = instance.lines[1], + outage_line = instance.lines[1], + amount = 300.0, ), ) _offer( filter, Violation( - time=1, - monitored_line=instance.lines[1], - outage_line=instance.lines[5], - amount=500., + time = 1, + monitored_line = instance.lines[1], + outage_line = instance.lines[5], + amount = 500.0, ), ) _offer( filter, Violation( - time=1, - monitored_line=instance.lines[1], - outage_line=instance.lines[4], - amount=400., + time = 1, + monitored_line = instance.lines[1], + outage_line = instance.lines[4], + amount = 400.0, ), ) _offer( filter, Violation( - time=1, - monitored_line=instance.lines[2], - outage_line=instance.lines[1], - amount=200., + time = 1, + monitored_line = instance.lines[2], + outage_line = instance.lines[1], + amount = 200.0, ), ) _offer( filter, Violation( - time=1, - monitored_line=instance.lines[2], - outage_line=instance.lines[8], - amount=100., - ) + time = 1, + monitored_line = instance.lines[2], + outage_line = instance.lines[8], + amount = 100.0, + ), ) actual = _query(filter) expected = [ Violation( - time=1, - monitored_line=instance.lines[2], - outage_line=instance.lines[1], - amount=200., + time = 1, + monitored_line = instance.lines[2], + outage_line = instance.lines[1], + amount = 200.0, ), Violation( - time=1, - monitored_line=instance.lines[1], - outage_line=instance.lines[5], - amount=500., + time = 1, + monitored_line = instance.lines[1], + outage_line = instance.lines[5], + amount = 500.0, ), ] @test actual == expected end - + @testset "find_violations" begin instance = UnitCommitment.read_benchmark("test/case14") for line in instance.lines, t in 1:instance.time @@ -90,22 +90,22 @@ import UnitCommitment: Violation, _offer, _query line.emergency_flow_limit[t] = 1.0 end isf = UnitCommitment._injection_shift_factors( - lines=instance.lines, - buses=instance.buses, + lines = instance.lines, + buses = instance.buses, ) lodf = UnitCommitment._line_outage_factors( - lines=instance.lines, - buses=instance.buses, - isf=isf, + lines = instance.lines, + buses = instance.buses, + isf = isf, ) inj = [1000.0 for b in 1:13, t in 1:instance.time] overflow = [0.0 for l in instance.lines, t in 1:instance.time] violations = UnitCommitment._find_violations( - instance=instance, - net_injections=inj, - overflow=overflow, - isf=isf, - lodf=lodf, + instance = instance, + net_injections = inj, + overflow = overflow, + isf = isf, + lodf = lodf, ) @test length(violations) == 20 end diff --git a/test/sensitivity_test.jl b/test/sensitivity_test.jl index 0676a5a..3032a93 100644 --- a/test/sensitivity_test.jl +++ b/test/sensitivity_test.jl @@ -9,117 +9,137 @@ using UnitCommitment, Test, LinearAlgebra instance = UnitCommitment.read_benchmark("test/case14") actual = UnitCommitment._susceptance_matrix(instance.lines) @test size(actual) == (20, 20) - expected = Diagonal([29.5, 7.83, 8.82, 9.9, 10.04, - 10.2, 41.45, 8.35, 3.14, 6.93, - 8.77, 6.82, 13.4, 9.91, 15.87, - 20.65, 6.46, 9.09, 8.73, 5.02]) - @test round.(actual, digits=2) == expected + expected = Diagonal([ + 29.5, + 7.83, + 8.82, + 9.9, + 10.04, + 10.2, + 41.45, + 8.35, + 3.14, + 6.93, + 8.77, + 6.82, + 13.4, + 9.91, + 15.87, + 20.65, + 6.46, + 9.09, + 8.73, + 5.02, + ]) + @test round.(actual, digits = 2) == expected end - + @testset "Reduced incidence matrix" begin instance = UnitCommitment.read_benchmark("test/case14") actual = UnitCommitment._reduced_incidence_matrix( - lines=instance.lines, - buses=instance.buses, + lines = instance.lines, + buses = instance.buses, ) @test size(actual) == (20, 13) - @test actual[1, 1] == -1.0 - @test actual[3, 1] == 1.0 - @test actual[4, 1] == 1.0 - @test actual[5, 1] == 1.0 - @test actual[3, 2] == -1.0 - @test actual[6, 2] == 1.0 - @test actual[4, 3] == -1.0 - @test actual[6, 3] == -1.0 - @test actual[7, 3] == 1.0 - @test actual[8, 3] == 1.0 - @test actual[9, 3] == 1.0 - @test actual[2, 4] == -1.0 - @test actual[5, 4] == -1.0 - @test actual[7, 4] == -1.0 - @test actual[10, 4] == 1.0 - @test actual[10, 5] == -1.0 - @test actual[11, 5] == 1.0 - @test actual[12, 5] == 1.0 - @test actual[13, 5] == 1.0 - @test actual[8, 6] == -1.0 - @test actual[14, 6] == 1.0 - @test actual[15, 6] == 1.0 - @test actual[14, 7] == -1.0 - @test actual[9, 8] == -1.0 - @test actual[15, 8] == -1.0 - @test actual[16, 8] == 1.0 - @test actual[17, 8] == 1.0 - @test actual[16, 9] == -1.0 - @test actual[18, 9] == 1.0 - @test actual[11, 10] == -1.0 - @test actual[18, 10] == -1.0 - @test actual[12, 11] == -1.0 - @test actual[19, 11] == 1.0 - @test actual[13, 12] == -1.0 - @test actual[19, 12] == -1.0 - @test actual[20, 12] == 1.0 - @test actual[17, 13] == -1.0 - @test actual[20, 13] == -1.0 + @test actual[1, 1] == -1.0 + @test actual[3, 1] == 1.0 + @test actual[4, 1] == 1.0 + @test actual[5, 1] == 1.0 + @test actual[3, 2] == -1.0 + @test actual[6, 2] == 1.0 + @test actual[4, 3] == -1.0 + @test actual[6, 3] == -1.0 + @test actual[7, 3] == 1.0 + @test actual[8, 3] == 1.0 + @test actual[9, 3] == 1.0 + @test actual[2, 4] == -1.0 + @test actual[5, 4] == -1.0 + @test actual[7, 4] == -1.0 + @test actual[10, 4] == 1.0 + @test actual[10, 5] == -1.0 + @test actual[11, 5] == 1.0 + @test actual[12, 5] == 1.0 + @test actual[13, 5] == 1.0 + @test actual[8, 6] == -1.0 + @test actual[14, 6] == 1.0 + @test actual[15, 6] == 1.0 + @test actual[14, 7] == -1.0 + @test actual[9, 8] == -1.0 + @test actual[15, 8] == -1.0 + @test actual[16, 8] == 1.0 + @test actual[17, 8] == 1.0 + @test actual[16, 9] == -1.0 + @test actual[18, 9] == 1.0 + @test actual[11, 10] == -1.0 + @test actual[18, 10] == -1.0 + @test actual[12, 11] == -1.0 + @test actual[19, 11] == 1.0 + @test actual[13, 12] == -1.0 + @test actual[19, 12] == -1.0 + @test actual[20, 12] == 1.0 + @test actual[17, 13] == -1.0 + @test actual[20, 13] == -1.0 end - + @testset "Injection Shift Factors (ISF)" begin instance = UnitCommitment.read_benchmark("test/case14") actual = UnitCommitment._injection_shift_factors( - lines=instance.lines, - buses=instance.buses, + lines = instance.lines, + buses = instance.buses, ) @test size(actual) == (20, 13) - @test round.(actual, digits=2) == [ - -0.84 -0.75 -0.67 -0.61 -0.63 -0.66 -0.66 -0.65 -0.65 -0.64 -0.63 -0.63 -0.64; - -0.16 -0.25 -0.33 -0.39 -0.37 -0.34 -0.34 -0.35 -0.35 -0.36 -0.37 -0.37 -0.36; - 0.03 -0.53 -0.15 -0.1 -0.12 -0.14 -0.14 -0.14 -0.13 -0.13 -0.12 -0.12 -0.13; - 0.06 -0.14 -0.32 -0.22 -0.25 -0.3 -0.3 -0.29 -0.28 -0.27 -0.25 -0.26 -0.27; - 0.08 -0.07 -0.2 -0.29 -0.26 -0.22 -0.22 -0.22 -0.23 -0.25 -0.26 -0.26 -0.24; - 0.03 0.47 -0.15 -0.1 -0.12 -0.14 -0.14 -0.14 -0.13 -0.13 -0.12 -0.12 -0.13; - 0.08 0.31 0.5 -0.3 -0.03 0.36 0.36 0.28 0.23 0.1 -0.0 0.02 0.17; - 0.0 0.01 0.02 -0.01 -0.22 -0.63 -0.63 -0.45 -0.41 -0.32 -0.24 -0.25 -0.36; - 0.0 0.01 0.01 -0.01 -0.12 -0.17 -0.17 -0.26 -0.24 -0.18 -0.14 -0.14 -0.21; - -0.0 -0.02 -0.03 0.02 -0.66 -0.2 -0.2 -0.29 -0.36 -0.5 -0.63 -0.61 -0.43; - -0.0 -0.01 -0.02 0.01 0.21 -0.12 -0.12 -0.17 -0.28 -0.53 0.18 0.15 -0.03; - -0.0 -0.0 -0.0 0.0 0.03 -0.02 -0.02 -0.03 -0.02 0.01 -0.52 -0.17 -0.09; - -0.0 -0.01 -0.01 0.01 0.11 -0.06 -0.06 -0.09 -0.05 0.02 -0.28 -0.59 -0.31; - -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 ; - 0.0 0.01 0.02 -0.01 -0.22 0.37 0.37 -0.45 -0.41 -0.32 -0.24 -0.25 -0.36; - 0.0 0.01 0.02 -0.01 -0.21 0.12 0.12 0.17 -0.72 -0.47 -0.18 -0.15 0.03; - 0.0 0.01 0.01 -0.01 -0.14 0.08 0.08 0.12 0.07 -0.03 -0.2 -0.24 -0.6 ; - 0.0 0.01 0.02 -0.01 -0.21 0.12 0.12 0.17 0.28 -0.47 -0.18 -0.15 0.03; - -0.0 -0.0 -0.0 0.0 0.03 -0.02 -0.02 -0.03 -0.02 0.01 0.48 -0.17 -0.09; - -0.0 -0.01 -0.01 0.01 0.14 -0.08 -0.08 -0.12 -0.07 0.03 0.2 0.24 -0.4 ] + @test round.(actual, digits = 2) == [ + -0.84 -0.75 -0.67 -0.61 -0.63 -0.66 -0.66 -0.65 -0.65 -0.64 -0.63 -0.63 -0.64 + -0.16 -0.25 -0.33 -0.39 -0.37 -0.34 -0.34 -0.35 -0.35 -0.36 -0.37 -0.37 -0.36 + 0.03 -0.53 -0.15 -0.1 -0.12 -0.14 -0.14 -0.14 -0.13 -0.13 -0.12 -0.12 -0.13 + 0.06 -0.14 -0.32 -0.22 -0.25 -0.3 -0.3 -0.29 -0.28 -0.27 -0.25 -0.26 -0.27 + 0.08 -0.07 -0.2 -0.29 -0.26 -0.22 -0.22 -0.22 -0.23 -0.25 -0.26 -0.26 -0.24 + 0.03 0.47 -0.15 -0.1 -0.12 -0.14 -0.14 -0.14 -0.13 -0.13 -0.12 -0.12 -0.13 + 0.08 0.31 0.5 -0.3 -0.03 0.36 0.36 0.28 0.23 0.1 -0.0 0.02 0.17 + 0.0 0.01 0.02 -0.01 -0.22 -0.63 -0.63 -0.45 -0.41 -0.32 -0.24 -0.25 -0.36 + 0.0 0.01 0.01 -0.01 -0.12 -0.17 -0.17 -0.26 -0.24 -0.18 -0.14 -0.14 -0.21 + -0.0 -0.02 -0.03 0.02 -0.66 -0.2 -0.2 -0.29 -0.36 -0.5 -0.63 -0.61 -0.43 + -0.0 -0.01 -0.02 0.01 0.21 -0.12 -0.12 -0.17 -0.28 -0.53 0.18 0.15 -0.03 + -0.0 -0.0 -0.0 0.0 0.03 -0.02 -0.02 -0.03 -0.02 0.01 -0.52 -0.17 -0.09 + -0.0 -0.01 -0.01 0.01 0.11 -0.06 -0.06 -0.09 -0.05 0.02 -0.28 -0.59 -0.31 + -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -1.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 + 0.0 0.01 0.02 -0.01 -0.22 0.37 0.37 -0.45 -0.41 -0.32 -0.24 -0.25 -0.36 + 0.0 0.01 0.02 -0.01 -0.21 0.12 0.12 0.17 -0.72 -0.47 -0.18 -0.15 0.03 + 0.0 0.01 0.01 -0.01 -0.14 0.08 0.08 0.12 0.07 -0.03 -0.2 -0.24 -0.6 + 0.0 0.01 0.02 -0.01 -0.21 0.12 0.12 0.17 0.28 -0.47 -0.18 -0.15 0.03 + -0.0 -0.0 -0.0 0.0 0.03 -0.02 -0.02 -0.03 -0.02 0.01 0.48 -0.17 -0.09 + -0.0 -0.01 -0.01 0.01 0.14 -0.08 -0.08 -0.12 -0.07 0.03 0.2 0.24 -0.4 + ] end @testset "Line Outage Distribution Factors (LODF)" begin instance = UnitCommitment.read_benchmark("test/case14") isf_before = UnitCommitment._injection_shift_factors( - lines=instance.lines, - buses=instance.buses, + lines = instance.lines, + buses = instance.buses, ) lodf = UnitCommitment._line_outage_factors( - lines=instance.lines, - buses=instance.buses, - isf=isf_before, + lines = instance.lines, + buses = instance.buses, + isf = isf_before, ) for contingency in instance.contingencies for lc in contingency.lines prev_susceptance = lc.susceptance lc.susceptance = 0.0 - isf_after = UnitCommitment._injection_shift_factors( - lines=instance.lines, - buses=instance.buses, + isf_after = UnitCommitment._injection_shift_factors( + lines = instance.lines, + buses = instance.buses, ) lc.susceptance = prev_susceptance for lm in instance.lines expected = isf_after[lm.offset, :] - actual = isf_before[lm.offset, :] + - lodf[lm.offset, lc.offset] * isf_before[lc.offset, :] + actual = + isf_before[lm.offset, :] + + lodf[lm.offset, lc.offset] * isf_before[lc.offset, :] @test norm(expected - actual) < 1e-6 end end end end -end \ No newline at end of file +end diff --git a/test/validate_test.jl b/test/validate_test.jl index 9ca49c8..c1f047e 100644 --- a/test/validate_test.jl +++ b/test/validate_test.jl @@ -4,36 +4,40 @@ using UnitCommitment, JSON, GZip, DataStructures -parse_case14() = JSON.parse(GZip.gzopen("../instances/test/case14.json.gz"), - dicttype=()->DefaultOrderedDict(nothing)) +function parse_case14() + return JSON.parse( + GZip.gzopen("../instances/test/case14.json.gz"), + dicttype = () -> DefaultOrderedDict(nothing), + ) +end @testset "Validation" begin @testset "repair!" begin - @testset "Cost curve should be convex" begin json = parse_case14() - json["Generators"]["g1"]["Production cost curve (MW)"] = [100, 150, 200] - json["Generators"]["g1"]["Production cost curve (\$)"] = [10, 25, 30] - instance = UnitCommitment._from_json(json, repair=false) + json["Generators"]["g1"]["Production cost curve (MW)"] = + [100, 150, 200] + json["Generators"]["g1"]["Production cost curve (\$)"] = + [10, 25, 30] + instance = UnitCommitment._from_json(json, repair = false) @test UnitCommitment.repair!(instance) == 4 end - + @testset "Startup limit must be greater than Pmin" begin json = parse_case14() json["Generators"]["g1"]["Production cost curve (MW)"] = [100, 150] json["Generators"]["g1"]["Production cost curve (\$)"] = [100, 150] json["Generators"]["g1"]["Startup limit (MW)"] = 80 - instance = UnitCommitment._from_json(json, repair=false) + instance = UnitCommitment._from_json(json, repair = false) @test UnitCommitment.repair!(instance) == 1 end - + @testset "Startup costs and delays must be increasing" begin json = parse_case14() json["Generators"]["g1"]["Startup costs (\$)"] = [300, 200, 100] json["Generators"]["g1"]["Startup delays (h)"] = [8, 4, 2] - instance = UnitCommitment._from_json(json, repair=false) + instance = UnitCommitment._from_json(json, repair = false) @test UnitCommitment.repair!(instance) == 4 end - end end