LOLD Implementation

PRAS.jl/examples/pras_adequacy_metrics.jl

@@ -0,0 +1,167 @@
+# # Interpreting Resource Adequacy Metrics
+#
+# In practice, no single metric fully captures system adequacy. Instead,
+# multiple complementary metrics should be considered together to understand
+# the frequency, distribution and severity of shortfall events.
+# ([NERC (2018)](https://www.nerc.com/globalassets/who-we-are/standing-committees/rstc/pawg/probabilistic_adequacy_and_measures_report.pdf),
+# [EPRI](https://www.epri.com/research/products/3002027833), 
+# [Stephen et al. 2022](https://doi.org/10.1109/PMAPS53380.2022.9810615)).
+#
+# For this reason, PRAS provides multiple result specifications and derived
+# metrics that allow different aspects of system risk to be evaluated
+# consistently.
+
+# ## Event-Based Interpretation
+#
+# A useful way to interpret adequacy metrics is through the concept of
+# **event-periods**.
+#
+# - An **event-period** occurs when shortfall exists in a given simulation time step
+# - An **event-day** occurs when at least one event-period occurs within a day
+#
+# An **adequacy event** is a set of event-periods that are contiguous at the
+# highest available temporal resolution
+# ([Stephen et al. 2022](https://doi.org/10.1109/PMAPS53380.2022.9810615)).
+#
+# This distinction is important because different metrics count different
+# quantities:
+#
+# - **LOLE** counts event-periods
+# - **LOLD** counts event-days
+# - **LOLEv** counts adequacy events
+#
+# These metrics are related, but they are not interchangeable.
+#
+#md # !!! note
+#md #     In PRAS the time resolution of LOLE is determined by the
+#md #     simulation timestamps of the system and is not assumed to always be hourly.
+
+# Another important reason to use multiple metrics, as described in
+# ([Stephen et al. 2022](https://doi.org/10.1109/PMAPS53380.2022.9810615)), 
+# is that systems with similar shortfall magnitudes or counts of event-periods
+# can exhibit very different temporal patterns.
+#
+# We can consider a simple example of two cases next:
+#
+# **Case A**: One day with 10 hours of shortfall
+#
+# **Case B**: Ten days with 1 hour of shortfall each
+#
+# | Metric | Case A | Case B |
+# |------|--------|--------|
+# | LOLE | same | same |
+# | EUE | same | same |
+# | LOLD | 1 | 10 |
+#
+# As we can see in the table above, even though LOLE and EUE are identical in this case, 
+# LOLD reveals that shortfall events are more dispersed in Case B.
+#
+
+# Because event-periods may be distributed across many days, a system with the
+# same number of shortfall periods can have very different numbers of event-days. 
+# As a result, exact conversions between hourly and daily adequacy
+# criteria are not generally possible
+# ([Stephen et al. 2022](https://doi.org/10.1109/PMAPS53380.2022.9810615)).
+
+# This behavior is reflected in PRAS results, where LOLE and LOLD provide
+# complementary views of how shortfall events are distributed in time.
+
+# ## Mathematical Interpretation
+#
+# In PRAS, adequacy metrics can be interpreted from Monte Carlo shortfall
+# samples.
+#
+# Using the following notation:
+#
+# - ``r`` indexes regions
+# - ``t`` indexes timestamps
+# - ``d`` indexes calendar days
+# - ``s`` indexes Monte Carlo samples
+# - ``e`` indexes adequacy events
+# - ``S_{r,t,s}`` denotes the shortfall in region ``r``, at timestamp ``t``,
+#   in Monte Carlo sample ``s``
+# - ``T(d)`` is the set of timestamps in day ``d``
+#
+# the adequacy metrics can be expressed as expectations over Monte Carlo samples:
+#
+# ### LOLE
+#
+# LOLE counts the expected number of event-periods with shortfall:
+
+# ```math
+# \mathrm{LOLE} =
+# \mathbb{E}\left[\sum_t
+# \mathbf{1}\left(\sum_r S_{r,t,s} > 0\right)\right]
+# ```
+#
+#
+# ### LOLD
+#
+# LOLD counts the expected number of days containing at least one shortfall:
+
+# ```math
+# \mathrm{LOLD} = \mathbb{E}\left[\sum_d I_{d,s}\right]
+# ```
+#
+# where:
+
+# ```math
+# I_{d,s} =
+# \begin{cases}
+# 1 & \text{if } \exists t \in T(d) \text{ such that } \sum_r S_{r,t,s} > 0 \\
+# 0 & \text{otherwise}
+# \end{cases}
+# ```
+#
+#
+# ### LOLEv
+#
+# LOLEv counts the expected number of adequacy events:
+
+# ```math
+# \mathrm{LOLEv} = \mathbb{E}\left[\sum_e J_{e,s}\right]
+# ```
+#
+# where:
+
+# ```math
+# J_{e,s} =
+# \begin{cases}
+# 1 & \text{if adequacy event } e \text{ occurs in sample } s \\
+# 0 & \text{otherwise}
+# \end{cases}
+# ```
+
+# ## Analysis with PRAS
+#
+# We revisit the RTS-GMLC with increased system load to induce shortfall
+# which was described in ref(@id pras_walkthrough)
+
+using PRAS
+sys = PRAS.rts_gmlc()
+sys.regions.load .+= 700.0
+
+shortfall_samples, = assess(
+    sys,
+    SequentialMonteCarlo(samples=100, seed=1),
+    ShortfallSamples(),
+)
+
+# And print the metrics we discussed above: 
+println(LOLE(shortfall_samples))
+println(LOLD(shortfall_samples))
+
+# LOLE describes how many simulation periods experience shortfall, while LOLD
+# describes how many days contain at least one such period.
+#
+# In the RTS example above, the system has approximately 85 shortfall hours
+# but only 25.8 shortfall days. This indicates that shortfall events are
+# temporally clustered, meaning that multiple shortfall hours tend to occur within the
+# same day rather than being evenly distributed across the year.
+
+
+# ## References
+#
+# - [NERC (2018), *Probabilistic Adequacy and Measures Technical Reference Report*](https://www.nerc.com/globalassets/who-we-are/standing-committees/rstc/pawg/probabilistic_adequacy_and_measures_report.pdf)
+# - [EPRI, *Resource Adequacy Gap Assessment: Resource Adequacy Assessment Framework*](https://www.epri.com/research/products/3002027833)
+# - [Stephen et al. (2022), *Clarifying the Interpretation and Use of the LOLE Resource Adequacy Metric*](https://doi.org/10.1109/PMAPS53380.2022.9810615)

PRAS.jl/examples/pras_walkthrough.jl

@@ -176,3 +176,28 @@ println("Surplus in")
 # performed on the subset of samples in which the event was observed, using the
 # `ShortfallSamples`, `UtilizationSamples`, and
 # `StorageEnergySamples` result specifications instead.
+
+
+# ## Export Aggregate Results
+
+# After exploring the simulation outputs, we may want to save the
+# aggregate results for reporting or further post-processing.
+
+# Rather than querying individual metrics (e.g., LOLE, EUE, NEUE)
+# one by one, PRAS provides a utility to export all aggregate
+# system-level and region-level results in a single step.
+
+using PRASFiles
+
+output_path = saveshortfall(shortfall, sys, "pras_output");
+println("Results exported to: ", output_path)
+
+# This creates a timestamped directory containing a `pras_results.json`
+# file with:
+# - system-level metrics (LOLE, EUE, NEUE)
+# - region-level metrics
+# - load and capacity summaries
+# - horizon timestamps
+
+# Note that only aggregate results are exported. Sample-level data
+# from the Monte Carlo simulation are not included.

PRASCore.jl/src/Results/Results.jl

@@ -5,14 +5,15 @@ import OnlineStats: Series
 import OnlineStatsBase: EqualWeight, Mean, Variance, value
 import Printf: @sprintf
 import StatsBase: mean, std, stderror
+import Dates: Date
 
 import ..Systems: SystemModel, ZonedDateTime, Period,
                   PowerUnit, EnergyUnit, conversionfactor,
                   unitsymbol, Regions
 export
 
     # Metrics
-    ReliabilityMetric, LOLE, EUE, NEUE,
+    ReliabilityMetric, LOLE, EUE, NEUE, LOLD,
     val, stderror,
 
     # Result specifications
@@ -80,6 +81,15 @@ NEUE(x::AbstractShortfallResult, r::AbstractString, ::Colon) =
 NEUE(x::AbstractShortfallResult, ::Colon, ::Colon) =
     NEUE.(x, x.regions.names, permutedims(x.timestamps))
 
+LOLD(x::AbstractShortfallResult, ::Colon, d::Date) =
+    LOLD.(x, x.regions.names, d)
+
+LOLD(x::AbstractShortfallResult, r::AbstractString, ::Colon) =
+    LOLD.(x, r, _unique_days(x.timestamps))
+
+LOLD(x::AbstractShortfallResult, ::Colon, ::Colon) =
+    LOLD.(x, x.regions.names, permutedims(_unique_days(x.timestamps)))
+
 include("Shortfall.jl")
 include("ShortfallSamples.jl")
 

PRASCore.jl/src/Results/ShortfallSamples.jl

@@ -15,6 +15,7 @@ shortfall, =
     assess(sys, SequentialMonteCarlo(samples=10), ShortfallSamples())
 
 period = ZonedDateTime(2020, 1, 1, 0, tz"UTC")
+day = Date(2020, 1, 1)
 
 samples = shortfall["Region A", period]
 
@@ -25,19 +26,23 @@ samples = shortfall["Region A", period]
 eue = EUE(shortfall)
 lole = LOLE(shortfall)
 neue = NEUE(shortfall)
+lold = LOLD(shortfall)
 
 # Regional risk metrics
 regional_eue = EUE(shortfall, "Region A")
 regional_lole = LOLE(shortfall, "Region A")
 regional_neue = NEUE(shortfall, "Region A")
+regional_lold = LOLD(shortfall, "Region A")
 
 # Period-specific risk metrics
 period_eue = EUE(shortfall, period)
 period_lolp = LOLE(shortfall, period)
+day_lold = LOLD(shortfall, day)
 
 # Region- and period-specific risk metrics
 period_eue = EUE(shortfall, "Region A", period)
 period_lolp = LOLE(shortfall, "Region A", period)
+regional_day_lold = LOLD(shortfall, "Region A", day)
 ```
 
 Note that this result specification requires large amounts of memory for
@@ -194,3 +199,68 @@ function finalize(
         system.regions, system.timestamps, acc.shortfall)
 
 end
+
+function _count_dropped_days_by_sample(flags_by_period_sample, day_ids)
+    nperiods, nsamples = size(flags_by_period_sample)
+    counts = zeros(Int, nsamples)
+
+    isempty(day_ids) && return counts
+
+    for s in 1:nsamples
+        current_day = day_ids[1]
+        day_has_shortfall = false
+        dropped_days = 0
+
+        for t in 1:nperiods
+            if day_ids[t] != current_day
+                dropped_days += day_has_shortfall
+                current_day = day_ids[t]
+                day_has_shortfall = false
+            end
+
+            day_has_shortfall |= flags_by_period_sample[t, s]
+        end
+
+        dropped_days += day_has_shortfall
+        counts[s] = dropped_days
+    end
+
+    return counts
+end
+
+function LOLD(x::ShortfallSamplesResult)
+    flags = dropdims(sum(x.shortfall, dims=1) .> 0, dims=1)
+    daycounts = _count_dropped_days_by_sample(flags, _day_ids(x.timestamps))
+    return LOLD{_ndays(x.timestamps)}(MeanEstimate(daycounts))
+end
+
+function LOLD(x::ShortfallSamplesResult, r::AbstractString)
+    i_r = findfirstunique(x.regions.names, r)
+    flags = Matrix(view(x.shortfall, i_r, :, :) .> 0)
+    daycounts = _count_dropped_days_by_sample(flags, _day_ids(x.timestamps))
+    return LOLD{_ndays(x.timestamps)}(MeanEstimate(daycounts))
+end
+
+function LOLD(x::ShortfallSamplesResult, d::Date)
+    dayrange = _day_range(x.timestamps, d)
+
+    # day_slice has shape: region × day_timestamps × sample
+    day_slice = view(x.shortfall, :, dayrange, :)
+
+    # For each sample, event occurs if any region has any shortfall in any timestep of the day
+    eventdays = vec(dropdims(any(day_slice .> 0, dims=(1, 2)), dims=(1, 2)))
+
+    return LOLD{1}(MeanEstimate(eventdays))
+end
+
+function LOLD(x::ShortfallSamplesResult, r::AbstractString, d::Date)
+    i_r = findfirstunique(x.regions.names, r)
+    dayrange = _day_range(x.timestamps, d)
+
+    region_day_slice = view(x.shortfall, i_r, dayrange, :)
+
+    # For each sample, did this region have any shortfall during the day?
+    eventdays = vec(dropdims(any(region_day_slice .> 0, dims=1), dims=1))
+
+    return LOLD{1}(MeanEstimate(eventdays))
+end

PRASCore.jl/src/Results/metrics.jl

@@ -164,3 +164,32 @@ function Base.show(io::IO, x::NEUE)
     print(io, "NEUE = ", x.neue, " ppm")
 
 end
+
+
+"""
+    LOLD
+
+`LOLD` reports loss of load days 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 LOLD{D} <: ReliabilityMetric
+    lold::MeanEstimate
+
+    function LOLD{D}(lold::MeanEstimate) where {D}
+        val(lold) >= 0 || throw(DomainError(val(lold),
+            "$(val(lold)) is not a valid expected count of event-days"))
+        new{D}(lold)
+    end
+end
+
+val(x::LOLD) = val(x.lold)
+stderror(x::LOLD) = stderror(x.lold)
+
+function Base.show(io::IO, x::LOLD{D}) where {D}
+    print(io, "LOLD = ", x.lold, " event-day/",
+          D == 1 ? "day" : string(D) * "days")
+end