From dd6ad3905b2372166b864771b547945b9f401868 Mon Sep 17 00:00:00 2001 From: akrivi <20168326+akrivi@users.noreply.github.com> Date: Thu, 30 Apr 2026 17:14:11 -0600 Subject: [PATCH 01/13] shortfall events implementation (originally in PR#103) --- PRASCore.jl/src/Results/Results.jl | 8 +- PRASCore.jl/src/Results/ShortfallEvents.jl | 240 +++++++++++++++++++++ PRASCore.jl/src/Results/metrics.jl | 145 ++++++++++++- PRASCore.jl/src/Simulations/recording.jl | 103 ++++++++- PRASCore.jl/src/Simulations/utils.jl | 3 +- docs/src/PRASCore/api.md | 6 + 6 files changed, 500 insertions(+), 5 deletions(-) create mode 100644 PRASCore.jl/src/Results/ShortfallEvents.jl diff --git a/PRASCore.jl/src/Results/Results.jl b/PRASCore.jl/src/Results/Results.jl index 059bae38..b7ba96bc 100644 --- a/PRASCore.jl/src/Results/Results.jl +++ b/PRASCore.jl/src/Results/Results.jl @@ -13,7 +13,8 @@ import ..Systems: SystemModel, ZonedDateTime, Period, export # Metrics - ReliabilityMetric, LOLE, EUE, NEUE, LOLD, + ReliabilityMetric, LOLE, EUE, NEUE, LOLD, LOLEv, + MeanEventDuration, MaxEventDuration, MeanEventEnergy, MaxEventEnergy, val, stderror, CVAR, NCVAR, # Result specifications @@ -26,7 +27,7 @@ export DemandResponseEnergy, DemandResponseEnergySamples, GeneratorAvailability, StorageAvailability, GeneratorStorageAvailability,DemandResponseAvailability, - LineAvailability + LineAvailability, ShortfallEvents include("metrics.jl") include("utils.jl") @@ -188,12 +189,15 @@ getindex(x::AbstractEnergyResult, name::String, ::Colon) = getindex(x::AbstractEnergyResult, ::Colon, ::Colon) = getindex.(x, names(x), permutedims(x.timestamps)) +abstract type AbstractShortfallEventResult{N,L,T} <: Result{N,L,T} end + include("StorageEnergy.jl") include("GeneratorStorageEnergy.jl") include("DemandResponseEnergy.jl") include("StorageEnergySamples.jl") include("GeneratorStorageEnergySamples.jl") include("DemandResponseEnergySamples.jl") +include("ShortfallEvents.jl") function resultchannel( results::T, nworkers::Int diff --git a/PRASCore.jl/src/Results/ShortfallEvents.jl b/PRASCore.jl/src/Results/ShortfallEvents.jl new file mode 100644 index 00000000..7328661c --- /dev/null +++ b/PRASCore.jl/src/Results/ShortfallEvents.jl @@ -0,0 +1,240 @@ +""" + ShortfallEvents + +The `ShortfallEvents` result specification reports sample-level shortfall +events, producing a `ShortfallEventsResult`. + +A shortfall event is a contiguous run of one or more simulation timesteps +with positive shortfall. + +This result can be used to inspect event start/end times and to compute +event-based reliability metrics such as [`LOLEv`](@ref). +""" +struct ShortfallEvents <: ResultSpec end + +usesamplepartitions(::ShortfallEvents) = true + +struct ShortfallEvent + start_idx::Int + end_idx::Int + energy::Int + + function ShortfallEvent(start_idx::Int, end_idx::Int, energy::Int) + start_idx > 0 || throw(DomainError(start_idx, "start_idx must be positive")) + end_idx >= start_idx || throw(DomainError(end_idx, "end_idx must be >= start_idx")) + energy >= 0 || throw(DomainError(energy, "energy must be non-negative")) + new(start_idx, end_idx, energy) + end +end + +duration_periods(ev::ShortfallEvent) = ev.end_idx - ev.start_idx + 1 +event_energy(ev::ShortfallEvent) = ev.energy + +mutable struct ShortfallEventsAccumulator{S} <: ResultAccumulator{ShortfallEvents} + + system_events::Vector{Vector{ShortfallEvent}} + region_events::Matrix{Vector{ShortfallEvent}} + + in_system_event::Bool + system_event_start::Int + system_event_energy::Int + + in_region_event::Vector{Bool} + region_event_start::Vector{Int} + region_event_energy::Vector{Int} + + nperiods::Int +end + +function accumulator( + sys::SystemModel{N}, nsamples::Int, ::S +) where {N,S<:ShortfallEvents} + + nregions = length(sys.regions) + + system_events = [ShortfallEvent[] for _ in 1:nsamples] + region_events = [ShortfallEvent[] for _ in 1:nregions, _ in 1:nsamples] + + in_system_event = false + system_event_start = 0 + system_event_energy = 0 + + in_region_event = falses(nregions) + region_event_start = zeros(Int, nregions) + region_event_energy = zeros(Int, nregions) + + return ShortfallEventsAccumulator{S}( + system_events, region_events, + in_system_event, system_event_start, system_event_energy, + in_region_event, region_event_start, region_event_energy, + N) +end + +function merge!( + x::ShortfallEventsAccumulator, y::ShortfallEventsAccumulator +) + foreach(append!, x.system_events, y.system_events) + foreach(append!, x.region_events, y.region_events) + return +end + +function copy_sample_partition!( + x::ShortfallEventsAccumulator, + y::ShortfallEventsAccumulator, + sampleids::UnitRange{Int}, +) + @views x.system_events[sampleids] .= y.system_events + @views x.region_events[:, sampleids] .= y.region_events + return +end + +accumulatortype(::S) where { + S<:ShortfallEvents + } = ShortfallEventsAccumulator{S} + +struct ShortfallEventsResult{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit,S} <: AbstractShortfallEventResult{N,L,T} + regions::Regions + timestamps::StepRange{ZonedDateTime,T} + + system_events::Vector{Vector{ShortfallEvent}} + region_events::Matrix{Vector{ShortfallEvent}} +end + +""" + getindex(x::ShortfallEventsResult, s::Int) + +Return the vector of system-wide shortfall events for sample `s`. +""" +function getindex(x::ShortfallEventsResult, s::Int) + return x.system_events[s] +end + +""" + getindex(x::ShortfallEventsResult, r::AbstractString) + +Return a vector whose `s`th entry is the vector of shortfall events +for region `r` in sample `s`. +""" +function getindex(x::ShortfallEventsResult, r::AbstractString) + i_r = findfirstunique(x.regions.names, r) + return [x.region_events[i_r, s] for s in axes(x.region_events, 2)] +end + +""" + getindex(x::ShortfallEventsResult, r::AbstractString, s::Int) + +Return the vector of shortfall events for region `r` in sample `s`. +""" +function getindex(x::ShortfallEventsResult, r::AbstractString, s::Int) + i_r = findfirstunique(x.regions.names, r) + return x.region_events[i_r, s] +end + +start_event_timestamp(x::ShortfallEventsResult, ev::ShortfallEvent) = x.timestamps[ev.start_idx] +end_event_timestamp(x::ShortfallEventsResult, ev::ShortfallEvent) = x.timestamps[ev.end_idx] + +LOLEv(x::ShortfallEventsResult{N,L,T}) where {N,L,T} = + LOLEv{N,L,T}(MeanEstimate(length.(x.system_events))) + +function LOLEv(x::ShortfallEventsResult{N,L,T}, r::AbstractString) where {N,L,T} + i_r = findfirstunique(x.regions.names, r) + counts = [length(x.region_events[i_r, s]) for s in axes(x.region_events, 2)] + return LOLEv{N,L,T}(MeanEstimate(counts)) +end + +function finalize( + acc::ShortfallEventsAccumulator{S}, + system::SystemModel{N,L,T,P,E}, +) where {N,L,T,P,E,S<:ShortfallEvents} + + + return ShortfallEventsResult{N,L,T,P,E,S}( + system.regions, system.timestamps, + acc.system_events, acc.region_events) +end + +function MeanEventDuration(x::ShortfallEventsResult{N,L,T}) where {N,L,T} + durations = Float64[ + isempty(events) ? 0.0 : mean(duration_periods.(events)) + for events in x.system_events + ] + return MeanEventDuration{N,L,T}(MeanEstimate(durations)) +end + +function MeanEventDuration(x::ShortfallEventsResult{N,L,T}, r::AbstractString) where {N,L,T} + i_r = findfirstunique(x.regions.names, r) + durations = Float64[ + isempty(x.region_events[i_r, s]) ? 0.0 : + mean(duration_periods.(x.region_events[i_r, s])) + for s in axes(x.region_events, 2) + ] + return MeanEventDuration{N,L,T}(MeanEstimate(durations)) +end + +function MaxEventDuration(x::ShortfallEventsResult{N,L,T}) where {N,L,T} + durations = Float64[ + isempty(events) ? 0.0 : maximum(duration_periods.(events)) + for events in x.system_events + ] + return MaxEventDuration{N,L,T}(MeanEstimate(durations)) +end + +function MaxEventDuration(x::ShortfallEventsResult{N,L,T}, r::AbstractString) where {N,L,T} + i_r = findfirstunique(x.regions.names, r) + durations = Float64[ + isempty(x.region_events[i_r, s]) ? 0.0 : + maximum(duration_periods.(x.region_events[i_r, s])) + for s in axes(x.region_events, 2) + ] + return MaxEventDuration{N,L,T}(MeanEstimate(durations)) +end + +function MeanEventEnergy(x::ShortfallEventsResult{N,L,T,P,E}) where {N,L,T,P,E} + p2e = conversionfactor(L, T, P, E) + energies = Float64[ + isempty(events) ? 0.0 : mean(p2e .* event_energy.(events)) + for events in x.system_events + ] + return MeanEventEnergy{N,L,T,E}(MeanEstimate(energies)) +end + +function MeanEventEnergy(x::ShortfallEventsResult{N,L,T,P,E}, r::AbstractString) where {N,L,T,P,E} + i_r = findfirstunique(x.regions.names, r) + p2e = conversionfactor(L, T, P, E) + energies = Float64[ + isempty(x.region_events[i_r, s]) ? 0.0 : + mean(p2e .* event_energy.(x.region_events[i_r, s])) + for s in axes(x.region_events, 2) + ] + return MeanEventEnergy{N,L,T,E}(MeanEstimate(energies)) +end + +function MaxEventEnergy(x::ShortfallEventsResult{N,L,T,P,E}) where {N,L,T,P,E} + p2e = conversionfactor(L, T, P, E) + energies = Float64[ + isempty(events) ? 0.0 : maximum(p2e .* event_energy.(events)) + for events in x.system_events + ] + return MaxEventEnergy{N,L,T,E}(MeanEstimate(energies)) +end + +function MaxEventEnergy(x::ShortfallEventsResult{N,L,T,P,E}, r::AbstractString) where {N,L,T,P,E} + i_r = findfirstunique(x.regions.names, r) + p2e = conversionfactor(L, T, P, E) + energies = Float64[ + isempty(x.region_events[i_r, s]) ? 0.0 : + maximum(p2e .* event_energy.(x.region_events[i_r, s])) + for s in axes(x.region_events, 2) + ] + return MaxEventEnergy{N,L,T,E}(MeanEstimate(energies)) +end + + +function totalevents(x::ShortfallEventsResult) + return sum(length, x.system_events) +end + +function totalevents(x::ShortfallEventsResult, r::AbstractString) + i_r = findfirstunique(x.regions.names, r) + return sum(length, view(x.region_events, i_r, :)) +end diff --git a/PRASCore.jl/src/Results/metrics.jl b/PRASCore.jl/src/Results/metrics.jl index da9b8df6..5f13fc97 100644 --- a/PRASCore.jl/src/Results/metrics.jl +++ b/PRASCore.jl/src/Results/metrics.jl @@ -217,6 +217,149 @@ function Base.show(io::IO, x::LOLD{D}) where {D} D == 1 ? "day" : string(D) * "days") end +""" + LOLEv + +`LOLEv` reports loss of load events over a particular time period +and regional extent. + +Contains both the estimated value itself as well as the standard error +of that estimate, which can be extracted with `val` and `stderror`, +respectively. +""" +struct LOLEv{N, L, T <: Period} <: ReliabilityMetric + lolev::MeanEstimate + + function LOLEv{N,L,T}(lolev::MeanEstimate) where {N,L,T<:Period} + val(lolev) >= 0 || throw(DomainError(val(lolev), + "$(val(lolev)) is not a valid expected count of events")) + new{N,L,T}(lolev) + end +end + +val(x::LOLEv) = val(x.lolev) +stderror(x::LOLEv) = stderror(x.lolev) + +function Base.show(io::IO, x::LOLEv{N,L,T}) where {N,L,T} + print(io, "LOLEv = ", x.lolev, " events") +end + + +""" + MeanEventDuration + +`MeanEventDuration` reports the mean duration across all observed shortfall events. +If no events are observed, the value is reported as zero. + +Contains both the estimated value itself as well as the standard error +of that estimate, which can be extracted with `val` and `stderror`, +respectively. +""" +struct MeanEventDuration{N, L, T <: Period} <: ReliabilityMetric + duration::MeanEstimate + + function MeanEventDuration{N,L,T}(duration::MeanEstimate) where {N,L,T<:Period} + val(duration) >= 0 || throw(DomainError(val(duration), + "$(val(duration)) is not a valid expected event duration")) + new{N,L,T}(duration) + end +end + +val(x::MeanEventDuration) = val(x.duration) +stderror(x::MeanEventDuration) = stderror(x.duration) + +function Base.show(io::IO, x::MeanEventDuration{N,L,T}) where {N,L,T} + t_symbol = unitsymbol(T) + print(io, "MeanEventDuration = ", x.duration, " ", + L == 1 ? t_symbol : "(" * string(L) * t_symbol * ")") +end + + +""" + MaxEventDuration + +`MaxEventDuration` reports the maximum duration across all observed shortfall events. +If no events are observed, the value is reported as zero. + +Contains both the estimated value itself as well as the standard error +of that estimate, which can be extracted with `val` and `stderror`, +respectively. +""" +struct MaxEventDuration{N, L, T <: Period} <: ReliabilityMetric + duration::MeanEstimate + + function MaxEventDuration{N,L,T}(duration::MeanEstimate) where {N,L,T<:Period} + val(duration) >= 0 || throw(DomainError(val(duration), + "$(val(duration)) is not a valid expected maximum event duration")) + new{N,L,T}(duration) + end +end + +val(x::MaxEventDuration) = val(x.duration) +stderror(x::MaxEventDuration) = stderror(x.duration) + +function Base.show(io::IO, x::MaxEventDuration{N,L,T}) where {N,L,T} + t_symbol = unitsymbol(T) + print(io, "MaxEventDuration = ", x.duration, " ", + L == 1 ? t_symbol : "(" * string(L) * t_symbol * ")") +end + + +""" + MeanEventEnergy + +`MeanEventEnergy` reports the mean unserved energy across all observed shortfall events. +If no events are observed, the value is reported as zero. + +Contains both the estimated value itself as well as the standard error +of that estimate, which can be extracted with `val` and `stderror`, +respectively. +""" +struct MeanEventEnergy{N,L,T<:Period,E<:EnergyUnit} <: ReliabilityMetric + energy::MeanEstimate + + function MeanEventEnergy{N,L,T,E}(energy::MeanEstimate) where {N,L,T<:Period,E<:EnergyUnit} + val(energy) >= 0 || throw(DomainError(val(energy), + "$(val(energy)) is not a valid expected event energy")) + new{N,L,T,E}(energy) + end +end + +val(x::MeanEventEnergy) = val(x.energy) +stderror(x::MeanEventEnergy) = stderror(x.energy) + +function Base.show(io::IO, x::MeanEventEnergy{N,L,T,E}) where {N,L,T,E} + print(io, "MeanEventEnergy = ", x.energy, " ", unitsymbol(E)) +end + + +""" + MaxEventEnergy + +`MaxEventEnergy` reports the maximum unserved energy across all observed shortfall events. +If no events are observed, the value is reported as zero. + +Contains both the estimated value itself as well as the standard error +of that estimate, which can be extracted with `val` and `stderror`, +respectively. +""" +struct MaxEventEnergy{N,L,T<:Period,E<:EnergyUnit} <: ReliabilityMetric + energy::MeanEstimate + + function MaxEventEnergy{N,L,T,E}(energy::MeanEstimate) where {N,L,T<:Period,E<:EnergyUnit} + val(energy) >= 0 || throw(DomainError(val(energy), + "$(val(energy)) is not a valid expected maximum event energy")) + new{N,L,T,E}(energy) + end +end + +val(x::MaxEventEnergy) = val(x.energy) +stderror(x::MaxEventEnergy) = stderror(x.energy) + +function Base.show(io::IO, x::MaxEventEnergy{N,L,T,E}) where {N,L,T,E} + print(io, "MaxEventEnergy = ", x.energy, " ", unitsymbol(E)) +end + const CVAR_QUANTITIES = (:energy,) _cvar_quantity_unitsymbol(::Val{:energy}, ::Type{E}, ::Type) where {E<:EnergyUnit} = unitsymbol(E) @@ -273,7 +416,7 @@ struct NCVAR <: ReliabilityMetric ncvar::MeanEstimate alpha::Float64 var::Float64 - + function NCVAR(quantity::Symbol, ncvar::MeanEstimate, alpha::Float64, var::Float64) val(ncvar) >= 0 || throw(DomainError(val(ncvar), diff --git a/PRASCore.jl/src/Simulations/recording.jl b/PRASCore.jl/src/Simulations/recording.jl index 44f6ff7d..5429523c 100644 --- a/PRASCore.jl/src/Simulations/recording.jl +++ b/PRASCore.jl/src/Simulations/recording.jl @@ -512,4 +512,105 @@ function record!( end -reset!(acc::Results.DemandResponseEnergySamplesAccumulator, sampleid::Int) = nothing \ No newline at end of file +reset!(acc::Results.DemandResponseEnergySamplesAccumulator, sampleid::Int) = nothing + +# ShortfallEvents + +function record!( + acc::Results.ShortfallEventsAccumulator{S}, + system::SystemModel{N,L,T,P,E}, + state::SystemState, problem::DispatchProblem, + sampleid::Int, t::Int +) where {N,L,T,P,E,S} + + isshortfall = false + totalshortfall = 0 + edges = problem.fp.edges + + for (r, dr_idxs) in zip(problem.region_unserved_edges, system.region_dr_idxs) + + regionshortfall = init_regionshortfall(S, edges, r) + + dr_shortfall = 0 + for i in dr_idxs + dr_shortfall += state.drs_unservedenergy[i] + end + + regionshortfall += dr_shortfall + isregionshortfall = regionshortfall > 0 + + if isregionshortfall + totalshortfall += regionshortfall + + if !acc.in_region_event[r] + acc.in_region_event[r] = true + acc.region_event_start[r] = t + acc.region_event_energy[r] = 0 + end + + acc.region_event_energy[r] += regionshortfall + + elseif acc.in_region_event[r] + push!(acc.region_events[r, sampleid], + Results.ShortfallEvent( + acc.region_event_start[r], + t - 1, + acc.region_event_energy[r])) + acc.in_region_event[r] = false + acc.region_event_start[r] = 0 + acc.region_event_energy[r] = 0 + end + + isshortfall |= isregionshortfall + end + + if isshortfall + if !acc.in_system_event + acc.in_system_event = true + acc.system_event_start = t + acc.system_event_energy = 0 + end + acc.system_event_energy += totalshortfall + + elseif acc.in_system_event + push!(acc.system_events[sampleid], + Results.ShortfallEvent( + acc.system_event_start, + t - 1, + acc.system_event_energy)) + acc.in_system_event = false + acc.system_event_start = 0 + acc.system_event_energy = 0 + end + + return +end + +function reset!(acc::Results.ShortfallEventsAccumulator, sampleid::Int) + + if acc.in_system_event + push!(acc.system_events[sampleid], + Results.ShortfallEvent( + acc.system_event_start, + acc.nperiods, + acc.system_event_energy)) + acc.in_system_event = false + acc.system_event_start = 0 + acc.system_event_energy = 0 + end + + for r in eachindex(acc.in_region_event) + if acc.in_region_event[r] + push!(acc.region_events[r, sampleid], + Results.ShortfallEvent( + acc.region_event_start[r], + acc.nperiods, + acc.region_event_energy[r])) + acc.in_region_event[r] = false + acc.region_event_start[r] = 0 + acc.region_event_energy[r] = 0 + end + end + + return +end \ No newline at end of file diff --git a/PRASCore.jl/src/Simulations/utils.jl b/PRASCore.jl/src/Simulations/utils.jl index bad7e38a..303e4f88 100644 --- a/PRASCore.jl/src/Simulations/utils.jl +++ b/PRASCore.jl/src/Simulations/utils.jl @@ -265,7 +265,8 @@ function init_regionshortfall( edges, region) where {S <: Union{ Results.Shortfall, - Results.ShortfallSamples}} + Results.ShortfallSamples, + Results.ShortfallEvents}} return edges[region].flow end diff --git a/docs/src/PRASCore/api.md b/docs/src/PRASCore/api.md index 1686994a..bc4f6908 100644 --- a/docs/src/PRASCore/api.md +++ b/docs/src/PRASCore/api.md @@ -47,4 +47,10 @@ PRASCore.Results.CVAR PRASCore.Results.NCVAR PRASCore.Results.val PRASCore.Results.stderror +PRASCore.Results.ShortfallEvents +PRASCore.Results.LOLEv +PRASCore.Results.MeanEventDuration +PRASCore.Results.MaxEventDuration +PRASCore.Results.MeanEventEnergy +PRASCore.Results.MaxEventEnergy ``` From 35617ed395af548537929d9ac5c6ee94747fb444 Mon Sep 17 00:00:00 2001 From: akrivi <20168326+akrivi@users.noreply.github.com> Date: Fri, 1 May 2026 19:15:53 -0600 Subject: [PATCH 02/13] tests (originally in PR#103) --- PRASCore.jl/test/Simulations/runtests.jl | 59 ++++++++++++++++++++++++ 1 file changed, 59 insertions(+) diff --git a/PRASCore.jl/test/Simulations/runtests.jl b/PRASCore.jl/test/Simulations/runtests.jl index 67c55ab6..97691a00 100644 --- a/PRASCore.jl/test/Simulations/runtests.jl +++ b/PRASCore.jl/test/Simulations/runtests.jl @@ -33,6 +33,8 @@ shortfall2_1a, _, flow2_1a, util2_1a, _ = assess(TestData.singlenode_a, simspec, resultspecs...) + events_1a, = assess(TestData.singlenode_a, simspec, ShortfallEvents()) + assess(TestData.singlenode_a_5min, smallsample, resultspecs...) shortfall_1a5, _, flow_1a5, util_1a5, shortfall2_1a5, _, flow2_1a5, util2_1a5, _ = @@ -54,6 +56,8 @@ StorageEnergy(), GeneratorStorageEnergy(),DemandResponseEnergy(), StorageEnergySamples(), GeneratorStorageEnergySamples(),DemandResponseEnergySamples()) + events_3, = assess(TestData.threenode, simspec, ShortfallEvents()) + @testset "Shortfall Results" begin # Single-region system A @@ -738,6 +742,61 @@ end + @testset "Shortfall Event Metrics" begin + # Single-region system + @test val(LOLEv(events_1a)) >= 0 + @test stderror(LOLEv(events_1a)) >= 0 + @test val(MeanEventDuration(events_1a)) >= 0 + @test stderror(MeanEventDuration(events_1a)) >= 0 + + @test LOLEv(events_1a) ≈ LOLEv(events_1a, "Region") + @test MeanEventDuration(events_1a) ≈ MeanEventDuration(events_1a, "Region") + @test MaxEventDuration(events_1a) ≈ MaxEventDuration(events_1a, "Region") + @test MeanEventEnergy(events_1a) ≈ MeanEventEnergy(events_1a, "Region") + @test MaxEventEnergy(events_1a) ≈ MaxEventEnergy(events_1a, "Region") + + manual_lolev_1a = mean(length.(events_1a.system_events)) + @test isapprox(val(LOLEv(events_1a)), manual_lolev_1a; rtol=1e-10) + + manual_meandur_1a = mean([ + isempty(evts) ? 0.0 : mean(Results.duration_periods.(evts)) + for evts in events_1a.system_events + ]) + @test isapprox(val(MeanEventDuration(events_1a)), manual_meandur_1a; rtol=1e-10) + + manual_maxdur_1a = mean([ + isempty(evts) ? 0.0 : maximum(Results.duration_periods.(evts)) + for evts in events_1a.system_events + ]) + @test isapprox(val(MaxEventDuration(events_1a)), manual_maxdur_1a; rtol=1e-10) + + p2e_1a = PRASCore.Systems.conversionfactor(1, Hour, PRASCore.Systems.MW, PRASCore.Systems.MWh) + + manual_meanenergy_1a = mean([ + isempty(evts) ? 0.0 : mean(p2e_1a .* Results.event_energy.(evts)) + for evts in events_1a.system_events + ]) + @test isapprox(val(MeanEventEnergy(events_1a)), manual_meanenergy_1a; rtol=1e-10) + + manual_maxenergy_1a = mean([ + isempty(evts) ? 0.0 : maximum(p2e_1a .* Results.event_energy.(evts)) + for evts in events_1a.system_events + ]) + @test isapprox(val(MaxEventEnergy(events_1a)), manual_maxenergy_1a; rtol=1e-10) + + # Multi-region system + @test val(LOLEv(events_3)) >= 0 + @test val(MeanEventDuration(events_3)) >= 0 + @test val(LOLEv(events_3, "Region A")) >= 0 + @test val(MeanEventDuration(events_3, "Region A")) >= 0 + + @test val(LOLEv(events_3)) >= val(LOLEv(events_3, "Region A")) + + @test Results.totalevents(events_1a) >= 0 + @test Results.totalevents(events_3, "Region A") >= 0 + + end + @testset "Threaded sample result partitioning" begin simspec_serial = SequentialMonteCarlo(samples=100, seed=123, threaded=false) From 0e4033cf730996998e3e989fe45ca6b91b65885e Mon Sep 17 00:00:00 2001 From: akrivi <20168326+akrivi@users.noreply.github.com> Date: Fri, 1 May 2026 20:35:52 -0600 Subject: [PATCH 03/13] PRASFiles modifications (originally in PR#103) --- PRASFiles.jl/src/PRASFiles.jl | 6 +- PRASFiles.jl/src/Results/utils.jl | 152 ++++++++++++++++++++++++++++++ PRASFiles.jl/src/Results/write.jl | 108 +++++++++++++++++++++ 3 files changed, 265 insertions(+), 1 deletion(-) diff --git a/PRASFiles.jl/src/PRASFiles.jl b/PRASFiles.jl/src/PRASFiles.jl index 614a92f1..4d2a0420 100644 --- a/PRASFiles.jl/src/PRASFiles.jl +++ b/PRASFiles.jl/src/PRASFiles.jl @@ -7,7 +7,10 @@ import PRASCore.Systems: SystemModel, Regions, Interfaces, import PRASCore.Results: EUE, LOLE, NEUE, LOLD, ShortfallResult, ShortfallSamplesResult, - AbstractShortfallResult, Result + AbstractShortfallResult, Result, ShortfallEventsResult, + ShortfallEvent, LOLEv, totalevents, + MeanEventDuration, MaxEventDuration, + MeanEventEnergy, MaxEventEnergy import StatsBase: mean import Dates: @dateformat_str, format, now import TimeZones: ZonedDateTime @@ -21,6 +24,7 @@ import JSON3: pretty export savemodel export saveshortfall +export saveevents export read_attrs include("Systems/read.jl") diff --git a/PRASFiles.jl/src/Results/utils.jl b/PRASFiles.jl/src/Results/utils.jl index db3eb3b5..a7b8eaea 100644 --- a/PRASFiles.jl/src/Results/utils.jl +++ b/PRASFiles.jl/src/Results/utils.jl @@ -72,6 +72,79 @@ function LOLDResult(shortfall::ShortfallSamplesResult; region::Union{Nothing, St ) end +struct LOLEvResult + mean::Float64 + stderror::Float64 +end + +function LOLEvResult(events::ShortfallEventsResult; region::Union{Nothing, String} = nothing) + lolev = (region === nothing) ? LOLEv(events) : LOLEv(events, region) + return LOLEvResult( + lolev.lolev.estimate, + lolev.lolev.standarderror, + ) +end + +struct MeanEventDurationResult + mean::Float64 + stderror::Float64 +end + +function MeanEventDurationResult(events::ShortfallEventsResult; region::Union{Nothing, String} = nothing) + duration = (region === nothing) ? MeanEventDuration(events) : MeanEventDuration(events, region) + return MeanEventDurationResult( + duration.duration.estimate, + duration.duration.standarderror, + ) +end + +struct MaxEventDurationResult + mean::Float64 + stderror::Float64 +end + +function MaxEventDurationResult(events::ShortfallEventsResult; region::Union{Nothing, String} = nothing) + duration = (region === nothing) ? MaxEventDuration(events) : MaxEventDuration(events, region) + return MaxEventDurationResult( + duration.duration.estimate, + duration.duration.standarderror, + ) +end + +struct MeanEventEnergyResult + mean::Float64 + stderror::Float64 +end + +function MeanEventEnergyResult(events::ShortfallEventsResult; region::Union{Nothing, String} = nothing) + energy = (region === nothing) ? MeanEventEnergy(events) : MeanEventEnergy(events, region) + return MeanEventEnergyResult( + energy.energy.estimate, + energy.energy.standarderror, + ) +end + +struct MaxEventEnergyResult + mean::Float64 + stderror::Float64 +end + +function MaxEventEnergyResult(events::ShortfallEventsResult; region::Union{Nothing, String} = nothing) + energy = (region === nothing) ? MaxEventEnergy(events) : MaxEventEnergy(events, region) + return MaxEventEnergyResult( + energy.energy.estimate, + energy.energy.standarderror, + ) +end + +struct EventRecord + sample_id::Int64 + start_timestamp::ZonedDateTime + end_timestamp::ZonedDateTime + duration_periods::Int64 + energy::Float64 +end + struct RegionResult name::String eue::EUEResult @@ -85,6 +158,17 @@ struct RegionResult shortfall_timestamps::Vector{ZonedDateTime} end +struct RegionEventResult + name::String + lolev::LOLEvResult + mean_event_duration::MeanEventDurationResult + max_event_duration::MaxEventDurationResult + mean_event_energy::MeanEventEnergyResult + max_event_energy::MaxEventEnergyResult + total_events::Int64 + events::Vector{EventRecord} +end + struct SystemResult num_samples::Int64 type_params::TypeParams @@ -97,6 +181,21 @@ struct SystemResult region_results::Vector{RegionResult} end +struct SystemEventResult + num_samples::Int64 + type_params::TypeParams + sys_attributes::Dict{String, String} + timestamps::Vector{ZonedDateTime} + lolev::LOLEvResult + mean_event_duration::MeanEventDurationResult + max_event_duration::MaxEventDurationResult + mean_event_energy::MeanEventEnergyResult + max_event_energy::MaxEventEnergyResult + total_events::Int64 + system_events::Vector{EventRecord} + region_results::Vector{RegionEventResult} +end + function get_shortfall_mean(shortfall::ShortfallResult) return shortfall.shortfall_mean end @@ -134,6 +233,51 @@ function get_lold_result( return LOLDResult(shortfall; region = region) end +function get_eventrecords( + events_by_sample::Vector{Vector{ShortfallEvent}}, + timestamps, + p2e, +) + records = EventRecord[] + + for (sample_id, evts) in enumerate(events_by_sample) + for ev in evts + push!(records, EventRecord( + sample_id, + timestamps[ev.start_idx], + timestamps[ev.end_idx], + ev.end_idx - ev.start_idx + 1, + p2e * ev.energy, + )) + end + end + + return records +end + +function get_eventrecords( + events::ShortfallEventsResult{N,L,T,P,E}, + region::String, +) where {N,L,T,P,E} + i_r = findfirst(isequal(region), events.regions.names) + p2e = conversionfactor(L, T, P, E) + + records = EventRecord[] + for sample_id in axes(events.region_events, 2) + for ev in events.region_events[i_r, sample_id] + push!(records, EventRecord( + sample_id, + events.timestamps[ev.start_idx], + events.timestamps[ev.end_idx], + ev.end_idx - ev.start_idx + 1, + p2e * ev.energy, + )) + end + end + + return records +end + # Define structtypes for different structs defined above StructType(::Type{TypeParams}) = Struct() StructType(::Type{EUEResult}) = Struct() @@ -142,3 +286,11 @@ StructType(::Type{LOLEResult}) = Struct() StructType(::Type{LOLDResult}) = Struct() StructType(::Type{RegionResult}) = OrderedStruct() StructType(::Type{SystemResult}) = OrderedStruct() +StructType(::Type{LOLEvResult}) = Struct() +StructType(::Type{MeanEventDurationResult}) = Struct() +StructType(::Type{MaxEventDurationResult}) = Struct() +StructType(::Type{MeanEventEnergyResult}) = Struct() +StructType(::Type{MaxEventEnergyResult}) = Struct() +StructType(::Type{EventRecord}) = Struct() +StructType(::Type{RegionEventResult}) = OrderedStruct() +StructType(::Type{SystemEventResult}) = OrderedStruct() diff --git a/PRASFiles.jl/src/Results/write.jl b/PRASFiles.jl/src/Results/write.jl index 5edd45d2..f930cc97 100644 --- a/PRASFiles.jl/src/Results/write.jl +++ b/PRASFiles.jl/src/Results/write.jl @@ -104,3 +104,111 @@ function saveshortfall( error("saveshortfall is not implemented for $(typeof(shortfall))") end + +function generate_eventresult( + events::ShortfallEventsResult{N,L,T,P,E}, + pras_sys::SystemModel; + include_events::Bool = false, +) where {N,L,T,P,E} + + p2e = conversionfactor(L, T, P, E) + + system_event_records = include_events ? + get_eventrecords(events.system_events, events.timestamps, p2e) : + EventRecord[] + + region_results = RegionEventResult[] + for reg_name in pras_sys.regions.names + region_event_records = include_events ? + get_eventrecords(events, reg_name) : + EventRecord[] + + push!(region_results, + RegionEventResult( + reg_name, + LOLEvResult(events, region = reg_name), + MeanEventDurationResult(events, region = reg_name), + MaxEventDurationResult(events, region = reg_name), + MeanEventEnergyResult(events, region = reg_name), + MaxEventEnergyResult(events, region = reg_name), + totalevents(events, reg_name), + region_event_records, + ) + ) + end + + sys_result = SystemEventResult( + length(events.system_events), + TypeParams(pras_sys), + pras_sys.attrs, + collect(events.timestamps), + LOLEvResult(events), + MeanEventDurationResult(events), + MaxEventDurationResult(events), + MeanEventEnergyResult(events), + MaxEventEnergyResult(events), + totalevents(events), + system_event_records, + region_results, + ) + + return sys_result +end + +""" + saveevents( + events::ShortfallEventsResult, + pras_sys::SystemModel, + outfile::String, + ) + +Save `ShortfallEventsResult` in JSON format, including both summary event +metrics and raw event records. + +# Arguments + + - `events::ShortfallEventsResult`: PRAS shortfall events result + - `pras_sys::SystemModel`: PRAS SystemModel + - `outfile::String`: Location to save the event results + +# Returns + + - Location where the event results are exported in JSON format. + +# Keywords + + - `include_events::Bool = false`: If `true`, include full raw event records + at the system and regional levels. If `false`, only summary event metrics + and counts are exported. +""" +function saveevents( + events::ShortfallEventsResult, + pras_sys::SystemModel, + outfile::String; + include_events::Bool = false, +) + + dt_now = format(now(), "dd-u-yy-H-M-S") + export_location = joinpath(outfile, dt_now) + if !(isdir(export_location)) + mkpath(export_location) + end + + event_result = generate_eventresult(events, pras_sys; include_events = include_events) + open(joinpath(export_location, "pras_event_results.json"), "w") do io + pretty(io, event_result) + end + + @info "Successfully exported PRAS ShortfallEventsResult here: $(export_location)" + return export_location +end + +function saveevents( + events::R, + pras_sys::SystemModel, + outfile::String; + include_events::Bool = false, +) where {R <: Result} + + error("saveevents is not implemented for $(typeof(events))") +end From 218d6d1fee3b4118e49ea3d7101ddafb470a9597 Mon Sep 17 00:00:00 2001 From: akrivi <20168326+akrivi@users.noreply.github.com> Date: Sun, 3 May 2026 13:01:22 -0600 Subject: [PATCH 04/13] Switch event metrics to event-based aggregation instead of sample-based --- PRASCore.jl/src/Results/ShortfallEvents.jl | 65 +++++++++++-------- PRASCore.jl/test/Simulations/runtests.jl | 74 ++++++++++++++++++---- PRASFiles.jl/src/PRASFiles.jl | 4 +- PRASFiles.jl/src/Results/utils.jl | 2 +- PRASFiles.jl/src/Results/write.jl | 3 +- 5 files changed, 101 insertions(+), 47 deletions(-) diff --git a/PRASCore.jl/src/Results/ShortfallEvents.jl b/PRASCore.jl/src/Results/ShortfallEvents.jl index 7328661c..59959391 100644 --- a/PRASCore.jl/src/Results/ShortfallEvents.jl +++ b/PRASCore.jl/src/Results/ShortfallEvents.jl @@ -29,6 +29,11 @@ end duration_periods(ev::ShortfallEvent) = ev.end_idx - ev.start_idx + 1 event_energy(ev::ShortfallEvent) = ev.energy +_event_meanestimate(xs::AbstractVector{<:Real}) = + isempty(xs) ? MeanEstimate(0.0) : MeanEstimate(xs) + +_event_maxestimate(xs::AbstractVector{<:Real}) = + isempty(xs) ? MeanEstimate(0.0) : MeanEstimate(maximum(xs)) mutable struct ShortfallEventsAccumulator{S} <: ResultAccumulator{ShortfallEvents} @@ -154,79 +159,83 @@ function finalize( end function MeanEventDuration(x::ShortfallEventsResult{N,L,T}) where {N,L,T} - durations = Float64[ - isempty(events) ? 0.0 : mean(duration_periods.(events)) + durations = [ + duration_periods(ev) for events in x.system_events + for ev in events ] - return MeanEventDuration{N,L,T}(MeanEstimate(durations)) + return MeanEventDuration{N,L,T}(_event_meanestimate(durations)) end function MeanEventDuration(x::ShortfallEventsResult{N,L,T}, r::AbstractString) where {N,L,T} i_r = findfirstunique(x.regions.names, r) - durations = Float64[ - isempty(x.region_events[i_r, s]) ? 0.0 : - mean(duration_periods.(x.region_events[i_r, s])) + durations = [ + duration_periods(ev) for s in axes(x.region_events, 2) + for ev in x.region_events[i_r, s] ] - return MeanEventDuration{N,L,T}(MeanEstimate(durations)) + return MeanEventDuration{N,L,T}(_event_meanestimate(durations)) end function MaxEventDuration(x::ShortfallEventsResult{N,L,T}) where {N,L,T} - durations = Float64[ - isempty(events) ? 0.0 : maximum(duration_periods.(events)) + durations = [ + duration_periods(ev) for events in x.system_events + for ev in events ] - return MaxEventDuration{N,L,T}(MeanEstimate(durations)) + return MaxEventDuration{N,L,T}(_event_maxestimate(durations)) end function MaxEventDuration(x::ShortfallEventsResult{N,L,T}, r::AbstractString) where {N,L,T} i_r = findfirstunique(x.regions.names, r) - durations = Float64[ - isempty(x.region_events[i_r, s]) ? 0.0 : - maximum(duration_periods.(x.region_events[i_r, s])) + durations = [ + duration_periods(ev) for s in axes(x.region_events, 2) + for ev in x.region_events[i_r, s] ] - return MaxEventDuration{N,L,T}(MeanEstimate(durations)) + return MaxEventDuration{N,L,T}(_event_maxestimate(durations)) end function MeanEventEnergy(x::ShortfallEventsResult{N,L,T,P,E}) where {N,L,T,P,E} p2e = conversionfactor(L, T, P, E) - energies = Float64[ - isempty(events) ? 0.0 : mean(p2e .* event_energy.(events)) + energies = [ + p2e * event_energy(ev) for events in x.system_events + for ev in events ] - return MeanEventEnergy{N,L,T,E}(MeanEstimate(energies)) + return MeanEventEnergy{N,L,T,E}(_event_meanestimate(energies)) end function MeanEventEnergy(x::ShortfallEventsResult{N,L,T,P,E}, r::AbstractString) where {N,L,T,P,E} i_r = findfirstunique(x.regions.names, r) p2e = conversionfactor(L, T, P, E) - energies = Float64[ - isempty(x.region_events[i_r, s]) ? 0.0 : - mean(p2e .* event_energy.(x.region_events[i_r, s])) + energies = [ + p2e * event_energy(ev) for s in axes(x.region_events, 2) + for ev in x.region_events[i_r, s] ] - return MeanEventEnergy{N,L,T,E}(MeanEstimate(energies)) + return MeanEventEnergy{N,L,T,E}(_event_meanestimate(energies)) end function MaxEventEnergy(x::ShortfallEventsResult{N,L,T,P,E}) where {N,L,T,P,E} p2e = conversionfactor(L, T, P, E) - energies = Float64[ - isempty(events) ? 0.0 : maximum(p2e .* event_energy.(events)) + energies = [ + p2e * event_energy(ev) for events in x.system_events + for ev in events ] - return MaxEventEnergy{N,L,T,E}(MeanEstimate(energies)) + return MaxEventEnergy{N,L,T,E}(_event_maxestimate(energies)) end function MaxEventEnergy(x::ShortfallEventsResult{N,L,T,P,E}, r::AbstractString) where {N,L,T,P,E} i_r = findfirstunique(x.regions.names, r) p2e = conversionfactor(L, T, P, E) - energies = Float64[ - isempty(x.region_events[i_r, s]) ? 0.0 : - maximum(p2e .* event_energy.(x.region_events[i_r, s])) + energies = [ + p2e * event_energy(ev) for s in axes(x.region_events, 2) + for ev in x.region_events[i_r, s] ] - return MaxEventEnergy{N,L,T,E}(MeanEstimate(energies)) + return MaxEventEnergy{N,L,T,E}(_event_maxestimate(energies)) end diff --git a/PRASCore.jl/test/Simulations/runtests.jl b/PRASCore.jl/test/Simulations/runtests.jl index 97691a00..b2a45608 100644 --- a/PRASCore.jl/test/Simulations/runtests.jl +++ b/PRASCore.jl/test/Simulations/runtests.jl @@ -758,30 +758,34 @@ manual_lolev_1a = mean(length.(events_1a.system_events)) @test isapprox(val(LOLEv(events_1a)), manual_lolev_1a; rtol=1e-10) - manual_meandur_1a = mean([ - isempty(evts) ? 0.0 : mean(Results.duration_periods.(evts)) + durations_1a = [ + Results.duration_periods(ev) for evts in events_1a.system_events - ]) + for ev in evts + ] + + manual_meandur_1a = isempty(durations_1a) ? 0.0 : mean(durations_1a) + @test isapprox(val(MeanEventDuration(events_1a)), manual_meandur_1a; rtol=1e-10) - manual_maxdur_1a = mean([ - isempty(evts) ? 0.0 : maximum(Results.duration_periods.(evts)) - for evts in events_1a.system_events - ]) + manual_maxdur_1a = isempty(durations_1a) ? 0.0 : maximum(durations_1a) + @test isapprox(val(MaxEventDuration(events_1a)), manual_maxdur_1a; rtol=1e-10) p2e_1a = PRASCore.Systems.conversionfactor(1, Hour, PRASCore.Systems.MW, PRASCore.Systems.MWh) - manual_meanenergy_1a = mean([ - isempty(evts) ? 0.0 : mean(p2e_1a .* Results.event_energy.(evts)) + energies_1a = [ + p2e_1a * Results.event_energy(ev) for evts in events_1a.system_events - ]) + for ev in evts + ] + + manual_meanenergy_1a = isempty(energies_1a) ? 0.0 : mean(energies_1a) + @test isapprox(val(MeanEventEnergy(events_1a)), manual_meanenergy_1a; rtol=1e-10) - manual_maxenergy_1a = mean([ - isempty(evts) ? 0.0 : maximum(p2e_1a .* Results.event_energy.(evts)) - for evts in events_1a.system_events - ]) + manual_maxenergy_1a = isempty(energies_1a) ? 0.0 : maximum(energies_1a) + @test isapprox(val(MaxEventEnergy(events_1a)), manual_maxenergy_1a; rtol=1e-10) # Multi-region system @@ -797,6 +801,48 @@ end + @testset "Event metrics return zero when no events exist" begin + sys = deepcopy(TestData.singlenode_a) + sys.regions.load .= 0 + + spec = SequentialMonteCarlo(samples=100, seed=42, threaded=false) + events, = assess(sys, spec, ShortfallEvents()) + + @test Results.totalevents(events) == 0 + + @test val(MeanEventDuration(events)) == 0.0 + @test val(MaxEventDuration(events)) == 0.0 + @test val(MeanEventEnergy(events)) == 0.0 + @test val(MaxEventEnergy(events)) == 0.0 + + @test stderror(MeanEventDuration(events)) == 0.0 + @test stderror(MaxEventDuration(events)) == 0.0 + @test stderror(MeanEventEnergy(events)) == 0.0 + @test stderror(MaxEventEnergy(events)) == 0.0 + end + + @testset "ShortfallEvents threaded and serial results match" begin + serial_spec = SequentialMonteCarlo(samples=100_000, seed=123, threaded=false) + threaded_spec = SequentialMonteCarlo(samples=100_000, seed=123, threaded=true) + + serial_events, = assess(TestData.threenode, serial_spec, ShortfallEvents()) + threaded_events, = assess(TestData.threenode, threaded_spec, ShortfallEvents()) + + @test LOLEv(serial_events) ≈ LOLEv(threaded_events) + @test MeanEventDuration(serial_events) ≈ MeanEventDuration(threaded_events) + @test MaxEventDuration(serial_events) ≈ MaxEventDuration(threaded_events) + @test MeanEventEnergy(serial_events) ≈ MeanEventEnergy(threaded_events) + @test MaxEventEnergy(serial_events) ≈ MaxEventEnergy(threaded_events) + + for r in serial_events.regions.names + @test LOLEv(serial_events, r) ≈ LOLEv(threaded_events, r) + @test MeanEventDuration(serial_events, r) ≈ MeanEventDuration(threaded_events, r) + @test MaxEventDuration(serial_events, r) ≈ MaxEventDuration(threaded_events, r) + @test MeanEventEnergy(serial_events, r) ≈ MeanEventEnergy(threaded_events, r) + @test MaxEventEnergy(serial_events, r) ≈ MaxEventEnergy(threaded_events, r) + end + end + @testset "Threaded sample result partitioning" begin simspec_serial = SequentialMonteCarlo(samples=100, seed=123, threaded=false) diff --git a/PRASFiles.jl/src/PRASFiles.jl b/PRASFiles.jl/src/PRASFiles.jl index 4d2a0420..bef1f819 100644 --- a/PRASFiles.jl/src/PRASFiles.jl +++ b/PRASFiles.jl/src/PRASFiles.jl @@ -2,7 +2,7 @@ module PRASFiles import PRASCore.Systems: SystemModel, Regions, Interfaces, Generators, Storages, GeneratorStorages, DemandResponses, Lines, - timeunits, powerunits, energyunits, unitsymbol + timeunits, powerunits, energyunits, unitsymbol, conversionfactor import PRASCore.Results: EUE, LOLE, NEUE, LOLD, @@ -10,7 +10,7 @@ import PRASCore.Results: AbstractShortfallResult, Result, ShortfallEventsResult, ShortfallEvent, LOLEv, totalevents, MeanEventDuration, MaxEventDuration, - MeanEventEnergy, MaxEventEnergy + MeanEventEnergy, MaxEventEnergy, findfirstunique import StatsBase: mean import Dates: @dateformat_str, format, now import TimeZones: ZonedDateTime diff --git a/PRASFiles.jl/src/Results/utils.jl b/PRASFiles.jl/src/Results/utils.jl index a7b8eaea..68fc2f19 100644 --- a/PRASFiles.jl/src/Results/utils.jl +++ b/PRASFiles.jl/src/Results/utils.jl @@ -259,7 +259,7 @@ function get_eventrecords( events::ShortfallEventsResult{N,L,T,P,E}, region::String, ) where {N,L,T,P,E} - i_r = findfirst(isequal(region), events.regions.names) + i_r = findfirstunique(events.regions.names, region) p2e = conversionfactor(L, T, P, E) records = EventRecord[] diff --git a/PRASFiles.jl/src/Results/write.jl b/PRASFiles.jl/src/Results/write.jl index f930cc97..9f6db1c2 100644 --- a/PRASFiles.jl/src/Results/write.jl +++ b/PRASFiles.jl/src/Results/write.jl @@ -162,8 +162,7 @@ end outfile::String, ) -Save `ShortfallEventsResult` in JSON format, including both summary event -metrics and raw event records. +Save `ShortfallEventsResult` summary metrics in JSON format, optionally raw event records. # Arguments From 14edc4b78c1c8a227869bd669afe9a33cc5c0b88 Mon Sep 17 00:00:00 2001 From: akrivi <20168326+akrivi@users.noreply.github.com> Date: Mon, 4 May 2026 08:57:57 -0600 Subject: [PATCH 05/13] fix tests --- PRASCore.jl/test/Simulations/runtests.jl | 22 ---------------------- 1 file changed, 22 deletions(-) diff --git a/PRASCore.jl/test/Simulations/runtests.jl b/PRASCore.jl/test/Simulations/runtests.jl index b2a45608..96e68a5e 100644 --- a/PRASCore.jl/test/Simulations/runtests.jl +++ b/PRASCore.jl/test/Simulations/runtests.jl @@ -821,28 +821,6 @@ @test stderror(MaxEventEnergy(events)) == 0.0 end - @testset "ShortfallEvents threaded and serial results match" begin - serial_spec = SequentialMonteCarlo(samples=100_000, seed=123, threaded=false) - threaded_spec = SequentialMonteCarlo(samples=100_000, seed=123, threaded=true) - - serial_events, = assess(TestData.threenode, serial_spec, ShortfallEvents()) - threaded_events, = assess(TestData.threenode, threaded_spec, ShortfallEvents()) - - @test LOLEv(serial_events) ≈ LOLEv(threaded_events) - @test MeanEventDuration(serial_events) ≈ MeanEventDuration(threaded_events) - @test MaxEventDuration(serial_events) ≈ MaxEventDuration(threaded_events) - @test MeanEventEnergy(serial_events) ≈ MeanEventEnergy(threaded_events) - @test MaxEventEnergy(serial_events) ≈ MaxEventEnergy(threaded_events) - - for r in serial_events.regions.names - @test LOLEv(serial_events, r) ≈ LOLEv(threaded_events, r) - @test MeanEventDuration(serial_events, r) ≈ MeanEventDuration(threaded_events, r) - @test MaxEventDuration(serial_events, r) ≈ MaxEventDuration(threaded_events, r) - @test MeanEventEnergy(serial_events, r) ≈ MeanEventEnergy(threaded_events, r) - @test MaxEventEnergy(serial_events, r) ≈ MaxEventEnergy(threaded_events, r) - end - end - @testset "Threaded sample result partitioning" begin simspec_serial = SequentialMonteCarlo(samples=100, seed=123, threaded=false) From 412db6436243ba40cdd6cf4af773841857e48408 Mon Sep 17 00:00:00 2001 From: akrivi <20168326+akrivi@users.noreply.github.com> Date: Thu, 3 Sep 2026 19:41:40 -0600 Subject: [PATCH 06/13] ShortfallEvents docs --- PRASCore.jl/src/Results/ShortfallEvents.jl | 38 +++++++++++++++++++--- 1 file changed, 34 insertions(+), 4 deletions(-) diff --git a/PRASCore.jl/src/Results/ShortfallEvents.jl b/PRASCore.jl/src/Results/ShortfallEvents.jl index 59959391..d5755c0c 100644 --- a/PRASCore.jl/src/Results/ShortfallEvents.jl +++ b/PRASCore.jl/src/Results/ShortfallEvents.jl @@ -5,10 +5,40 @@ The `ShortfallEvents` result specification reports sample-level shortfall events, producing a `ShortfallEventsResult`. A shortfall event is a contiguous run of one or more simulation timesteps -with positive shortfall. - -This result can be used to inspect event start/end times and to compute -event-based reliability metrics such as [`LOLEv`](@ref). +with positive shortfall. A `ShortfallEventsResult` can be indexed by sample +number to retrieve system-wide events or by region name and sample number to +retrieve regional events. Each event records its starting timestep, ending +timestep and unserved energy. + +Example: + +```julia +events, = + assess(sys, SequentialMonteCarlo(samples=1000), ShortfallEvents()) + +# Events for the first sample +system_events = events[1] +regional_events = events["Region A", 1] + +# Each event has start_idx, end_idx and energy fields +first_system_event = first(system_events) +start_idx = first_system_event.start_idx +end_idx = first_system_event.end_idx +energy = first_system_event.energy + +# System-wide event metrics +lolev = LOLEv(events) +mean_duration = MeanEventDuration(events) +max_duration = MaxEventDuration(events) +mean_energy = MeanEventEnergy(events) +max_energy = MaxEventEnergy(events) + +# Regional event metrics +regional_lolev = LOLEv(events, "Region A") +``` + +This result stores every shortfall event for every sample and can require +significant memory for simulations with many samples or events. """ struct ShortfallEvents <: ResultSpec end From 6756140337bc0c16c0a0f4dbef75cdc5e6a23065 Mon Sep 17 00:00:00 2001 From: akrivi <20168326+akrivi@users.noreply.github.com> Date: Thu, 3 Sep 2026 21:04:20 -0600 Subject: [PATCH 07/13] cleanup --- PRASCore.jl/src/Results/Results.jl | 2 - PRASCore.jl/src/Results/ShortfallEvents.jl | 20 +++--- PRASCore.jl/src/Simulations/recording.jl | 8 +-- PRASFiles.jl/src/PRASFiles.jl | 4 +- PRASFiles.jl/src/Results/utils.jl | 54 +++++++--------- PRASFiles.jl/src/Results/write.jl | 8 +-- PRASFiles.jl/test/runtests.jl | 75 ++++++++++++++++++++++ 7 files changed, 118 insertions(+), 53 deletions(-) diff --git a/PRASCore.jl/src/Results/Results.jl b/PRASCore.jl/src/Results/Results.jl index b7ba96bc..738a7381 100644 --- a/PRASCore.jl/src/Results/Results.jl +++ b/PRASCore.jl/src/Results/Results.jl @@ -189,8 +189,6 @@ getindex(x::AbstractEnergyResult, name::String, ::Colon) = getindex(x::AbstractEnergyResult, ::Colon, ::Colon) = getindex.(x, names(x), permutedims(x.timestamps)) -abstract type AbstractShortfallEventResult{N,L,T} <: Result{N,L,T} end - include("StorageEnergy.jl") include("GeneratorStorageEnergy.jl") include("DemandResponseEnergy.jl") diff --git a/PRASCore.jl/src/Results/ShortfallEvents.jl b/PRASCore.jl/src/Results/ShortfallEvents.jl index d5755c0c..7fbc084c 100644 --- a/PRASCore.jl/src/Results/ShortfallEvents.jl +++ b/PRASCore.jl/src/Results/ShortfallEvents.jl @@ -65,7 +65,7 @@ _event_meanestimate(xs::AbstractVector{<:Real}) = _event_maxestimate(xs::AbstractVector{<:Real}) = isempty(xs) ? MeanEstimate(0.0) : MeanEstimate(maximum(xs)) -mutable struct ShortfallEventsAccumulator{S} <: ResultAccumulator{ShortfallEvents} +mutable struct ShortfallEventsAccumulator <: ResultAccumulator{ShortfallEvents} system_events::Vector{Vector{ShortfallEvent}} region_events::Matrix{Vector{ShortfallEvent}} @@ -82,8 +82,8 @@ mutable struct ShortfallEventsAccumulator{S} <: ResultAccumulator{ShortfallEvent end function accumulator( - sys::SystemModel{N}, nsamples::Int, ::S -) where {N,S<:ShortfallEvents} + sys::SystemModel{N}, nsamples::Int, ::ShortfallEvents +) where {N} nregions = length(sys.regions) @@ -98,7 +98,7 @@ function accumulator( region_event_start = zeros(Int, nregions) region_event_energy = zeros(Int, nregions) - return ShortfallEventsAccumulator{S}( + return ShortfallEventsAccumulator( system_events, region_events, in_system_event, system_event_start, system_event_energy, in_region_event, region_event_start, region_event_energy, @@ -123,11 +123,9 @@ function copy_sample_partition!( return end -accumulatortype(::S) where { - S<:ShortfallEvents - } = ShortfallEventsAccumulator{S} +accumulatortype(::ShortfallEvents) = ShortfallEventsAccumulator -struct ShortfallEventsResult{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit,S} <: AbstractShortfallEventResult{N,L,T} +struct ShortfallEventsResult{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} <: Result{N,L,T} regions::Regions timestamps::StepRange{ZonedDateTime,T} @@ -178,12 +176,12 @@ function LOLEv(x::ShortfallEventsResult{N,L,T}, r::AbstractString) where {N,L,T} end function finalize( - acc::ShortfallEventsAccumulator{S}, + acc::ShortfallEventsAccumulator, system::SystemModel{N,L,T,P,E}, -) where {N,L,T,P,E,S<:ShortfallEvents} +) where {N,L,T,P,E} - return ShortfallEventsResult{N,L,T,P,E,S}( + return ShortfallEventsResult{N,L,T,P,E}( system.regions, system.timestamps, acc.system_events, acc.region_events) end diff --git a/PRASCore.jl/src/Simulations/recording.jl b/PRASCore.jl/src/Simulations/recording.jl index 5429523c..2034b22d 100644 --- a/PRASCore.jl/src/Simulations/recording.jl +++ b/PRASCore.jl/src/Simulations/recording.jl @@ -517,11 +517,11 @@ reset!(acc::Results.DemandResponseEnergySamplesAccumulator, sampleid::Int) = not # ShortfallEvents function record!( - acc::Results.ShortfallEventsAccumulator{S}, + acc::Results.ShortfallEventsAccumulator, system::SystemModel{N,L,T,P,E}, state::SystemState, problem::DispatchProblem, sampleid::Int, t::Int -) where {N,L,T,P,E,S} +) where {N,L,T,P,E} isshortfall = false totalshortfall = 0 @@ -529,7 +529,7 @@ function record!( for (r, dr_idxs) in zip(problem.region_unserved_edges, system.region_dr_idxs) - regionshortfall = init_regionshortfall(S, edges, r) + regionshortfall = init_regionshortfall(Results.ShortfallEvents, edges, r) dr_shortfall = 0 for i in dr_idxs @@ -613,4 +613,4 @@ function reset!(acc::Results.ShortfallEventsAccumulator, sampleid::Int) end return -end \ No newline at end of file +end diff --git a/PRASFiles.jl/src/PRASFiles.jl b/PRASFiles.jl/src/PRASFiles.jl index bef1f819..58c0a56e 100644 --- a/PRASFiles.jl/src/PRASFiles.jl +++ b/PRASFiles.jl/src/PRASFiles.jl @@ -10,7 +10,9 @@ import PRASCore.Results: AbstractShortfallResult, Result, ShortfallEventsResult, ShortfallEvent, LOLEv, totalevents, MeanEventDuration, MaxEventDuration, - MeanEventEnergy, MaxEventEnergy, findfirstunique + MeanEventEnergy, MaxEventEnergy, findfirstunique, + duration_periods, event_energy, + start_event_timestamp, end_event_timestamp import StatsBase: mean import Dates: @dateformat_str, format, now import TimeZones: ZonedDateTime diff --git a/PRASFiles.jl/src/Results/utils.jl b/PRASFiles.jl/src/Results/utils.jl index 68fc2f19..08d5c20b 100644 --- a/PRASFiles.jl/src/Results/utils.jl +++ b/PRASFiles.jl/src/Results/utils.jl @@ -233,49 +233,43 @@ function get_lold_result( return LOLDResult(shortfall; region = region) end -function get_eventrecords( - events_by_sample::Vector{Vector{ShortfallEvent}}, - timestamps, - p2e, -) - records = EventRecord[] +function _get_eventrecords( + events::ShortfallEventsResult{N,L,T,P,E}, + events_by_sample::AbstractVector{<:AbstractVector{ShortfallEvent}}, +) where {N,L,T,P,E} + p2e = conversionfactor(L, T, P, E) + nrecords = sum(length, events_by_sample) + iszero(nrecords) && return EventRecord[] + + records = Vector{EventRecord}(undef, nrecords) + record_idx = 1 for (sample_id, evts) in enumerate(events_by_sample) for ev in evts - push!(records, EventRecord( + records[record_idx] = EventRecord( sample_id, - timestamps[ev.start_idx], - timestamps[ev.end_idx], - ev.end_idx - ev.start_idx + 1, - p2e * ev.energy, - )) + start_event_timestamp(events, ev), + end_event_timestamp(events, ev), + duration_periods(ev), + p2e * event_energy(ev), + ) + record_idx += 1 end end return records end +function get_eventrecords(events::ShortfallEventsResult) + return _get_eventrecords(events, events.system_events) +end + function get_eventrecords( - events::ShortfallEventsResult{N,L,T,P,E}, + events::ShortfallEventsResult, region::String, -) where {N,L,T,P,E} +) i_r = findfirstunique(events.regions.names, region) - p2e = conversionfactor(L, T, P, E) - - records = EventRecord[] - for sample_id in axes(events.region_events, 2) - for ev in events.region_events[i_r, sample_id] - push!(records, EventRecord( - sample_id, - events.timestamps[ev.start_idx], - events.timestamps[ev.end_idx], - ev.end_idx - ev.start_idx + 1, - p2e * ev.energy, - )) - end - end - - return records + return _get_eventrecords(events, view(events.region_events, i_r, :)) end # Define structtypes for different structs defined above diff --git a/PRASFiles.jl/src/Results/write.jl b/PRASFiles.jl/src/Results/write.jl index 9f6db1c2..39030519 100644 --- a/PRASFiles.jl/src/Results/write.jl +++ b/PRASFiles.jl/src/Results/write.jl @@ -106,15 +106,13 @@ function saveshortfall( end function generate_eventresult( - events::ShortfallEventsResult{N,L,T,P,E}, + events::ShortfallEventsResult, pras_sys::SystemModel; include_events::Bool = false, -) where {N,L,T,P,E} - - p2e = conversionfactor(L, T, P, E) +) system_event_records = include_events ? - get_eventrecords(events.system_events, events.timestamps, p2e) : + get_eventrecords(events) : EventRecord[] region_results = RegionEventResult[] diff --git a/PRASFiles.jl/test/runtests.jl b/PRASFiles.jl/test/runtests.jl index eb6fdc90..a9aa7275 100644 --- a/PRASFiles.jl/test/runtests.jl +++ b/PRASFiles.jl/test/runtests.jl @@ -159,4 +159,79 @@ using JSON3 @test_throws "saveshortfall is not implemented for" PRASFiles.saveshortfall(surplus, rts_sys, path) end + @testset "Save Event Results" begin + rts_sys = PRASFiles.rts_gmlc() + # Make load in all regions 10 times the original load for meaningful results + for i in 1:length(rts_sys.regions.names) + rts_sys.regions.load[i, :] = 10 * rts_sys.regions.load[i, :] + end + + events, = assess( + rts_sys, + SequentialMonteCarlo(samples = 10, threaded = false, seed = 1), + ShortfallEvents(), + ) + + summary_result = PRASFiles.generate_eventresult(events, rts_sys) + event_result = PRASFiles.generate_eventresult( + events, + rts_sys; + include_events = true, + ) + + @test summary_result.total_events == event_result.total_events + @test isempty(summary_result.system_events) + @test all(isempty(region.events) for region in summary_result.region_results) + @test length(event_result.system_events) == event_result.total_events + @test all( + length(region.events) == region.total_events + for region in event_result.region_results + ) + @test event_result.lolev.mean == PRASCore.LOLEv(events).lolev.estimate + @test event_result.mean_event_duration.mean == + PRASCore.MeanEventDuration(events).duration.estimate + @test event_result.mean_event_energy.mean == + PRASCore.MeanEventEnergy(events).energy.estimate + + @test !isempty(event_result.system_events) + first_record = first(event_result.system_events) + first_event = first(events.system_events[first_record.sample_id]) + @test first_record.start_timestamp == events.timestamps[first_event.start_idx] + @test first_record.end_timestamp == events.timestamps[first_event.end_idx] + @test first_record.duration_periods == + PRASCore.Results.duration_periods(first_event) + (_, period_length, period_unit, power_unit, energy_unit) = get_params(rts_sys) + @test first_record.energy == + conversionfactor(period_length, period_unit, power_unit, energy_unit) * + PRASCore.Results.event_energy(first_event) + + path = joinpath(dirname(@__FILE__), "PRAS_Results_Export") + exp_location = PRASFiles.saveevents( + events, + rts_sys, + path; + include_events = true, + ) + result_path = joinpath(exp_location, "pras_event_results.json") + @test isfile(result_path) + exp_result = JSON3.read(result_path, PRASFiles.SystemEventResult) + @test exp_result.total_events == event_result.total_events + @test length(exp_result.system_events) == event_result.total_events + + empty_sys = PRASFiles.rts_gmlc() + fill!(empty_sys.regions.load, 0) + empty_events, = assess( + empty_sys, + SequentialMonteCarlo(samples = 2, threaded = false, seed = 1), + ShortfallEvents(), + ) + empty_result = PRASFiles.generate_eventresult( + empty_events, + empty_sys; + include_events = true, + ) + @test isempty(empty_result.system_events) + @test all(isempty(region.events) for region in empty_result.region_results) + end + end From 8342650297f2cd1458ac6eb4131ed1ddba94518e Mon Sep 17 00:00:00 2001 From: akrivi <20168326+akrivi@users.noreply.github.com> Date: Sat, 5 Sep 2026 17:41:50 -0600 Subject: [PATCH 08/13] test update --- PRASCore.jl/test/Simulations/runtests.jl | 26 ++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/PRASCore.jl/test/Simulations/runtests.jl b/PRASCore.jl/test/Simulations/runtests.jl index 96e68a5e..ae717f38 100644 --- a/PRASCore.jl/test/Simulations/runtests.jl +++ b/PRASCore.jl/test/Simulations/runtests.jl @@ -794,13 +794,35 @@ @test val(LOLEv(events_3, "Region A")) >= 0 @test val(MeanEventDuration(events_3, "Region A")) >= 0 - @test val(LOLEv(events_3)) >= val(LOLEv(events_3, "Region A")) - @test Results.totalevents(events_1a) >= 0 @test Results.totalevents(events_3, "Region A") >= 0 end + @testset "System events remain continuous across regions" begin + sys = deepcopy(TestData.threenode) + sys.generators.capacity .= 0 + sys.interfaces.limit_forward .= 0 + sys.interfaces.limit_backward .= 0 + sys.regions.load .= [1 0 1 0; 0 1 0 0; 0 0 0 0] + + spec = SequentialMonteCarlo(samples=2, seed=42, threaded=false) + events, = assess(sys, spec, ShortfallEvents()) + + # Region B bridges the gap between Region A's two events. + for s in 1:spec.nsamples + @test [(ev.start_idx, ev.end_idx) for ev in events[s]] == [(1, 3)] + @test [(ev.start_idx, ev.end_idx) for ev in events["Region A", s]] == + [(1, 1), (3, 3)] + @test [(ev.start_idx, ev.end_idx) for ev in events["Region B", s]] == + [(2, 2)] + end + + @test val(LOLEv(events)) == 1.0 + @test val(LOLEv(events, "Region A")) == 2.0 + @test val(LOLEv(events, "Region B")) == 1.0 + end + @testset "Event metrics return zero when no events exist" begin sys = deepcopy(TestData.singlenode_a) sys.regions.load .= 0 From c46ba2e9a61fea67b5c786672350ac43f31ad414 Mon Sep 17 00:00:00 2001 From: akrivi <20168326+akrivi@users.noreply.github.com> Date: Sat, 5 Sep 2026 17:59:29 -0600 Subject: [PATCH 09/13] simplify max calculation --- PRASCore.jl/src/Results/ShortfallEvents.jl | 27 ++++++++++------------ PRASFiles.jl/src/Results/utils.jl | 4 ++++ PRASFiles.jl/src/Results/write.jl | 2 +- 3 files changed, 17 insertions(+), 16 deletions(-) diff --git a/PRASCore.jl/src/Results/ShortfallEvents.jl b/PRASCore.jl/src/Results/ShortfallEvents.jl index 7fbc084c..1587b8bd 100644 --- a/PRASCore.jl/src/Results/ShortfallEvents.jl +++ b/PRASCore.jl/src/Results/ShortfallEvents.jl @@ -62,9 +62,6 @@ event_energy(ev::ShortfallEvent) = ev.energy _event_meanestimate(xs::AbstractVector{<:Real}) = isempty(xs) ? MeanEstimate(0.0) : MeanEstimate(xs) -_event_maxestimate(xs::AbstractVector{<:Real}) = - isempty(xs) ? MeanEstimate(0.0) : MeanEstimate(maximum(xs)) - mutable struct ShortfallEventsAccumulator <: ResultAccumulator{ShortfallEvents} system_events::Vector{Vector{ShortfallEvent}} @@ -206,22 +203,22 @@ function MeanEventDuration(x::ShortfallEventsResult{N,L,T}, r::AbstractString) w end function MaxEventDuration(x::ShortfallEventsResult{N,L,T}) where {N,L,T} - durations = [ + duration = maximum(( duration_periods(ev) for events in x.system_events for ev in events - ] - return MaxEventDuration{N,L,T}(_event_maxestimate(durations)) + ); init=0) + return MaxEventDuration{N,L,T}(MeanEstimate(duration)) end function MaxEventDuration(x::ShortfallEventsResult{N,L,T}, r::AbstractString) where {N,L,T} i_r = findfirstunique(x.regions.names, r) - durations = [ + duration = maximum(( duration_periods(ev) for s in axes(x.region_events, 2) for ev in x.region_events[i_r, s] - ] - return MaxEventDuration{N,L,T}(_event_maxestimate(durations)) + ); init=0) + return MaxEventDuration{N,L,T}(MeanEstimate(duration)) end function MeanEventEnergy(x::ShortfallEventsResult{N,L,T,P,E}) where {N,L,T,P,E} @@ -247,23 +244,23 @@ end function MaxEventEnergy(x::ShortfallEventsResult{N,L,T,P,E}) where {N,L,T,P,E} p2e = conversionfactor(L, T, P, E) - energies = [ + energy = maximum(( p2e * event_energy(ev) for events in x.system_events for ev in events - ] - return MaxEventEnergy{N,L,T,E}(_event_maxestimate(energies)) + ); init=0.0) + return MaxEventEnergy{N,L,T,E}(MeanEstimate(energy)) end function MaxEventEnergy(x::ShortfallEventsResult{N,L,T,P,E}, r::AbstractString) where {N,L,T,P,E} i_r = findfirstunique(x.regions.names, r) p2e = conversionfactor(L, T, P, E) - energies = [ + energy = maximum(( p2e * event_energy(ev) for s in axes(x.region_events, 2) for ev in x.region_events[i_r, s] - ] - return MaxEventEnergy{N,L,T,E}(_event_maxestimate(energies)) + ); init=0.0) + return MaxEventEnergy{N,L,T,E}(MeanEstimate(energy)) end diff --git a/PRASFiles.jl/src/Results/utils.jl b/PRASFiles.jl/src/Results/utils.jl index 08d5c20b..6ffab949 100644 --- a/PRASFiles.jl/src/Results/utils.jl +++ b/PRASFiles.jl/src/Results/utils.jl @@ -212,6 +212,10 @@ function get_nsamples(shortfall::ShortfallSamplesResult) return size(shortfall.shortfall,3) end +function get_nsamples(events::ShortfallEventsResult) + return length(events.system_events) +end + function get_lold_result( shortfall::ShortfallResult, log_lold_info::Bool; diff --git a/PRASFiles.jl/src/Results/write.jl b/PRASFiles.jl/src/Results/write.jl index 39030519..66e5e8f3 100644 --- a/PRASFiles.jl/src/Results/write.jl +++ b/PRASFiles.jl/src/Results/write.jl @@ -136,7 +136,7 @@ function generate_eventresult( end sys_result = SystemEventResult( - length(events.system_events), + get_nsamples(events), TypeParams(pras_sys), pras_sys.attrs, collect(events.timestamps), From fcc8dbb9acd884068e9b8e3aeaa3a92028293837 Mon Sep 17 00:00:00 2001 From: akrivi <20168326+akrivi@users.noreply.github.com> Date: Sun, 6 Sep 2026 14:19:29 -0600 Subject: [PATCH 10/13] fix subhourly json --- PRASFiles.jl/src/Results/utils.jl | 6 +- PRASFiles.jl/test/runtests.jl | 129 +++++++++++++++++++----------- 2 files changed, 86 insertions(+), 49 deletions(-) diff --git a/PRASFiles.jl/src/Results/utils.jl b/PRASFiles.jl/src/Results/utils.jl index 6ffab949..dda63520 100644 --- a/PRASFiles.jl/src/Results/utils.jl +++ b/PRASFiles.jl/src/Results/utils.jl @@ -200,8 +200,10 @@ function get_shortfall_mean(shortfall::ShortfallResult) return shortfall.shortfall_mean end -function get_shortfall_mean(shortfall::ShortfallSamplesResult) - return mean(shortfall.shortfall, dims = 3) +function get_shortfall_mean(shortfall::ShortfallSamplesResult{N,L,T,P,E}) where {N,L,T,P,E} + shortfall_mean = mean(shortfall.shortfall, dims = 3) + shortfall_mean .*= conversionfactor(L, T, P, E) + return shortfall_mean end function get_nsamples(shortfall::ShortfallResult) diff --git a/PRASFiles.jl/test/runtests.jl b/PRASFiles.jl/test/runtests.jl index a9aa7275..30bc29e1 100644 --- a/PRASFiles.jl/test/runtests.jl +++ b/PRASFiles.jl/test/runtests.jl @@ -110,53 +110,88 @@ using JSON3 for i in 1:length(rts_sys.regions.names) rts_sys.regions.load[i, :] = 10 * rts_sys.regions.load[i, :] end - results = assess(rts_sys, SequentialMonteCarlo(samples=10, threaded = false, seed = 1), Shortfall(), ShortfallSamples(), Surplus()); - shortfall = results[1]; - lold_message = r"LOLD is not implemented for ShortfallResult" - @test_logs (:info, lold_message) PRASFiles.generate_systemresult(shortfall, rts_sys) - @test_logs (:info, lold_message) PRASFiles.generate_systemresult(shortfall, rts_sys) - path = joinpath(dirname(@__FILE__),"PRAS_Results_Export"); - exp_location_1 = PRASFiles.saveshortfall(shortfall, rts_sys, path); - @test isfile(joinpath(exp_location_1, "pras_results.json")) - exp_results_1 = JSON3.read(joinpath(exp_location_1, "pras_results.json"), PRASFiles.SystemResult) - @test exp_results_1.lole.mean == PRASCore.LOLE(shortfall).lole.estimate - @test exp_results_1.eue.mean == PRASCore.EUE(shortfall).eue.estimate - @test exp_results_1.neue.mean == PRASCore.NEUE(shortfall).neue.estimate - @test exp_results_1.region_results[1].lole.mean == PRASCore.LOLE(shortfall, exp_results_1.region_results[1].name).lole.estimate - @test exp_results_1.region_results[1].eue.mean == PRASCore.EUE(shortfall, exp_results_1.region_results[1].name).eue.estimate - @test exp_results_1.region_results[1].neue.mean == PRASCore.NEUE(shortfall, exp_results_1.region_results[1].name).neue.estimate - @test exp_results_1.lold === nothing - @test exp_results_1.region_results[1].lold === nothing - - shortfall_samples = results[2]; - exp_location_2 = PRASFiles.saveshortfall(shortfall_samples, rts_sys, path); - @test isfile(joinpath(exp_location_2, "pras_results.json")) - exp_results_2 = JSON3.read(joinpath(exp_location_2, "pras_results.json"), PRASFiles.SystemResult) - @test exp_results_2.lole.mean == PRASCore.LOLE(shortfall_samples).lole.estimate - @test exp_results_2.eue.mean == PRASCore.EUE(shortfall_samples).eue.estimate - @test exp_results_2.neue.mean == PRASCore.NEUE(shortfall_samples).neue.estimate - @test exp_results_2.region_results[1].lole.mean == PRASCore.LOLE(shortfall_samples, exp_results_2.region_results[1].name).lole.estimate - @test exp_results_2.region_results[1].eue.mean == PRASCore.EUE(shortfall_samples, exp_results_2.region_results[1].name).eue.estimate - @test exp_results_2.region_results[1].neue.mean == PRASCore.NEUE(shortfall_samples, exp_results_2.region_results[1].name).neue.estimate - @test exp_results_2.lold.mean == PRASCore.LOLD(shortfall_samples).lold.estimate - @test exp_results_2.lold.stderror == PRASCore.LOLD(shortfall_samples).lold.standarderror - region_name = exp_results_2.region_results[1].name - @test exp_results_2.region_results[1].lold.mean == PRASCore.LOLD(shortfall_samples, region_name).lold.estimate - @test exp_results_2.region_results[1].lold.stderror == PRASCore.LOLD(shortfall_samples, region_name).lold.standarderror - - @test exp_results_1.lole.mean ≈ exp_results_2.lole.mean - @test exp_results_1.eue.mean ≈ exp_results_2.eue.mean - @test exp_results_1.neue.mean ≈ exp_results_2.neue.mean - @test exp_results_1.region_results[1].lole.mean ≈ exp_results_2.region_results[1].lole.mean - @test exp_results_1.region_results[1].eue.mean ≈ exp_results_2.region_results[1].eue.mean - @test exp_results_1.region_results[1].neue.mean ≈ exp_results_2.region_results[1].neue.mean - - @test exp_results_1.lold === nothing - @test exp_results_2.lold !== nothing - @test exp_results_2.region_results[1].lold !== nothing - - surplus = results[3] - @test_throws "saveshortfall is not implemented for" PRASFiles.saveshortfall(surplus, rts_sys, path) + export_cases = ( + ("Hourly", rts_sys, 10), + ("Subhourly", PRASCore.Systems.TestData.singlenode_a_5min, 100), + ) + + for (name, system, nsamples) in export_cases + @testset "$name" begin + results = assess(system, SequentialMonteCarlo(samples=nsamples, threaded = false, seed = 1), Shortfall(), ShortfallSamples(), Surplus()); + shortfall = results[1]; + lold_message = r"LOLD is not implemented for ShortfallResult" + @test_logs (:info, lold_message) PRASFiles.generate_systemresult(shortfall, system) + @test_logs (:info, lold_message) PRASFiles.generate_systemresult(shortfall, system) + path = joinpath(dirname(@__FILE__), "PRAS_Results_Export", name); + exp_location_1 = PRASFiles.saveshortfall(shortfall, system, joinpath(path, "shortfall")); + @test isfile(joinpath(exp_location_1, "pras_results.json")) + exp_results_1 = JSON3.read(joinpath(exp_location_1, "pras_results.json"), PRASFiles.SystemResult) + @test exp_results_1.lole.mean == PRASCore.LOLE(shortfall).lole.estimate + @test exp_results_1.eue.mean == PRASCore.EUE(shortfall).eue.estimate + @test exp_results_1.neue.mean == PRASCore.NEUE(shortfall).neue.estimate + @test exp_results_1.region_results[1].lole.mean == PRASCore.LOLE(shortfall, exp_results_1.region_results[1].name).lole.estimate + @test exp_results_1.region_results[1].eue.mean == PRASCore.EUE(shortfall, exp_results_1.region_results[1].name).eue.estimate + @test exp_results_1.region_results[1].neue.mean == PRASCore.NEUE(shortfall, exp_results_1.region_results[1].name).neue.estimate + @test exp_results_1.lold === nothing + @test exp_results_1.region_results[1].lold === nothing + + shortfall_samples = results[2]; + exp_location_2 = PRASFiles.saveshortfall(shortfall_samples, system, joinpath(path, "samples")); + @test isfile(joinpath(exp_location_2, "pras_results.json")) + exp_results_2 = JSON3.read(joinpath(exp_location_2, "pras_results.json"), PRASFiles.SystemResult) + @test exp_results_2.lole.mean == PRASCore.LOLE(shortfall_samples).lole.estimate + @test exp_results_2.eue.mean == PRASCore.EUE(shortfall_samples).eue.estimate + @test exp_results_2.neue.mean == PRASCore.NEUE(shortfall_samples).neue.estimate + @test exp_results_2.region_results[1].lole.mean == PRASCore.LOLE(shortfall_samples, exp_results_2.region_results[1].name).lole.estimate + @test exp_results_2.region_results[1].eue.mean == PRASCore.EUE(shortfall_samples, exp_results_2.region_results[1].name).eue.estimate + @test exp_results_2.region_results[1].neue.mean == PRASCore.NEUE(shortfall_samples, exp_results_2.region_results[1].name).neue.estimate + @test exp_results_2.lold.mean == PRASCore.LOLD(shortfall_samples).lold.estimate + @test exp_results_2.lold.stderror == PRASCore.LOLD(shortfall_samples).lold.standarderror + region_name = exp_results_2.region_results[1].name + @test exp_results_2.region_results[1].lold.mean == PRASCore.LOLD(shortfall_samples, region_name).lold.estimate + @test exp_results_2.region_results[1].lold.stderror == PRASCore.LOLD(shortfall_samples, region_name).lold.standarderror + + @test any(>(0), exp_results_1.region_results[1].shortfall_mean) + @test exp_results_1.num_samples == exp_results_2.num_samples + @test exp_results_1.type_params.N == exp_results_2.type_params.N + @test exp_results_1.type_params.L == exp_results_2.type_params.L + @test exp_results_1.type_params.T == exp_results_2.type_params.T + @test exp_results_1.type_params.P == exp_results_2.type_params.P + @test exp_results_1.type_params.E == exp_results_2.type_params.E + @test exp_results_1.sys_attributes == exp_results_2.sys_attributes + @test exp_results_1.timestamps == exp_results_2.timestamps + @test exp_results_1.lole.mean ≈ exp_results_2.lole.mean + @test exp_results_1.lole.stderror ≈ exp_results_2.lole.stderror + @test exp_results_1.eue.mean ≈ exp_results_2.eue.mean + @test exp_results_1.eue.stderror ≈ exp_results_2.eue.stderror + @test exp_results_1.neue.mean ≈ exp_results_2.neue.mean + @test exp_results_1.neue.stderror ≈ exp_results_2.neue.stderror + + @test length(exp_results_1.region_results) == length(exp_results_2.region_results) + for (region_1, region_2) in zip(exp_results_1.region_results, exp_results_2.region_results) + @test region_1.name == region_2.name + @test region_1.load == region_2.load + @test region_1.peak_load == region_2.peak_load + @test region_1.capacity == region_2.capacity + @test region_1.shortfall_timestamps == region_2.shortfall_timestamps + @test region_1.shortfall_mean ≈ region_2.shortfall_mean + @test region_1.lole.mean ≈ region_2.lole.mean + @test region_1.lole.stderror ≈ region_2.lole.stderror + @test region_1.eue.mean ≈ region_2.eue.mean + @test region_1.eue.stderror ≈ region_2.eue.stderror + @test region_1.neue.mean ≈ region_2.neue.mean + @test region_1.neue.stderror ≈ region_2.neue.stderror + @test region_1.lold === nothing + @test region_2.lold !== nothing + end + + @test exp_results_1.lold === nothing + @test exp_results_2.lold !== nothing + + surplus = results[3] + @test_throws "saveshortfall is not implemented for" PRASFiles.saveshortfall(surplus, system, path) + end + end end @testset "Save Event Results" begin From abc31dfc026fcba353fdbe5f961c7e91064db6a5 Mon Sep 17 00:00:00 2001 From: akrivi <20168326+akrivi@users.noreply.github.com> Date: Wed, 9 Sep 2026 11:49:57 -0600 Subject: [PATCH 11/13] eventsinterval implementation --- PRASCore.jl/src/Results/Results.jl | 5 +- PRASCore.jl/src/Results/ShortfallEvents.jl | 53 +++++++++++++++ PRASCore.jl/test/Results/shortfall.jl | 75 ++++++++++++++++++++++ docs/src/PRASCore/api.md | 1 + 4 files changed, 132 insertions(+), 2 deletions(-) diff --git a/PRASCore.jl/src/Results/Results.jl b/PRASCore.jl/src/Results/Results.jl index 738a7381..01347aa7 100644 --- a/PRASCore.jl/src/Results/Results.jl +++ b/PRASCore.jl/src/Results/Results.jl @@ -17,7 +17,7 @@ export MeanEventDuration, MaxEventDuration, MeanEventEnergy, MaxEventEnergy, val, stderror, CVAR, NCVAR, - # Result specifications + # Result specifications and accessors Shortfall, ShortfallSamples, DemandResponseShortfall, DemandResponseShortfallSamples, Surplus, SurplusSamples, @@ -27,7 +27,8 @@ export DemandResponseEnergy, DemandResponseEnergySamples, GeneratorAvailability, StorageAvailability, GeneratorStorageAvailability,DemandResponseAvailability, - LineAvailability, ShortfallEvents + LineAvailability, ShortfallEvents, + eventsinterval include("metrics.jl") include("utils.jl") diff --git a/PRASCore.jl/src/Results/ShortfallEvents.jl b/PRASCore.jl/src/Results/ShortfallEvents.jl index 1587b8bd..4bc1629d 100644 --- a/PRASCore.jl/src/Results/ShortfallEvents.jl +++ b/PRASCore.jl/src/Results/ShortfallEvents.jl @@ -163,6 +163,59 @@ end start_event_timestamp(x::ShortfallEventsResult, ev::ShortfallEvent) = x.timestamps[ev.start_idx] end_event_timestamp(x::ShortfallEventsResult, ev::ShortfallEvent) = x.timestamps[ev.end_idx] +""" + eventsinterval(x::ShortfallEventsResult, t::StepRange{ZonedDateTime}) + eventsinterval(x::ShortfallEventsResult, r::AbstractString, t::StepRange{ZonedDateTime}) + +Return system-wide events, or events in region `r`, fully contained between the first and last timestamps of `t`, inclusive. +The interval must be nonempty and forward in time with endpoints present in `x.timestamps`. +All simulation timesteps between the endpoints are considered, regardless of the step of `t`. + +The returned vector contains a new event list for each original sample, including empty lists. +Event boundaries and stored energy are unchanged. +An `@info` message reports the number of overlapping events excluded because they start before or end after the interval, including events spanning the entire interval. +An empty list therefore means no fully contained events, not necessarily no shortfall. + +For example, `eventsinterval(x, x.timestamps[2:7])` excludes events at timesteps 1–3 and 6–8 but includes an event at timesteps 3–6. +""" +function eventsinterval(x::ShortfallEventsResult, t::StepRange{ZonedDateTime}) + return _eventsinterval(x, t, x.system_events) +end + +function eventsinterval( + x::ShortfallEventsResult, r::AbstractString, t::StepRange{ZonedDateTime} +) + i_r = findfirstunique(x.regions.names, r) + return _eventsinterval(x, t, view(x.region_events, i_r, :)) +end + +function _eventsinterval( + x::ShortfallEventsResult, t::StepRange{ZonedDateTime}, + sample_events::AbstractVector{Vector{ShortfallEvent}} +) + (isempty(t) || first(t) > last(t)) && + throw(ArgumentError("The event interval must be nonempty and forward in time")) + i_t0 = findfirstunique(x.timestamps, first(t)) + i_tf = findlastunique(x.timestamps, last(t)) + + selected = [ShortfallEvent[] for _ in eachindex(sample_events)] + nexcluded = 0 + for (included, events) in zip(selected, sample_events) + for ev in events + if i_t0 <= ev.start_idx && ev.end_idx <= i_tf + push!(included, ev) + elseif ev.start_idx <= i_tf && i_t0 <= ev.end_idx + nexcluded += 1 + end + end + end + + if nexcluded > 0 + @info "Excluded $nexcluded overlapping events that extend beyond the requested interval. Only fully contained events are returned." + end + return selected +end + LOLEv(x::ShortfallEventsResult{N,L,T}) where {N,L,T} = LOLEv{N,L,T}(MeanEstimate(length.(x.system_events))) diff --git a/PRASCore.jl/test/Results/shortfall.jl b/PRASCore.jl/test/Results/shortfall.jl index 503d6b72..27c72b22 100644 --- a/PRASCore.jl/test/Results/shortfall.jl +++ b/PRASCore.jl/test/Results/shortfall.jl @@ -261,3 +261,78 @@ end @test_throws BoundsError CVAR(:energy, result, alpha, r_bad, t_bad) end + +@testset "ShortfallEventsResult" begin + N = DD.nperiods + r, r_bad = DD.testresource, DD.notaresource + event = PRASCore.Results.ShortfallEvent + system_events = [ + [event(1, 3, 6), event(6, 8, 12)], + [event(3, 6, 20)], + event[], + ] + region_events = [event[] for _ in 1:DD.nresources, _ in 1:3] + region_events[1, 1] = copy(system_events[1]) + region_events[1, 2] = [event(3, 4, 6)] + region_events[2, 2] = [event(5, 6, 14)] + result = PRASCore.Results.ShortfallEventsResult{N,1,Hour,MW,MWh}( + Regions{N,MW}(DD.resourcenames, DD.resource_vals), + DD.periods, system_events, region_events) + + # System events and sample identity + + selected = @test_logs eventsinterval(result, DD.periods) + @test selected == system_events + @test length(selected) == 3 + @test all(selected[s] !== system_events[s] for s in eachindex(selected)) + empty!(selected[1]) + @test result[1] == [event(1, 3, 6), event(6, 8, 12)] + + excluded_two = r"Excluded 2 overlapping events" + selected = @test_logs (:info, excluded_two) eventsinterval(result, DD.periods[2:7]) + @test selected == [event[], [event(3, 6, 20)], event[]] + + # Inclusive boundaries and an event spanning the entire interval + + selected = @test_logs (:info, excluded_two) eventsinterval(result, DD.periods[3:6]) + @test selected == [event[], [event(3, 6, 20)], event[]] + selected = @test_logs (:info, r"Excluded 1 overlapping events") eventsinterval( + result, DD.periods[4:5]) + @test selected == [event[], event[], event[]] + selected = @test_logs (:info, excluded_two) eventsinterval(result, DD.periods[3:3]) + @test selected == [event[], event[], event[]] + + # Events entirely outside the interval do not generate a message. + + selected = @test_logs eventsinterval(result, DD.periods[9:10]) + @test selected == [event[], event[], event[]] + # The range specifies interval endpoints, not a sparse timestep selection. + selected = @test_logs (:info, excluded_two) eventsinterval(result, DD.periods[3:3:6]) + @test selected == [event[], [event(3, 6, 20)], event[]] + + # Region-specific selection uses only that region's event boundaries. + + selected = @test_logs eventsinterval(result, r, DD.periods[3:6]) + @test selected == [event[], [event(5, 6, 14)], event[]] + selected = @test_logs (:info, excluded_two) eventsinterval( + result, DD.resourcenames[1], DD.periods[3:6]) + @test selected == [event[], [event(3, 4, 6)], event[]] + selected = @test_logs eventsinterval(result, DD.resourcenames[3], DD.periods) + @test selected == [event[], event[], event[]] + + # Invalid selections follow the existing timestamp and region lookup errors. + + @test_throws BoundsError eventsinterval(result, r_bad, DD.periods) + @test_throws BoundsError eventsinterval(result, DD.badperiods) + @test_throws BoundsError eventsinterval(result, r, DD.badperiods) + @test_throws BoundsError eventsinterval( + result, first(DD.periods) - Hour(1):Hour(1):DD.periods[3]) + @test_throws BoundsError eventsinterval( + result, DD.periods[3]:Hour(1):last(DD.periods) + Hour(1)) + @test_throws InexactError eventsinterval( + result, DD.periods[3] + Minute(30):Hour(1):DD.periods[6] + Minute(30)) + @test_throws ArgumentError eventsinterval( + result, DD.periods[2]:Hour(1):DD.periods[1]) + @test_throws ArgumentError eventsinterval( + result, DD.periods[8]:Hour(-1):DD.periods[3]) +end diff --git a/docs/src/PRASCore/api.md b/docs/src/PRASCore/api.md index bc4f6908..3d68fb05 100644 --- a/docs/src/PRASCore/api.md +++ b/docs/src/PRASCore/api.md @@ -48,6 +48,7 @@ PRASCore.Results.NCVAR PRASCore.Results.val PRASCore.Results.stderror PRASCore.Results.ShortfallEvents +PRASCore.Results.eventsinterval PRASCore.Results.LOLEv PRASCore.Results.MeanEventDuration PRASCore.Results.MaxEventDuration From d628083686086c94ec166e2dd89a6ea9bb56e849 Mon Sep 17 00:00:00 2001 From: akrivi <20168326+akrivi@users.noreply.github.com> Date: Wed, 9 Sep 2026 12:14:55 -0600 Subject: [PATCH 12/13] simplify tests --- PRASCore.jl/test/Results/shortfall.jl | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/PRASCore.jl/test/Results/shortfall.jl b/PRASCore.jl/test/Results/shortfall.jl index 27c72b22..af10fce7 100644 --- a/PRASCore.jl/test/Results/shortfall.jl +++ b/PRASCore.jl/test/Results/shortfall.jl @@ -283,17 +283,12 @@ end selected = @test_logs eventsinterval(result, DD.periods) @test selected == system_events - @test length(selected) == 3 - @test all(selected[s] !== system_events[s] for s in eachindex(selected)) empty!(selected[1]) @test result[1] == [event(1, 3, 6), event(6, 8, 12)] - excluded_two = r"Excluded 2 overlapping events" - selected = @test_logs (:info, excluded_two) eventsinterval(result, DD.periods[2:7]) - @test selected == [event[], [event(3, 6, 20)], event[]] - # Inclusive boundaries and an event spanning the entire interval + excluded_two = r"Excluded 2 overlapping events" selected = @test_logs (:info, excluded_two) eventsinterval(result, DD.periods[3:6]) @test selected == [event[], [event(3, 6, 20)], event[]] selected = @test_logs (:info, r"Excluded 1 overlapping events") eventsinterval( @@ -329,8 +324,6 @@ end result, first(DD.periods) - Hour(1):Hour(1):DD.periods[3]) @test_throws BoundsError eventsinterval( result, DD.periods[3]:Hour(1):last(DD.periods) + Hour(1)) - @test_throws InexactError eventsinterval( - result, DD.periods[3] + Minute(30):Hour(1):DD.periods[6] + Minute(30)) @test_throws ArgumentError eventsinterval( result, DD.periods[2]:Hour(1):DD.periods[1]) @test_throws ArgumentError eventsinterval( From ce21b061c4257b7af63bcf830ba0f1b0c6cdfc51 Mon Sep 17 00:00:00 2001 From: akrivi <20168326+akrivi@users.noreply.github.com> Date: Sun, 6 Sep 2026 15:41:59 -0600 Subject: [PATCH 13/13] add walkthrough for shortfallevents --- PRAS.jl/examples/pras_adequacy_metrics.jl | 58 +- docs/make.jl | 1 + docs/src/PRAS/results.md | 103 ++- docs/src/PRASFiles/index.md | 48 ++ docs/src/images/shortfall-events.svg | 795 ++++++++++++++++++++++ docs/src/resourceadequacy.md | 6 +- 6 files changed, 986 insertions(+), 25 deletions(-) create mode 100644 docs/src/images/shortfall-events.svg diff --git a/PRAS.jl/examples/pras_adequacy_metrics.jl b/PRAS.jl/examples/pras_adequacy_metrics.jl index 81210819..38957d8c 100644 --- a/PRAS.jl/examples/pras_adequacy_metrics.jl +++ b/PRAS.jl/examples/pras_adequacy_metrics.jl @@ -23,10 +23,11 @@ # - An **adequacy event** is a set of event-periods that are contiguous at the highest available temporal resolution. # # These distinctions are important because different metrics count different temporal quantities. -# LOLE and LOLD correspond to the first two concepts: +# LOLE, LOLD and LOLEv correspond to these three concepts: # # - **LOLE** is the expected number of event-periods # - **LOLD** is the expected number of event-days +# - **LOLEv** is the expected number of adequacy events # # These metrics are related, but they are not interchangeable. # @@ -36,7 +37,7 @@ # ## Shortfall Severity # -# LOLE and LOLD describe when shortfalls occur, but they do not describe their magnitude. +# LOLE, LOLD and LOLEv describe the temporal occurrence of shortfalls but not their magnitude. # EUE complements these metrics by measuring the expected total amount of unserved energy over the study horizon. # ## Why Multiple Metrics Matter @@ -49,18 +50,20 @@ # We can consider a simple example of two cases next, for which we assume that # every shortfall hour has the same amount of unserved energy. # -# **Case A**: One day with 10 hours of shortfall +# **Case A**: One day with 10 consecutive hours of shortfall # -# **Case B**: Ten days with 1 hour of shortfall each +# **Case B**: Ten nonconsecutive days with 1 hour of shortfall each # # | Metric | Case A | Case B | # |------|--------|--------| # | LOLE | 10 | 10 | # | EUE | same | same | # | LOLD | 1 | 10 | +# | LOLEv | 1 | 10 | # # As we can see in the table above, even though LOLE and EUE are identical in this case, # LOLD reveals that shortfall events are more dispersed in Case B. +# Here LOLEv also distinguishes one long event from ten separate events. # # Because event-periods may be distributed across many days, a system with the @@ -124,6 +127,22 @@ # \end{cases} # ``` +# ### LOLEv +# +# LOLEv counts the expected number of transitions from no system shortfall to positive system shortfall: +# +# ```math +# \mathrm{LOLEv} = +# \mathbb{E}\left[\sum_t +# \mathbf{1}\left( +# \sum_r S_{r,t,s} > 0 +# \;\land\; +# \sum_r S_{r,t-1,s} = 0 +# \right)\right] +# ``` +# +# Here ``S_{r,0,s}=0`` so shortfall at the first timestep of each sample starts an event. +# # ### EUE # # EUE measures expected total unserved energy across the Monte Carlo samples: @@ -133,6 +152,15 @@ # \mathbb{E}\left[\sum_t \sum_r S_{r,t,s}\,\Delta t\right] # ``` # +# ## Event Duration and Energy +# +# `MeanEventDuration` and `MeanEventEnergy` average across all events pooled from all samples. +# Each event receives equal weight, so a sample containing several events contributes several observations. +# Samples without events contribute zero to the event count used by LOLEv but do not contribute a zero-duration or zero-energy event to these means. +# `MaxEventDuration` and `MaxEventEnergy` return the largest observed duration and energy across the pooled events, not the mean of each sample's maximum. +# +# See [Shortfall Events](@ref shortfall_events) for units, standard errors and result accessors. +# # ## Analysis with PRAS # # We revisit the [RTS-GMLC](https://github.com/GridMod/RTS-GMLC) system with increased load to induce shortfall, @@ -142,20 +170,36 @@ using PRAS sys = PRAS.rts_gmlc() sys.regions.load .+= 700.0 -shortfall_samples, = assess( +# We request both result types to calculate the metrics from the same simulated samples. +# `ShortfallSamples()` provides the sample-level results needed for LOLD while `ShortfallEvents()` provides the event records needed for the event metrics. +shortfall_samples, events = assess( sys, SequentialMonteCarlo(samples=100, seed=1), ShortfallSamples(), + ShortfallEvents(), ) -# and we calculate the metrics we discussed above: +# We then calculate the metrics discussed above: system_lole = LOLE(shortfall_samples) system_lold = LOLD(shortfall_samples) system_eue = EUE(shortfall_samples) +system_lolev = LOLEv(events) println(system_lole) println(system_lold) println(system_eue) +println(system_lolev) + +# We also calculate system-level event duration and energy statistics: +system_mean_event_duration = MeanEventDuration(events) +system_max_event_duration = MaxEventDuration(events) +system_mean_event_energy = MeanEventEnergy(events) +system_max_event_energy = MaxEventEnergy(events) + +println(system_mean_event_duration) +println(system_max_event_duration) +println(system_mean_event_energy) +println(system_max_event_energy) # We can also evaluate upper-tail severity by selecting a CVAR confidence level: alpha = 0.95 @@ -169,6 +213,8 @@ println(system_cvar) # EUE summarizes the average total unserved energy, while CVAR (``\alpha = 0.95``) summarizes # unserved energy in outcomes beyond the 95th-percentile threshold. +# See [Exporting shortfall events](@ref exporting_shortfall_events) to save the event summaries and individual records as JSON. + # ## References # diff --git a/docs/make.jl b/docs/make.jl index 827ceb64..ecd5eb48 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -44,6 +44,7 @@ makedocs( "Capacity Credit Calculation" => "PRAS/capacitycredit.md", ], ".pras File Format" => "SystemModel_HDF5_spec.md", + "PRAS Files" => "PRASFiles/index.md", "Tutorials" => [ "PRAS 101 Walkthrough" => "examples/pras_walkthrough.md", "Demand Response Walkthrough" => "examples/demand_response_walkthrough.md", diff --git a/docs/src/PRAS/results.md b/docs/src/PRAS/results.md index cfeac887..02ca3b57 100644 --- a/docs/src/PRAS/results.md +++ b/docs/src/PRAS/results.md @@ -79,6 +79,7 @@ they support. |------------------------------------|-------|-----|--------|--------|---------|------------------| | `Shortfall` | Energy | • | | • | • | • | | `ShortfallSamples` | Energy | • | • | • | • | • | +| `ShortfallEvents` | [See below](@ref shortfall_events) | • | • | • | | | | `DemandResponseShortfall` | Energy | • | | • | • | • | | `DemandResponseShortfallSamples` | Energy | • | • | • | • | • | | `Surplus` | Power | • | | | • | • | @@ -87,14 +88,10 @@ they support. *Table: Regional result specification characteristics.* ### Shortfall -The Shortfall family of result specifications (`Shortfall` and -`ShortfallSamples`) reports on unserved energy -occuring during simulations. As quantifying unserved energy is the core aspect -of resource adequacy analysis, in practice almost every assessment requests a -Shortfall-related result. The basic `Shortfall` specification is most -commonly used and reports average shortfall results, while -`ShortfallSamples` provides more detailed results at the level of -individual simulations (samples). + +The Shortfall family of result specifications reports on unserved energy occurring during simulations. +As quantifying unserved energy is the core aspect of resource adequacy analysis, in practice almost every assessment requests a Shortfall-related result. +The basic `Shortfall` specification reports average shortfall results while `ShortfallSamples` provides results for individual simulation samples and timesteps. If the system includes demand response units, demand response specific shortfall results can also be queried: (`DemandResponseShortfall` and `DemandResponseShortfallSamples`). @@ -106,16 +103,12 @@ Note that the system shortfall object is inclusive of any attributable demand re and therefore the two should never be added together. Similar period and region specific metrics can be obtained to the metrics mentioned below in the system shortfall object. -Shortfall result objects can be indexed into by region, timestep, both region -and timestep, or neither. Indexing on neither (via `result[]`) reports -the total shortfall across all regions and time periods. +`Shortfall` and `ShortfallSamples` result objects can be indexed into by region, timestep, both region and timestep or neither. +Indexing on neither (via `result[]`) reports the total shortfall across all regions and time periods. -Shortfall results are unique among the built-in result types in that the raw -results can also be converted to specific probabilistic risk metrics -(**EUE**, **LOLE**, or **CVAR**). For sampling-based methods, both metric -estimates and the standard error of those estimates are provided. For example, -after assessing the system, metrics across all regions and the full simulation -horizon can be extracted as: +`Shortfall` and `ShortfallSamples` results can also be converted to probabilistic risk metrics such as **EUE**, **LOLE** or **CVAR**. +For sampling-based methods, both metric estimates and their standard errors are provided. +For example, metrics across all regions and the full simulation horizon can be extracted as: ```julia shortfall, = assess(sys, SequentialMonteCarlo(), Shortfall()) @@ -135,6 +128,80 @@ lole_region = LOLE(shortfall, region) eue_region_period = EUE(shortfall, region, period) ``` +#### [Shortfall Events](@id shortfall_events) + +`ShortfallEvents` records each contiguous run of simulation timesteps with positive shortfall as a separate event. +Events are identified independently within each Monte Carlo sample at both system and regional levels. +A timestep without shortfall ends an event; crossing midnight does not. +Events still active at the end of a sample are closed at its final timestep and never continue into another sample. + +Unlike `ShortfallResult`, `ShortfallEventsResult` retains individual event timing and energy. +It records events directly instead of requiring postprocessing of the full `ShortfallSamplesResult` array, reducing storage when events are sparse. + +![Shortfall time series for two samples and two regions with the corresponding system and regional event records](../images/shortfall-events.svg) + +*Illustrative shortfall values and their corresponding internal event records; events are recorded directly during simulation without first storing the displayed time series.* + +System events remain continuous while any region has shortfall, even if the affected region changes. +In sample 1, timesteps 1–3 therefore form one system event despite Region A having two separate events; system event counts are not sums of regional event counts. + +##### Requesting and Accessing Events + +```julia +events, = assess(sys, SequentialMonteCarlo(samples=1000), ShortfallEvents()) +region = first(sys.regions.names) + +system_events = events[1] # System events in sample 1 +regional_events = events[region, 1] # Regional events in sample 1 +events_by_sample = events[region] # Regional event lists for every sample +``` + +Samples without events have empty event lists. +Unlike `ShortfallSamplesResult`, this result does not support timestep indexing or contain the within-event shortfall profile or peak shortfall. +Each event has `start_idx`, `end_idx` and `energy` fields. + +```julia +if !isempty(system_events) + event = first(system_events) + start_timestamp = events.timestamps[event.start_idx] + last_timestamp = events.timestamps[event.end_idx] + duration_periods = event.end_idx - event.start_idx + 1 + duration = duration_periods * step(events.timestamps) +end +``` + +Both boundary indices are inclusive. +`last_timestamp` labels the final timestep included in the event; it is not the timestamp of the following timestep. +For example, an event covering three 15-minute timesteps lasts 45 minutes. + +##### Event Metrics and Units + +```julia +lolev = LOLEv(events) +mean_duration = MeanEventDuration(events) +max_duration = MaxEventDuration(events) +mean_energy = MeanEventEnergy(events) +max_energy = MaxEventEnergy(events) + +regional_lolev = LOLEv(events, region) +regional_mean_energy = MeanEventEnergy(events, region) +``` + +All five metrics support either the whole system or a region name. +`LOLEv` averages event counts across samples, including samples with no events. +The mean duration and energy metrics pool all observed events; the maximum metrics report the largest observed event rather than an average of per-sample maxima. +See [Multi-Metric Resource Adequacy Analyses with PRAS](@ref multi_metric_resource_adequacy) for their interpretation. + +If no events are observed, all five event metrics return zero with zero standard error. +The maximum metrics also report zero standard error because they return observed maxima, not estimates of an expected maximum. + +Duration metric values count simulation timesteps and their display identifies the timestep length and unit. +Energy metrics use the system's configured energy unit, such as MWh or kWh. +The internal `event.energy` field is an unconverted sum of shortfall power values, so it must not be interpreted directly as energy in that unit. +For example, 10 MW of shortfall over two 15-minute timesteps gives an internal sum of 20 and an event energy of 5 MWh after conversion. + +For JSON output, see [Exporting shortfall events](@ref exporting_shortfall_events). + ### Surplus The Surplus family of result specifications (`Surplus` and `SurplusSamples`) reports on excess grid injection capacity (via @@ -260,4 +327,4 @@ variants are available for selecting the category of energy-limited resource sample-level disaggregation. Energy is always disaggregated by timestep and may also be disaggregated by unit (get the state of charge or borrowed load of a single unit) or aggregated across the system (get the sum of states of charge or borrowed loads -of all component devices in the system). \ No newline at end of file +of all component devices in the system). diff --git a/docs/src/PRASFiles/index.md b/docs/src/PRASFiles/index.md index 9ed9f195..335fc698 100644 --- a/docs/src/PRASFiles/index.md +++ b/docs/src/PRASFiles/index.md @@ -3,8 +3,56 @@ PRASFiles.jl provides functionality for reading and writing PRAS-specific file formats, allowing you to save and load resource adequacy data structures. This package enables: + - Loading and saving PRAS system models - Importing data from various formats - Exporting results and system information to standard formats For detailed information on available methods and types, see the [API Reference](./api.md). + +## [Exporting shortfall events](@id exporting_shortfall_events) + +`saveevents` writes a `ShortfallEventsResult` to `pras_event_results.json` inside a timestamped output directory. +See [Shortfall Events](@ref shortfall_events) for event definitions and guidance on choosing a result specification. + +```julia +using PRAS + +sys = PRAS.rts_gmlc() +events, = assess(sys, SequentialMonteCarlo(samples=100, seed=1), ShortfallEvents()) + +summary_directory = saveevents(events, sys, "event_results_summary") +records_directory = saveevents( + events, sys, "event_results_with_records"; + include_events=true, +) +``` + +Both calls export the sample count, system type parameters, system attributes and timestamps. +They also export `lolev`, `mean_event_duration`, `max_event_duration`, `mean_event_energy`, `max_event_energy` and `total_events` at system and regional levels. +Each metric contains `mean` and `stderror` fields; for the maximum metrics, `mean` holds the observed maximum and `stderror` is zero. +`total_events` counts all observed events rather than averaging their counts across samples. + +With the default `include_events=false`, `system_events` and each region's `events` array are empty even when `total_events` is positive. +This omits records from the JSON file only; it does not change what `ShortfallEvents()` records during simulation. +With `include_events=true`, the same summary fields are accompanied by flattened event records. +The JSON uses `system_events` for system records and `region_results`, each containing a region `name` and an `events` array, for regional records. + +For example, an individual record for 10 MW of shortfall over two 15-minute timesteps in a system with energy units of MWh would be: + +```json +{ + "sample_id": 1, + "start_timestamp": "2020-01-01T00:00:00.000+00:00", + "end_timestamp": "2020-01-01T00:15:00.000+00:00", + "duration_periods": 2, + "energy": 5.0 +} +``` + +`sample_id` is one-based and identifies the Monte Carlo sample containing the event. +The timestamps label the first and last included timesteps, so `end_timestamp` is not an exclusive end boundary. +`duration_periods` is a timestep count; its physical duration follows from the timestep length `L` and time unit `T` in `type_params`. +`energy` is total unserved energy converted to the energy unit `E` in `type_params`, not the raw sum stored in the internal event object. + +Event exports are separate from the aggregate `pras_results.json` files written by `saveshortfall` and do not include EUE, LOLE or NEUE fields. diff --git a/docs/src/images/shortfall-events.svg b/docs/src/images/shortfall-events.svg new file mode 100644 index 00000000..3ce6fa3c --- /dev/null +++ b/docs/src/images/shortfall-events.svg @@ -0,0 +1,795 @@ +ShortfallEvents for two samples and two regions + +Single ShortfallEventsResult boundary + + + +Shortfall title + +Shortfall (MW) + +Result title + +ShortfallEventsResult + +Regions label + +regions: Region A, Region B + +Timestamps label + +timestamps (Δt = 1 h): 00:00, 01:00, …, 05:00 + +A heading + +Region A + +A timestep samples + + + + + +Timestep + + + + + +Sample 1 + + + + + +Sample 2 + + + + + +1 + + + + + +2 + + + + + +0 + + + + + +2 + + + + + +0 + + + + + +1 + + + + + +3 + + + + + +3 + + + + + +2 + + + + + +4 + + + + + +0 + + + + + +0 + + + + + +5 + + + + + +0 + + + + + +0 + + + + + +6 + + + + + +0 + + + + + +0 + +B heading + +Region B + +B timestep samples + + + + + +Timestep + + + + + +Sample 1 + + + + + +Sample 2 + + + + + +1 + + + + + +0 + + + + + +0 + + + + + +2 + + + + + +1 + + + + + +0 + + + + + +3 + + + + + +0 + + + + + +3 + + + + + +4 + + + + + +0 + + + + + +0 + + + + + +5 + + + + + +2 + + + + + +4 + + + + + +6 + + + + + +0 + + + + + +4 + +System heading + +System + +System timestep samples + + + + + +Timestep + + + + + +Sample 1 + + + + + +Sample 2 + + + + + +1 + + + + + +2 + + + + + +0 + + + + + +2 + + + + + +1 + + + + + +1 + + + + + +3 + + + + + +3 + + + + + +5 + + + + + +4 + + + + + +0 + + + + + +0 + + + + + +5 + + + + + +2 + + + + + +4 + + + + + +6 + + + + + +0 + + + + + +4 + +System events heading + +system_events[sample] + +System event records + + + + + +Sample + + + + + +start_idx + + + + + +end_idx + + + + + +energy + + + + + +1 + + + + + +1 + + + + + +3 + + + + + +6 + + + + + +1 + + + + + +5 + + + + + +5 + + + + + +2 + + + + + +2 + + + + + +2 + + + + + +3 + + + + + +6 + + + + + +2 + + + + + +5 + + + + + +6 + + + + + +8 + +Regional events heading + +region_events[region, sample] + +Regional event records + + + + + +Region + + + + + +Sample + + + + + +start_idx + + + + + +end_idx + + + + + +energy + + + + + +Region A + + + + + +1 + + + + + +1 + + + + + +1 + + + + + +2 + + + + + +Region A + + + + + +1 + + + + + +3 + + + + + +3 + + + + + +3 + + + + + +Region A + + + + + +2 + + + + + +2 + + + + + +3 + + + + + +3 + + + + + +Region B + + + + + +1 + + + + + +2 + + + + + +2 + + + + + +1 + + + + + +Region B + + + + + +1 + + + + + +5 + + + + + +5 + + + + + +2 + + + + + +Region B + + + + + +2 + + + + + +3 + + + + + +3 + + + + + +3 + + + + + +Region B + + + + + +2 + + + + + +5 + + + + + +6 + + + + + +8 + +Right Arrow 19 + + + + diff --git a/docs/src/resourceadequacy.md b/docs/src/resourceadequacy.md index fd91f598..e0c777d7 100644 --- a/docs/src/resourceadequacy.md +++ b/docs/src/resourceadequacy.md @@ -34,7 +34,11 @@ of event-periods (e.g. event-hours per year, event-days per year). When reported in terms of event-hours, LOLE is sometimes referred to as LOLH (loss-of-load hours). -For a practical comparison of resource adequacy metrics, see [Multi-Metric Resource Adequacy Analyses with PRAS](@ref multi_metric_resource_adequacy). +**Loss-of-Load Events (LOLEv)** is the expected number of contiguous shortfall events over the study period. +An event can contain several consecutive event-periods and can cross a calendar-day boundary. +LOLEv therefore measures event frequency rather than the number of shortfall periods or days, complementing LOLE and EUE. + +For definitions, references and a practical comparison of resource adequacy metrics, see [Multi-Metric Resource Adequacy Analyses with PRAS](@ref multi_metric_resource_adequacy). While a system's shortfall risk can never be eliminated entirely, if these risk metrics are assessed to be lower than some predetermined threshold, the