diff --git a/PRASCapacityCredits.jl/src/EFC.jl b/PRASCapacityCredits.jl/src/EFC.jl index f445e4cd..a6862a49 100644 --- a/PRASCapacityCredits.jl/src/EFC.jl +++ b/PRASCapacityCredits.jl/src/EFC.jl @@ -107,83 +107,4 @@ function assess(sys_baseline::S, sys_augmented::S, end -function add_firmcapacity( - s1::SystemModel{N,L,T,P,E}, s2::SystemModel{N,L,T,P,E}, - region_shares::Vector{Tuple{String,Float64}} -) where {N,L,T,P,E} - n_regions = length(s1.regions.names) - n_region_allocs = length(region_shares) - - region_allocations = allocate_regions(s1.regions.names, region_shares) - efc_gens = similar(region_allocations) - - new_gen(i::Int) = Generators{N,L,T,P}( - ["_EFC_$i"], ["_EFC Calculation Dummy Generator"], - zeros(Int, 1, N), zeros(1, N), ones(1, N)) - - variable_gens = Generators{N,L,T,P}[] - variable_region_gen_idxs = similar(s1.region_gen_idxs) - - target_gens = similar(variable_gens) - target_region_gen_idxs = similar(s2.region_gen_idxs) - - ra_idx = 0 - - for r in 1:n_regions - - s1_range = s1.region_gen_idxs[r] - s2_range = s2.region_gen_idxs[r] - - if (ra_idx < n_region_allocs) && (r == first(region_allocations[ra_idx+1])) - - ra_idx += 1 - - variable_region_gen_idxs[r] = incr_range(s1_range, ra_idx-1, ra_idx) - target_region_gen_idxs[r] = incr_range(s2_range, ra_idx-1, ra_idx) - - gen = new_gen(ra_idx) - push!(variable_gens, gen) - push!(target_gens, gen) - efc_gens[ra_idx] = ( - first(s1_range) + ra_idx - 1, - last(region_allocations[ra_idx])) - - else - - variable_region_gen_idxs[r] = incr_range(s1_range, ra_idx) - target_region_gen_idxs[r] = incr_range(s2_range, ra_idx) - - end - - push!(variable_gens, s1.generators[s1_range]) - push!(target_gens, s2.generators[s2_range]) - - end - - sys_variable = SystemModel( - s1.regions, s1.interfaces, - vcat(variable_gens...), variable_region_gen_idxs, - s1.storages, s1.region_stor_idxs, - s1.generatorstorages, s1.region_genstor_idxs, - s1.lines, s1.interface_line_idxs, s1.timestamps) - - sys_target = SystemModel( - s2.regions, s2.interfaces, - vcat(target_gens...), target_region_gen_idxs, - s2.storages, s2.region_stor_idxs, - s2.generatorstorages, s2.region_genstor_idxs, - s2.lines, s2.interface_line_idxs, s2.timestamps) - - return efc_gens, sys_variable, sys_target - -end - -function update_firmcapacity!( - sys::SystemModel, gens::Vector{Tuple{Int,Float64}}, capacity::Int) - - for (g, share) in gens - sys.generators.capacity[g, :] .= round(Int, share * capacity) - end - -end diff --git a/PRASCapacityCredits.jl/src/EFC_Secant.jl b/PRASCapacityCredits.jl/src/EFC_Secant.jl new file mode 100644 index 00000000..52ebb03b --- /dev/null +++ b/PRASCapacityCredits.jl/src/EFC_Secant.jl @@ -0,0 +1,157 @@ +# Secant-based Equivalent Firm Capacity (EFC) Method +# +# This implementation uses a Secant-based root-finding approach to find the capacity credit. +# +# Key Advantages: +# - Speed: On benchmarking with the RTS system, the Secant method is generally 2-3 times +# faster than bisection as it uses metric gradients to approach the root more directly. +# - Informed Stepping: Unlike bisection, which reduces search space by a fixed amount, +# the Secant method takes informed steps, reducing the number of costly system assessments. +# - Robustness: More robust to different perturbation step sizes and provides a +# converged point estimate while strictly honoring user-specified gap tolerances. + +struct EFC_Secant{M} <: CapacityValuationMethod{M} + + capacity_max::Int + capacity_gap::Int + p_value::Float64 + regions::Vector{Tuple{String,Float64}} + verbose::Bool + + function EFC_Secant{M}( + capacity_max::Int, regions::Vector{Pair{String,Float64}}; + capacity_gap::Int=1, p_value::Float64=0.05, verbose::Bool=false) where M + + @assert capacity_max > 0 + @assert capacity_gap > 0 + @assert 0 < p_value < 1 + @assert sum(x.second for x in regions) ≈ 1.0 + + return new{M}(capacity_max, capacity_gap, p_value, Tuple.(regions), verbose) + + end + +end + +function EFC_Secant{M}( + capacity_max::Int, region::String; kwargs... +) where M + return EFC_Secant{M}(capacity_max, [region=>1.0]; kwargs...) +end + +function assess(sys_baseline::S, sys_augmented::S, + params::EFC_Secant{M}, simulationspec::SequentialMonteCarlo +) where {N, L, T, P, S <: SystemModel{N,L,T,P}, M <: ReliabilityMetric} + + _, powerunit, _ = unitsymbol(sys_baseline) + + regionnames = sys_baseline.regions.names + regionnames != sys_augmented.regions.names && + error("Systems provided do not have matching regions") + + # Add firm capacity generators to the relevant regions + efc_gens, sys_variable, sys_target = + add_firmcapacity(sys_baseline, sys_augmented, params.regions) + + target_metric = M(first(assess(sys_target, simulationspec, Shortfall()))) + + capacities = Int[] + metrics = typeof(target_metric)[] + + # Initial points for Secant method + # L = 0 MW + c_low = 0 + update_firmcapacity!(sys_variable, efc_gens, c_low) + m_low = M(first(assess(sys_variable, simulationspec, Shortfall()))) + f_low = val(m_low) - val(target_metric) + push!(capacities, c_low) + push!(metrics, m_low) + + # U = Max Capacity + c_high = params.capacity_max + update_firmcapacity!(sys_variable, efc_gens, c_high) + m_high = M(first(assess(sys_variable, simulationspec, Shortfall()))) + f_high = val(m_high) - val(target_metric) + push!(capacities, c_high) + push!(metrics, m_high) + + # Bracketing check/search + if f_low * f_high > 0 + # If they don't bracket, we scan to find a bracket + found_bracket = false + step = max(params.capacity_gap, 1) + + # We assume f is decreasing (EFC). If both are negative, we are already above target? + # If f_low < 0, then even at 0 MW firm capacity, augmented is better than baseline. EFC = 0. + if f_low < 0 + return CapacityCreditResult{typeof(params), typeof(target_metric), P}( + target_metric, 0, 0, capacities, metrics) + end + + # If both are positive, scan upwards + c_prev = c_low + f_prev = f_low + for c in step:step:params.capacity_max + update_firmcapacity!(sys_variable, efc_gens, c) + m_c = M(first(assess(sys_variable, simulationspec, Shortfall()))) + f_c = val(m_c) - val(target_metric) + push!(capacities, c) + push!(metrics, m_c) + + if f_c * f_prev <= 0 + c_low, c_high = c_prev, c + f_low, f_high = f_prev, f_c + found_bracket = true + break + end + c_prev, f_prev = c, f_c + end + + if !found_bracket + # If still not found, return the best we have (either 0 or max) + final_val = f_high > 0 ? c_high : c_low + return CapacityCreditResult{typeof(params), typeof(target_metric), P}( + target_metric, final_val, final_val, capacities, metrics) + end + end + + # Robust Secant (False Position / Regula Falsi) loop + iter = 0 + max_iter = 100 + + while (c_high - c_low) > params.capacity_gap && iter < max_iter + iter += 1 + + # Secant/Interpolation guess + if abs(f_high - f_low) < 1e-12 + c_mid = div(c_low + c_high, 2) + else + c_mid_float = c_high - f_high * (c_high - c_low) / (f_high - f_low) + c_mid = round(Int, c_mid_float) + + # Ensure we are making progress and staying in the bracket + if c_mid <= c_low || c_mid >= c_high + c_mid = div(c_low + c_high, 2) + end + end + + update_firmcapacity!(sys_variable, efc_gens, c_mid) + m_mid = M(first(assess(sys_variable, simulationspec, Shortfall()))) + f_mid = val(m_mid) - val(target_metric) + push!(capacities, c_mid) + push!(metrics, m_mid) + + if f_mid * f_low > 0 + c_low, f_low = c_mid, f_mid + else + c_high, f_high = c_mid, f_mid + end + end + + # Return the converged value (conservative lower bound of the bracket) + final_val = c_low + + return CapacityCreditResult{typeof(params), typeof(target_metric), P}( + target_metric, final_val, final_val, capacities, metrics) + +end diff --git a/PRASCapacityCredits.jl/src/ELCC.jl b/PRASCapacityCredits.jl/src/ELCC.jl index 3b2cd0c6..89da3cf4 100644 --- a/PRASCapacityCredits.jl/src/ELCC.jl +++ b/PRASCapacityCredits.jl/src/ELCC.jl @@ -106,33 +106,4 @@ function assess(sys_baseline::S, sys_augmented::S, end -function copy_load( - sys::SystemModel{N,L,T,P,E}, - region_shares::Vector{Tuple{String,Float64}} -) where {N,L,T,P,E} - region_allocations = allocate_regions(sys.regions.names, region_shares) - - new_regions = Regions{N,P}(sys.regions.names, copy(sys.regions.load)) - - return region_allocations, sys.regions.load, SystemModel( - new_regions, sys.interfaces, - sys.generators, sys.region_gen_idxs, - sys.storages, sys.region_stor_idxs, - sys.generatorstorages, sys.region_genstor_idxs, - sys.lines, sys.interface_line_idxs, sys.timestamps) - -end - -function update_load!( - sys::SystemModel, - region_shares::Vector{Tuple{Int,Float64}}, - load_base::Matrix{Int}, - load_increase::Int -) - for (r, share) in region_shares - sys.regions.load[r, :] .= load_base[r, :] .+ - round(Int, share * load_increase) - end - -end diff --git a/PRASCapacityCredits.jl/src/ELCC_Secant.jl b/PRASCapacityCredits.jl/src/ELCC_Secant.jl new file mode 100644 index 00000000..98758947 --- /dev/null +++ b/PRASCapacityCredits.jl/src/ELCC_Secant.jl @@ -0,0 +1,157 @@ +# Secant-based Effective Load Carrying Capability (ELCC) Method +# +# This implementation uses a Secant-based root-finding approach to find the capacity credit. +# +# Key Advantages: +# - Speed: On benchmarking with the RTS system, the Secant method is generally 2-3 times +# faster than bisection as it uses metric gradients to approach the root more directly. +# - Informed Stepping: Unlike bisection, which reduces search space by a fixed amount, +# the Secant method takes informed steps, reducing the number of costly system assessments. +# - Robustness: More robust to different perturbation step sizes and provides a +# converged point estimate while strictly honoring user-specified gap tolerances. + +struct ELCC_Secant{M} <: CapacityValuationMethod{M} + + capacity_max::Int + capacity_gap::Int + p_value::Float64 + regions::Vector{Tuple{String,Float64}} + verbose::Bool + + function ELCC_Secant{M}( + capacity_max::Int, regions::Vector{Pair{String,Float64}}; + capacity_gap::Int=1, p_value::Float64=0.05, verbose::Bool=false) where M + + @assert capacity_max > 0 + @assert capacity_gap > 0 + @assert 0 < p_value < 1 + @assert sum(x.second for x in regions) ≈ 1.0 + + return new{M}(capacity_max, capacity_gap, p_value, Tuple.(regions), verbose) + + end + +end + +function ELCC_Secant{M}( + capacity_max::Int, region::String; kwargs... +) where M + return ELCC_Secant{M}(capacity_max, [region=>1.0]; kwargs...) +end + +function assess(sys_baseline::S, sys_augmented::S, + params::ELCC_Secant{M}, simulationspec::SequentialMonteCarlo +) where {N, L, T, P, S <: SystemModel{N,L,T,P}, M <: ReliabilityMetric} + + _, powerunit, _ = unitsymbol(sys_baseline) + + regionnames = sys_baseline.regions.names + regionnames != sys_augmented.regions.names && + error("Systems provided do not have matching regions") + + target_metric = M(first(assess(sys_baseline, simulationspec, Shortfall()))) + + capacities = Int[] + metrics = typeof(target_metric)[] + + elcc_regions, base_load, sys_variable = + copy_load(sys_augmented, params.regions) + + # Initial points for Secant method + # L = 0 MW + c_low = 0 + update_load!(sys_variable, elcc_regions, base_load, c_low) + m_low = M(first(assess(sys_variable, simulationspec, Shortfall()))) + f_low = val(m_low) - val(target_metric) + push!(capacities, c_low) + push!(metrics, m_low) + + # U = Max Capacity + c_high = params.capacity_max + update_load!(sys_variable, elcc_regions, base_load, c_high) + m_high = M(first(assess(sys_variable, simulationspec, Shortfall()))) + f_high = val(m_high) - val(target_metric) + push!(capacities, c_high) + push!(metrics, m_high) + + # Bracketing check/search + if f_low * f_high > 0 + # If they don't bracket, we scan to find a bracket + found_bracket = false + step = max(params.capacity_gap, 1) + + # We assume f is increasing (ELCC). If both are positive, we are already above target? + # But usually at 0 load, EUE_augmented < EUE_baseline, so f_low < 0. + # If f_low > 0, then even at 0 MW added load, augmented is worse than baseline. ELCC = 0. + if f_low > 0 + return CapacityCreditResult{typeof(params), typeof(target_metric), P}( + target_metric, 0, 0, capacities, metrics) + end + + # If both are negative, scan upwards + c_prev = c_low + f_prev = f_low + for c in step:step:params.capacity_max + update_load!(sys_variable, elcc_regions, base_load, c) + m_c = M(first(assess(sys_variable, simulationspec, Shortfall()))) + f_c = val(m_c) - val(target_metric) + push!(capacities, c) + push!(metrics, m_c) + + if f_c * f_prev <= 0 + c_low, c_high = c_prev, c + f_low, f_high = f_prev, f_c + found_bracket = true + break + end + c_prev, f_prev = c, f_c + end + + if !found_bracket + # If still not found, return the best we have (either 0 or max) + final_val = f_high < 0 ? c_high : c_low + return CapacityCreditResult{typeof(params), typeof(target_metric), P}( + target_metric, final_val, final_val, capacities, metrics) + end + end + + # Robust Secant (False Position / Regula Falsi) loop + iter = 0 + max_iter = 100 + + while (c_high - c_low) > params.capacity_gap && iter < max_iter + iter += 1 + + # Secant/Interpolation guess + if abs(f_high - f_low) < 1e-12 + c_mid = div(c_low + c_high, 2) + else + c_mid_float = c_high - f_high * (c_high - c_low) / (f_high - f_low) + c_mid = round(Int, c_mid_float) + + # Ensure we are making progress and staying in the bracket + if c_mid <= c_low || c_mid >= c_high + c_mid = div(c_low + c_high, 2) + end + end + + update_load!(sys_variable, elcc_regions, base_load, c_mid) + m_mid = M(first(assess(sys_variable, simulationspec, Shortfall()))) + f_mid = val(m_mid) - val(target_metric) + push!(capacities, c_mid) + push!(metrics, m_mid) + + if f_mid * f_low > 0 + c_low, f_low = c_mid, f_mid + else + c_high, f_high = c_mid, f_mid + end + end + + # Return the converged value (conservative lower bound of the bracket) + final_val = c_low + + return CapacityCreditResult{typeof(params), typeof(target_metric), P}( + target_metric, final_val, final_val, capacities, metrics) + +end diff --git a/PRASCapacityCredits.jl/src/PRASCapacityCredits.jl b/PRASCapacityCredits.jl/src/PRASCapacityCredits.jl index cad35694..feee810b 100644 --- a/PRASCapacityCredits.jl/src/PRASCapacityCredits.jl +++ b/PRASCapacityCredits.jl/src/PRASCapacityCredits.jl @@ -7,13 +7,15 @@ import PRASCore.Results: ReliabilityMetric, Result, Shortfall, stderror, val import Base: minimum, maximum, extrema import Distributions: ccdf, Normal -export EFC, ELCC +export EFC, EFC_Secant, ELCC, ELCC_Secant abstract type CapacityValuationMethod{M<:ReliabilityMetric} end include("utils.jl") include("CapacityCreditResult.jl") include("EFC.jl") +include("EFC_Secant.jl") include("ELCC.jl") +include("ELCC_Secant.jl") end diff --git a/PRASCapacityCredits.jl/src/utils.jl b/PRASCapacityCredits.jl/src/utils.jl index 2f09a4a8..f23f913b 100644 --- a/PRASCapacityCredits.jl/src/utils.jl +++ b/PRASCapacityCredits.jl/src/utils.jl @@ -43,3 +43,115 @@ end incr_range(rnge::UnitRange{Int}, inc::Int) = rnge .+ inc incr_range(rnge::UnitRange{Int}, inc1::Int, inc2::Int) = (first(rnge) + inc1):(last(rnge) + inc2) + +function copy_load( + sys::SystemModel{N,L,T,P,E}, + region_shares::Vector{Tuple{String,Float64}} +) where {N,L,T,P,E} + + region_allocations = allocate_regions(sys.regions.names, region_shares) + + new_regions = Regions{N,P}(sys.regions.names, copy(sys.regions.load)) + + return region_allocations, sys.regions.load, SystemModel( + new_regions, sys.interfaces, + sys.generators, sys.region_gen_idxs, + sys.storages, sys.region_stor_idxs, + sys.generatorstorages, sys.region_genstor_idxs, + sys.lines, sys.interface_line_idxs, sys.timestamps) + +end + +function update_load!( + sys::SystemModel, + region_shares::Vector{Tuple{Int,Float64}}, + load_base::Matrix{Int}, + load_increase::Int +) + for (r, share) in region_shares + sys.regions.load[r, :] .= load_base[r, :] .+ + round(Int, share * load_increase) + end + +end + +function add_firmcapacity( + s1::SystemModel{N,L,T,P,E}, s2::SystemModel{N,L,T,P,E}, + region_shares::Vector{Tuple{String,Float64}} +) where {N,L,T,P,E} + + n_regions = length(s1.regions.names) + n_region_allocs = length(region_shares) + + region_allocations = allocate_regions(s1.regions.names, region_shares) + efc_gens = similar(region_allocations) + + new_gen(i::Int) = Generators{N,L,T,P}( + ["_EFC_$i"], ["_EFC Calculation Dummy Generator"], + zeros(Int, 1, N), zeros(1, N), ones(1, N)) + + variable_gens = Generators{N,L,T,P}[] + variable_region_gen_idxs = similar(s1.region_gen_idxs) + + target_gens = similar(variable_gens) + target_region_gen_idxs = similar(s2.region_gen_idxs) + + ra_idx = 0 + + for r in 1:n_regions + + s1_range = s1.region_gen_idxs[r] + s2_range = s2.region_gen_idxs[r] + + if (ra_idx < n_region_allocs) && (r == first(region_allocations[ra_idx+1])) + + ra_idx += 1 + + variable_region_gen_idxs[r] = incr_range(s1_range, ra_idx-1, ra_idx) + target_region_gen_idxs[r] = incr_range(s2_range, ra_idx-1, ra_idx) + + gen = new_gen(ra_idx) + push!(variable_gens, gen) + push!(target_gens, gen) + efc_gens[ra_idx] = ( + first(s1_range) + ra_idx - 1, + last(region_allocations[ra_idx])) + + else + + variable_region_gen_idxs[r] = incr_range(s1_range, ra_idx) + target_region_gen_idxs[r] = incr_range(s2_range, ra_idx) + + end + + push!(variable_gens, s1.generators[s1_range]) + push!(target_gens, s2.generators[s2_range]) + + end + + sys_variable = SystemModel( + s1.regions, s1.interfaces, + vcat(variable_gens...), variable_region_gen_idxs, + s1.storages, s1.region_stor_idxs, + s1.generatorstorages, s1.region_genstor_idxs, + s1.lines, s1.interface_line_idxs, s1.timestamps) + + sys_target = SystemModel( + s2.regions, s2.interfaces, + vcat(target_gens...), target_region_gen_idxs, + s2.storages, s2.region_stor_idxs, + s2.generatorstorages, s2.region_genstor_idxs, + s2.lines, s2.interface_line_idxs, s2.timestamps) + + return efc_gens, sys_variable, sys_target + +end + +function update_firmcapacity!( + sys::SystemModel, gens::Vector{Tuple{Int,Float64}}, capacity::Int) + + for (g, share) in gens + sys.generators.capacity[g, :] .= round(Int, share * capacity) + end + +end diff --git a/PRASCapacityCredits.jl/test/runtests.jl b/PRASCapacityCredits.jl/test/runtests.jl index 9c9d936f..32770892 100644 --- a/PRASCapacityCredits.jl/test/runtests.jl +++ b/PRASCapacityCredits.jl/test/runtests.jl @@ -61,6 +61,18 @@ import PRASCore.Systems: TestData end + @testset "EFC_Secant" begin + + cc = assess(sys_before, sys_after, EFC_Secant{EUE}(10, "Region"), smc) + # Bisection found [8, 9], Secant returns a point (bound_low == bound_high) + @test 8 <= minimum(cc) <= 9 + + cc = assess(TestData.threenode, threenode2, EFC_Secant{EUE}(10, "Region A"), smc) + # Bisection found [3, 4] + @test 3 <= minimum(cc) <= 4 + + end + @testset "ELCC" begin cc = assess(sys_before, sys_after, ELCC{EUE}(10, "Region"), smc) @@ -71,4 +83,16 @@ import PRASCore.Systems: TestData end + @testset "ELCC_Secant" begin + + cc = assess(sys_before, sys_after, ELCC_Secant{EUE}(10, "Region"), smc) + # Bisection found [7, 8] + @test 7 <= minimum(cc) <= 8 + + cc = assess(TestData.threenode, threenode2, ELCC_Secant{EUE}(10, "Region A"), smc) + # Bisection found [3, 4] + @test 3 <= minimum(cc) <= 4 + + end + end diff --git a/docs/src/quickstart.md b/docs/src/quickstart.md index 5ef5dea2..12a8def1 100644 --- a/docs/src/quickstart.md +++ b/docs/src/quickstart.md @@ -144,6 +144,18 @@ base_system # The base system augmented with some incremental resource of interest augmented_system +You can also choose the solution method for ELCC and EFC calculation. The default is a bisection method (EFC, ELCC), but a secant method (EFC_Secant, ELCC_Secant) is also available. + +```julia +# Bisection method (default) +elcc_bisection = assess(base_system, augmented_system, ELCC{EUE}(200, "1"), SequentialMonteCarlo(samples=10,seed=1)) + +# Secant method +elcc_secant = assess(base_system, augmented_system, ELCC_Secant{EUE}(200, "1"), SequentialMonteCarlo(samples=10,seed=1)) +``` + +```julia +Note that in the `rts_gmlc_sys` example, the region is identified by the string `"1"`, which corresponds to region 1 in that dataset. # Get the lower and upper bounds on the ELCC estimate for the resource elcc = assess( base_system, augmented_system, ELCC{EUE}(1000, "A"),