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/12] 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/12] 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/12] 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/12] 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/12] 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/12] 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/12] 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/12] 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/12] 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/12] 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/12] 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/12] 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(