Jf/add demand response component

PRAS.jl/examples/demand_response_walkthrough.jl

@@ -0,0 +1,157 @@
+# # [Demand Response Walkthrough](@id demand_response_walkthrough)  
+
+# This is a complete example of adding demand response to a system,
+# using the [RTS-GMLC](https://github.com/GridMod/RTS-GMLC) system
+
+# Load the PRAS package and other tools necessary for analysis
+using PRAS
+using Plots
+using DataFrames
+using Dates
+using Measures
+
+# ## Add Demand Response to the SystemModel
+
+# For the purposes of this example, we'll just use the built-in RTS-GMLC model. For
+# further information on loading in systems and exploring please see the [PRAS walkthrough](@ref pras_walkthrough).
+rts_gmlc_sys = PRAS.rts_gmlc();
+
+# Lets overview the system information and make sure the demand response we are creating is correct unit wise.
+rts_gmlc_sys
+
+# First, we  define our new demand response component. In accordance with the broader system
+# the component will have a simulation length of 8784 timesteps, hourly interval, and MW/MWh power/energy units.
+# We will then have a single demand response resource of type `"DR_TYPE1"` with a 50 MW borrow and payback capacity,
+# 200 MWh energy capacity, 0% borrowed energy interest, 6 hour allowable payback time periods,
+# 10% outage probability, and 90% recovery probability. The setup below uses the `fill` function to create matrices
+# with the correct dimensions for each of the parameters, which can be extended to multiple demand response resources
+# by changing the `number_of_drs` variable and adjusting names and types accordingly.
+(timesteps,periodlen,periodunit,powerunit,energyunit) = get_params(rts_gmlc_sys);
+number_of_drs = 1;
+new_drs = DemandResponses{timesteps,periodlen,periodunit,powerunit,energyunit}(
+    ["DR1"],
+    ["DR_TYPE1"],
+    fill(50, number_of_drs, timesteps),   # borrow power capacity
+    fill(50, number_of_drs, timesteps),   # payback power capacity
+    fill(200, number_of_drs, timesteps),  # load energy capacity
+    fill(0.0, number_of_drs, timesteps),  # 0% borrowed energy interest
+    fill(6, number_of_drs, timesteps),    # 6 hour allowable payback time periods
+    fill(0.1, number_of_drs, timesteps),  # 10% outage probability
+    fill(0.9, number_of_drs, timesteps),  # 90% recovery probability
+    );
+
+# We will also assign the demand response to region "2" of the system.
+dr_region_indices = [1:0,1:1,2:0];
+
+# We also want to increase the load in the system to see the effect of demand response being utilized.
+# We do this by creating a new load matrix that is 25% higher than the original load.
+updated_load = Int.(round.(rts_gmlc_sys.regions.load .* 1.25));
+
+# Lets define our new regions with the updated load.
+
+new_regions = Regions{timesteps,powerunit}(["1","2","3"],updated_load);
+
+# Finally, we create two new system models, one with dr and one without.
+modified_rts_with_dr  = SystemModel(
+    new_regions, rts_gmlc_sys.interfaces,
+    rts_gmlc_sys.generators, rts_gmlc_sys.region_gen_idxs,
+    rts_gmlc_sys.storages, rts_gmlc_sys.region_stor_idxs,
+    rts_gmlc_sys.generatorstorages, rts_gmlc_sys.region_genstor_idxs,
+    new_drs, dr_region_indices,
+    rts_gmlc_sys.lines, rts_gmlc_sys.interface_line_idxs,
+    rts_gmlc_sys.timestamps);
+
+modified_rts_without_dr  = SystemModel(
+    new_regions, rts_gmlc_sys.interfaces,
+    rts_gmlc_sys.generators, rts_gmlc_sys.region_gen_idxs,
+    rts_gmlc_sys.storages, rts_gmlc_sys.region_stor_idxs,
+    rts_gmlc_sys.generatorstorages, rts_gmlc_sys.region_genstor_idxs,
+    rts_gmlc_sys.lines, rts_gmlc_sys.interface_line_idxs,
+    rts_gmlc_sys.timestamps);
+
+# For validation, we can check that one new demand response device is in the system, and the other system has none.
+println("System with DR\n ",modified_rts_with_dr.demandresponses)
+println("\nSystem without DR\n ",modified_rts_without_dr.demandresponses)
+
+# ## Run a Sequential Monte Carlo Simulation with and without Demand Response
+# We can now run a sequential monte carlo simulation with and without the demand response to see the effect it has on the system.
+# We will query the shortfall attributable to demand response (load that was borrowed and never able to be paid back) and total system shortfall.
+simspec = SequentialMonteCarlo(samples=100, seed=112);
+resultspecs =   (Shortfall(),DemandResponseShortfall(),DemandResponseEnergy());
+
+
+shortfall_with_dr, dr_shortfall_with_dr,dr_energy_with_dr = assess(modified_rts_with_dr, simspec, resultspecs...);
+shortfall_without_dr, dr_shortfall_without_dr,dr_energy_without_dr = assess(modified_rts_without_dr, simspec, resultspecs...);
+
+# Querying the results, we can see that total system shortfall is lower with demand response, across EUE and LOLE metrics.
+println("LOLE Shortfall with DR: ", LOLE(shortfall_with_dr))
+println("LOLE Shortfall without DR: ", LOLE(shortfall_without_dr))
+
+println("\nEUE Shortfall with DR: ", EUE(shortfall_with_dr))
+println("EUE Shortfall without DR: ", EUE(shortfall_without_dr))
+
+# We can also collect the same reliability metrics with the demand response shortfall, which is the amount of load that was borrowed and never able to be paid back.
+println("EUE Demand Response Shortfall: ", EUE(dr_shortfall_with_dr))
+println("LOLE Demand Response Shortfall: ", LOLE(dr_shortfall_with_dr))
+
+# This means over the simulation, load borrowed and unable to be paid back was 250MWh plus or minus 30 MWh.
+# Similarly, we have a loss of load expectation from demand response of 2.1 event hours per year.
+
+# Lets plot the borrowed load of the demand response device over the simulation.
+borrowed_load = [x[1] for x in dr_energy_with_dr["DR1",:]];
+plot(rts_gmlc_sys.timestamps, borrowed_load, xlabel="Timestamp", ylabel="DR1 Borrowed Load", title="DR1 Borrowed Load vs Time", label="")
+
+# We can see that the demand response device was utilized during the summer months, however never borrowing up to its full 200MWh capacity.
+
+# Lets plot the demand response borrowed load across the month and hour of day for greater granularity on when load is being borrowed.
+months = month.(rts_gmlc_sys.timestamps);
+hours = hour.(rts_gmlc_sys.timestamps) .+ 1;
+
+heatmap_matrix = zeros(Float64, 24, 12);
+for (val, m, h) in zip(borrowed_load, months, hours)
+    heatmap_matrix[h, m] += val;
+end
+
+heatmap(
+    1:12, 0:23, heatmap_matrix;
+    xlabel="Month", ylabel="Hour of Day", title="Total DR1 Borrowed Load (MWh)",
+    xticks=(1:12, ["Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"]),
+    colorbar_title="Borrowed Load", color=cgrad([:white, :red], scale = :linear),
+    left_margin = 5mm, right_margin = 5mm
+)
+
+# Maximum borrowed load occurs in the late afternoon July month, with a different peaking pattern as greater surplus exists in August, with reduced load constraints.
+
+# We can also back calculate the borrow power and payback power, by calculating timestep to timestep differences. 
+# Note, payback power here will not capture any dr attributable shortfall or the impact of `borrowed_energy_interest`.
+
+borrow_power = zeros(Float64, timesteps);
+payback_power= zeros(Float64, timesteps);
+borrow_power = max.(0.0, borrowed_load[2:end] .- borrowed_load[1:end-1]);
+payback_power = max.(0.0, borrowed_load[1:end-1] .- borrowed_load[2:end]);
+
+
+# And then plotting the two heatmaps to identify when key borrowing and payback periods are occuring.
+borrow_heatmap = zeros(Float64, 24, 12)
+payback_heatmap = zeros(Float64, 24, 12)
+
+for (b, p, m, h) in zip(borrow_power, payback_power, months[2:end], hours[2:end])
+    borrow_heatmap[h, m] += b
+    payback_heatmap[h, m] += p
+end
+
+p1 = heatmap(1:12, 0:23, borrow_heatmap; ylabel = "Hour of Day", title="DR1 Borrow",
+    xticks=(1:12, ["Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"]),
+    xtickfont=font(7),
+    colorbar_title="Borrow (MW)", color=cgrad([:white, :red]),
+    left_margin = 5mm, right_margin = 3mm);
+p2 = heatmap(1:12, 0:23, payback_heatmap; ylabel = "Hour of Day", title="DR1 Payback",
+    xticks=(1:12, ["Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"]),
+    xtickfont=font(7),
+    colorbar_title="Payback (MW)", color=cgrad([:white, :blue]),
+    left_margin = 3mm, right_margin = 5mm);
+
+plot(p1, p2; layout=(1,2), size=(1000, 500), link = :all)
+
+# We can see peak borrowing occurs around 4-6pm, shifting earlier in the following month, with payback,
+# occurring almost immediately after borrowing, peaking around 7-9pm in July.

PRAS.jl/examples/pras_walkthrough.jl

@@ -1,4 +1,4 @@
-# # [PRAS walkthrough](@id pras_walkthrough)  
+# # [PRAS Walkthrough](@id pras_walkthrough)  
 
 # This is a complete example of running a PRAS assessment,
 # using the [RTS-GMLC](https://github.com/GridMod/RTS-GMLC) system
@@ -9,7 +9,7 @@ using Plots
 using DataFrames
 using Printf
 
-# ## [Read and explore a SystemModel](@id explore_systemmodel)
+# ## [Read and Explore a SystemModel](@id explore_systemmodel)
 
 # You can load in a system model from a [.pras file](@ref prasfile) if you have one like so:
 # ```julia
@@ -30,7 +30,7 @@ sys
 
 # This system has 3 regions, with multiple Generators, one GenerationStorage in 
 # region "2" and one Storage in region "3". We can see regional information by 
-# indexing the system with the region name:
+# indexing the system by region name:
 sys["2"]
 
 # We can visualize a time series of the regional load in region "2":
@@ -45,7 +45,7 @@ system_generators = sys.generators
 
 # This returns an object of the asset type [Generators](@ref PRASCore.Systems.Generators)
 # and we can retrieve capacities of all generators in the system, which returns 
-# a Matrix with the shape (number of generators) x (number of timepoints):
+# a Matrix with the shape (number of generators) x (number of timesteps):
 system_generators.capacity
 
 # We can visualize a time series of the total system capacity 
@@ -77,7 +77,7 @@ region_3_storage = sys["3", Storages]
 # and the generation-storage device in region "2" like so:
 region_2_genstorage = sys["2", GeneratorStorages]
 
-# ## Run a Sequential Monte Carlo simulation
+# ## Run a Sequential Monte Carlo Simulation
 
 # We can run a Sequential Monte Carlo simulation on this system using the
 # [assess](@ref PRASCore.Simulations.assess) function. 
@@ -151,20 +151,20 @@ utilization_str = join([@sprintf("Interface between regions 1 and 2 utilization:
                             utilization["1" => "3", max_lole_ts][1])], "\n");
 println(utilization_str)                            
 
-# We see that the interfaces are not fully utilized, meaning there is 
-# no excess generation in the system that could be wheeled into region "2"
+# We see that the interfaces are not fully utilized, which means there is 
+# no excess generation in the system that could be transferred into region "2"
 # and we can confirm this by looking at the surplus generation in each region
 println("Surplus in")
-println(@sprintf("  region 1: %0.2f",surplus["1",max_lole_ts][1]))
-println(@sprintf("  region 2: %0.2f",surplus["2",max_lole_ts][1]))
-println(@sprintf("  region 3: %0.2f",surplus["3",max_lole_ts][1]))
+@printf("  region 1: %0.2f\n",surplus["1",max_lole_ts][1])
+@printf("  region 2: %0.2f\n",surplus["2",max_lole_ts][1])
+@printf("  region 3: %0.2f\n",surplus["3",max_lole_ts][1])
 
 # Is local storage another alternative for region 3? One can check on the average 
 # state-of-charge of the existing battery in region "3", both in the 
 # hour before and during the problematic period:
 
-println(@sprintf("Storage energy T-1: %0.2f",storage["313_STORAGE_1", max_lole_ts-Hour(1)][1]))
-println(@sprintf("Storage energy T: %0.2f",storage["313_STORAGE_1", max_lole_ts][1]))
+@printf("Storage energy T-1: %0.2f\n",storage["313_STORAGE_1", max_lole_ts-Hour(1)][1])
+@printf("Storage energy T: %0.2f\n",storage["313_STORAGE_1", max_lole_ts][1])
 
 # It may be that the battery is on average charged going in to the event,
 # and perhaps retains some energy during the event, even as load is being

PRASCore.jl/src/Results/DemandResponseAvailability.jl

@@ -0,0 +1,80 @@
+"""
+    DemandResponseAvailability
+
+The `DemandResponseAvailability` result specification reports the sample-level
+discrete availability of `DemandResponses`, producing a `DemandResponseAvailabilityResult`.
+
+A `DemandResponseAvailabilityResult` can be indexed by demand response device name and
+a timestamp to retrieve a vector of sample-level availability states for
+the unit in the given timestep. States are provided as a boolean with
+`true` indicating that the unit is available and `false` indicating that
+it's unavailable.
+
+Example:
+
+```julia
+dravail, =
+    assess(sys, SequentialMonteCarlo(samples=10), DemandResponseAvailability())
+
+samples = dravail["MyDR123", ZonedDateTime(2020, 1, 1, 0, tz"UTC")]
+
+@assert samples isa Vector{Bool}
+@assert length(samples) == 10
+```
+"""
+struct DemandResponseAvailability <: ResultSpec end
+
+struct DRAvailabilityAccumulator <: ResultAccumulator{DemandResponseAvailability}
+
+    available::Array{Bool,3}
+
+end
+
+function accumulator(
+    sys::SystemModel{N}, nsamples::Int, ::DemandResponseAvailability
+) where {N}
+
+    ndrs = length(sys.demandresponses)
+    available = zeros(Bool, ndrs, N, nsamples)
+
+    return DRAvailabilityAccumulator(available)
+
+end
+
+function merge!(
+    x::DRAvailabilityAccumulator, y::DRAvailabilityAccumulator
+)
+
+    x.available .|= y.available
+    return
+
+end
+
+accumulatortype(::DemandResponseAvailability) = DRAvailabilityAccumulator
+
+struct DemandResponseAvailabilityResult{N,L,T<:Period} <: AbstractAvailabilityResult{N,L,T}
+
+    demandresponses::Vector{String}
+    timestamps::StepRange{ZonedDateTime,T}
+
+    available::Array{Bool,3}
+
+end
+
+names(x::DemandResponseAvailabilityResult) = x.demandresponses
+
+function getindex(x::DemandResponseAvailabilityResult, s::AbstractString, t::ZonedDateTime)
+    i_dr = findfirstunique(x.demandresponses, s)
+    i_t = findfirstunique(x.timestamps, t)
+    return vec(x.available[i_dr, i_t, :])
+end
+
+function finalize(
+    acc::DRAvailabilityAccumulator,
+    system::SystemModel{N,L,T,P,E},
+) where {N,L,T,P,E}
+
+    return DemandResponseAvailabilityResult{N,L,T}(
+        system.demandresponses.names, system.timestamps, acc.available)
+
+end