add Secant-based methods for ELCC and EFC calculations

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

PRASCapacityCredits.jl/src/EFC_Secant.jl

@@ -0,0 +1,123 @@
+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
+    # Point 0: 0 MW
+    c_prev = 0
+    update_firmcapacity!(sys_variable, efc_gens, c_prev)
+    m_prev = M(first(assess(sys_variable, simulationspec, Shortfall())))
+    push!(capacities, c_prev)
+    push!(metrics, m_prev)
+
+    # Point 1: Max Capacity
+    c_curr = params.capacity_max
+    update_firmcapacity!(sys_variable, efc_gens, c_curr)
+    m_curr = M(first(assess(sys_variable, simulationspec, Shortfall())))
+    push!(capacities, c_curr)
+    push!(metrics, m_curr)
+
+    iter = 0
+    max_iter = 100 
+
+    final_val = c_curr
+
+    while iter < max_iter
+        iter += 1
+
+        params.verbose && println(
+            "Iteration $iter: c_prev=$c_prev, c_curr=$c_curr, m_prev=$(val(m_prev)), m_curr=$(val(m_curr)), target=$(val(target_metric))"
+        )
+
+        # g(c) = Metric(c) - Target
+        g_prev = val(m_prev) - val(target_metric)
+        g_curr = val(m_curr) - val(target_metric)
+
+        if abs(g_curr - g_prev) < 1e-9
+            params.verbose && @info "Denominator too small in Secant method, stopping."
+            final_val = c_curr
+            break
+        end
+
+        # Secant update
+        c_next_float = c_curr - g_curr * (c_curr - c_prev) / (g_curr - g_prev)
+        c_next = round(Int, c_next_float)
+
+        # Check stopping criteria: Capacity gap
+        if abs(c_next - c_curr) <= params.capacity_gap
+             params.verbose && @info "Capacity change within tolerance ($(params.capacity_gap)), stopping."
+             final_val = c_next
+
+             # Evaluate final point if different
+             if c_next != c_curr
+                update_firmcapacity!(sys_variable, efc_gens, c_next)
+                m_next = M(first(assess(sys_variable, simulationspec, Shortfall())))
+                push!(capacities, c_next)
+                push!(metrics, m_next)
+             end
+             break
+        end
+
+        # Evaluate new point
+        update_firmcapacity!(sys_variable, efc_gens, c_next)
+        m_next = M(first(assess(sys_variable, simulationspec, Shortfall())))
+        push!(capacities, c_next)
+        push!(metrics, m_next)
+
+        # Setup for next iteration
+        c_prev = c_curr
+        m_prev = m_curr
+        c_curr = c_next
+        m_curr = m_next
+        final_val = c_curr
+
+    end
+
+    return CapacityCreditResult{typeof(params), typeof(target_metric), P}(
+        target_metric, final_val, final_val, capacities, metrics)
+
+end

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

PRASCapacityCredits.jl/src/ELCC_Secant.jl

@@ -0,0 +1,122 @@
+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
+    # Point 0: 0 MW
+    c_prev = 0
+    update_load!(sys_variable, elcc_regions, base_load, c_prev)
+    m_prev = M(first(assess(sys_variable, simulationspec, Shortfall())))
+    push!(capacities, c_prev)
+    push!(metrics, m_prev)
+
+    # Point 1: Max Capacity
+    c_curr = params.capacity_max
+    update_load!(sys_variable, elcc_regions, base_load, c_curr)
+    m_curr = M(first(assess(sys_variable, simulationspec, Shortfall())))
+    push!(capacities, c_curr)
+    push!(metrics, m_curr)
+
+    iter = 0
+    max_iter = 100 
+
+    final_val = c_curr
+
+    while iter < max_iter
+        iter += 1
+
+        params.verbose && println(
+            "Iteration $iter: c_prev=$c_prev, c_curr=$c_curr, m_prev=$(val(m_prev)), m_curr=$(val(m_curr)), target=$(val(target_metric))"
+        )
+
+        # g(c) = Metric(c) - Target
+        g_prev = val(m_prev) - val(target_metric)
+        g_curr = val(m_curr) - val(target_metric)
+
+        if abs(g_curr - g_prev) < 1e-9
+            params.verbose && @info "Denominator too small in Secant method, stopping."
+            final_val = c_curr
+            break
+        end
+
+        # Secant update
+        c_next_float = c_curr - g_curr * (c_curr - c_prev) / (g_curr - g_prev)
+        c_next = round(Int, c_next_float)
+
+        # Check stopping criteria: Capacity gap
+        if abs(c_next - c_curr) <= params.capacity_gap
+             params.verbose && @info "Capacity change within tolerance ($(params.capacity_gap)), stopping."
+             final_val = c_next
+
+             # Evaluate final point if different
+             if c_next != c_curr
+                update_load!(sys_variable, elcc_regions, base_load, c_next)
+                m_next = M(first(assess(sys_variable, simulationspec, Shortfall())))
+                push!(capacities, c_next)
+                push!(metrics, m_next)
+             end
+             break
+        end
+
+        # Evaluate new point
+        update_load!(sys_variable, elcc_regions, base_load, c_next)
+        m_next = M(first(assess(sys_variable, simulationspec, Shortfall())))
+        push!(capacities, c_next)
+        push!(metrics, m_next)
+
+        # Setup for next iteration
+        c_prev = c_curr
+        m_prev = m_curr
+        c_curr = c_next
+        m_curr = m_next
+        final_val = c_curr
+
+    end
+
+    return CapacityCreditResult{typeof(params), typeof(target_metric), P}(
+        target_metric, final_val, final_val, capacities, metrics)
+
+end

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