diff --git a/PRASCore.jl/src/Results/Results.jl b/PRASCore.jl/src/Results/Results.jl index 059bae38..01347aa7 100644 --- a/PRASCore.jl/src/Results/Results.jl +++ b/PRASCore.jl/src/Results/Results.jl @@ -13,10 +13,11 @@ 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 + # Result specifications and accessors Shortfall, ShortfallSamples, DemandResponseShortfall, DemandResponseShortfallSamples, Surplus, SurplusSamples, @@ -26,7 +27,8 @@ export DemandResponseEnergy, DemandResponseEnergySamples, GeneratorAvailability, StorageAvailability, GeneratorStorageAvailability,DemandResponseAvailability, - LineAvailability + LineAvailability, ShortfallEvents, + eventsinterval include("metrics.jl") include("utils.jl") @@ -194,6 +196,7 @@ 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..4bc1629d --- /dev/null +++ b/PRASCore.jl/src/Results/ShortfallEvents.jl @@ -0,0 +1,327 @@ +""" + 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. 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 + +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 +_event_meanestimate(xs::AbstractVector{<:Real}) = + isempty(xs) ? MeanEstimate(0.0) : MeanEstimate(xs) + +mutable struct ShortfallEventsAccumulator <: 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, ::ShortfallEvents +) where {N} + + 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( + 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(::ShortfallEvents) = ShortfallEventsAccumulator + +struct ShortfallEventsResult{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} <: Result{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] + +""" + 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))) + +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, + system::SystemModel{N,L,T,P,E}, +) where {N,L,T,P,E} + + + return ShortfallEventsResult{N,L,T,P,E}( + system.regions, system.timestamps, + acc.system_events, acc.region_events) +end + +function MeanEventDuration(x::ShortfallEventsResult{N,L,T}) where {N,L,T} + durations = [ + duration_periods(ev) + for events in x.system_events + for ev in events + ] + 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 = [ + 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}(_event_meanestimate(durations)) +end + +function MaxEventDuration(x::ShortfallEventsResult{N,L,T}) where {N,L,T} + duration = maximum(( + duration_periods(ev) + for events in x.system_events + for ev in events + ); 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) + duration = maximum(( + duration_periods(ev) + for s in axes(x.region_events, 2) + for ev in x.region_events[i_r, s] + ); 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} + p2e = conversionfactor(L, T, P, E) + energies = [ + p2e * event_energy(ev) + for events in x.system_events + for ev in events + ] + 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 = [ + 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}(_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) + energy = maximum(( + p2e * event_energy(ev) + for events in x.system_events + for ev in events + ); 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) + energy = maximum(( + p2e * event_energy(ev) + for s in axes(x.region_events, 2) + for ev in x.region_events[i_r, s] + ); init=0.0) + return MaxEventEnergy{N,L,T,E}(MeanEstimate(energy)) +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..2034b22d 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, + system::SystemModel{N,L,T,P,E}, + state::SystemState, problem::DispatchProblem, + sampleid::Int, t::Int +) where {N,L,T,P,E} + + 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(Results.ShortfallEvents, 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 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/PRASCore.jl/test/Results/shortfall.jl b/PRASCore.jl/test/Results/shortfall.jl index 503d6b72..af10fce7 100644 --- a/PRASCore.jl/test/Results/shortfall.jl +++ b/PRASCore.jl/test/Results/shortfall.jl @@ -261,3 +261,71 @@ 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 + empty!(selected[1]) + @test result[1] == [event(1, 3, 6), event(6, 8, 12)] + + # 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( + 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 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/PRASCore.jl/test/Simulations/runtests.jl b/PRASCore.jl/test/Simulations/runtests.jl index 67c55ab6..ae717f38 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,107 @@ 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) + + 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 = 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) + + 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 = isempty(energies_1a) ? 0.0 : maximum(energies_1a) + + @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 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 + + 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 "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 614a92f1..58c0a56e 100644 --- a/PRASFiles.jl/src/PRASFiles.jl +++ b/PRASFiles.jl/src/PRASFiles.jl @@ -2,12 +2,17 @@ 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, ShortfallResult, ShortfallSamplesResult, - AbstractShortfallResult, Result + AbstractShortfallResult, Result, ShortfallEventsResult, + ShortfallEvent, LOLEv, totalevents, + MeanEventDuration, MaxEventDuration, + 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 @@ -21,6 +26,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..dda63520 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,12 +181,29 @@ 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 -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) @@ -113,6 +214,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; @@ -134,6 +239,45 @@ function get_lold_result( return LOLDResult(shortfall; region = region) end +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 + records[record_idx] = EventRecord( + sample_id, + 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, + region::String, +) + i_r = findfirstunique(events.regions.names, region) + return _get_eventrecords(events, view(events.region_events, i_r, :)) +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..66e5e8f3 100644 --- a/PRASFiles.jl/src/Results/write.jl +++ b/PRASFiles.jl/src/Results/write.jl @@ -104,3 +104,108 @@ function saveshortfall( error("saveshortfall is not implemented for $(typeof(shortfall))") end + +function generate_eventresult( + events::ShortfallEventsResult, + pras_sys::SystemModel; + include_events::Bool = false, +) + + system_event_records = include_events ? + get_eventrecords(events) : + 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( + get_nsamples(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` summary metrics in JSON format, optionally 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 diff --git a/PRASFiles.jl/test/runtests.jl b/PRASFiles.jl/test/runtests.jl index eb6fdc90..30bc29e1 100644 --- a/PRASFiles.jl/test/runtests.jl +++ b/PRASFiles.jl/test/runtests.jl @@ -110,53 +110,163 @@ 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 + 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 diff --git a/docs/src/PRASCore/api.md b/docs/src/PRASCore/api.md index 1686994a..3d68fb05 100644 --- a/docs/src/PRASCore/api.md +++ b/docs/src/PRASCore/api.md @@ -47,4 +47,11 @@ PRASCore.Results.CVAR PRASCore.Results.NCVAR PRASCore.Results.val PRASCore.Results.stderror +PRASCore.Results.ShortfallEvents +PRASCore.Results.eventsinterval +PRASCore.Results.LOLEv +PRASCore.Results.MeanEventDuration +PRASCore.Results.MaxEventDuration +PRASCore.Results.MeanEventEnergy +PRASCore.Results.MaxEventEnergy ```