diff --git a/PRAS.jl/examples/demand_response_walkthrough.jl b/PRAS.jl/examples/demand_response_walkthrough.jl new file mode 100644 index 00000000..9d3fe643 --- /dev/null +++ b/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. \ No newline at end of file diff --git a/PRAS.jl/examples/pras_walkthrough.jl b/PRAS.jl/examples/pras_walkthrough.jl index 359de2ce..081f982a 100644 --- a/PRAS.jl/examples/pras_walkthrough.jl +++ b/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 diff --git a/PRASCore.jl/src/Results/DemandResponseAvailability.jl b/PRASCore.jl/src/Results/DemandResponseAvailability.jl new file mode 100644 index 00000000..fce28d8e --- /dev/null +++ b/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 diff --git a/PRASCore.jl/src/Results/DemandResponseEnergy.jl b/PRASCore.jl/src/Results/DemandResponseEnergy.jl new file mode 100644 index 00000000..0480d688 --- /dev/null +++ b/PRASCore.jl/src/Results/DemandResponseEnergy.jl @@ -0,0 +1,100 @@ +""" + DemandResponseEnergy + +The `DemandResponseEnergy` result specification reports the average energy borrowed +of `DemandResponses`, producing a `DemandResponseEnergyResult`. + +A `DemandResponseEnergyResult` can be indexed by demand response device name and a timestamp to +retrieve a tuple of sample mean and standard deviation, estimating the average +energy borrowed level for the given demand response device in that timestep. + +Example: + +```julia +drenergy, = + assess(sys, SequentialMonteCarlo(samples=1000), DemandResponseEnergy()) + +soc_mean, soc_std = + drenergy["MyDemandResponse123", ZonedDateTime(2020, 1, 1, 0, tz"UTC")] +``` + +See [`DemandResponseEnergySamples`](@ref) for sample-level demand response states of borrowed energy. +""" +struct DemandResponseEnergy <: ResultSpec end + +mutable struct DemandResponseEnergyAccumulator <: ResultAccumulator{DemandResponseEnergy} + + # Cross-simulation energy mean/variances + energy_period::Vector{MeanVariance} + energy_demandresponseperiod::Matrix{MeanVariance} + +end + +function accumulator( + sys::SystemModel{N}, nsamples::Int, ::DemandResponseEnergy +) where {N} + + ndemandresponses = length(sys.demandresponses) + + energy_period = [meanvariance() for _ in 1:N] + energy_demandresponseperiod = [meanvariance() for _ in 1:ndemandresponses, _ in 1:N] + + return DemandResponseEnergyAccumulator( + energy_period, energy_demandresponseperiod) + +end + +function merge!( + x::DemandResponseEnergyAccumulator, y::DemandResponseEnergyAccumulator +) + + foreach(merge!, x.energy_period, y.energy_period) + foreach(merge!, x.energy_demandresponseperiod, y.energy_demandresponseperiod) + + return + +end + +accumulatortype(::DemandResponseEnergy) = DemandResponseEnergyAccumulator + +struct DemandResponseEnergyResult{N,L,T<:Period,E<:EnergyUnit} <: AbstractEnergyResult{N,L,T} + + nsamples::Union{Int,Nothing} + demandresponses::Vector{String} + timestamps::StepRange{ZonedDateTime,T} + + energy_mean::Matrix{Float64} + + energy_period_std::Vector{Float64} + energy_regionperiod_std::Matrix{Float64} + +end + +names(x::DemandResponseEnergyResult) = x.demandresponses + +function getindex(x::DemandResponseEnergyResult, t::ZonedDateTime) + i_t = findfirstunique(x.timestamps, t) + return sum(view(x.energy_mean, :, i_t)), x.energy_period_std[i_t] +end + +function getindex(x::DemandResponseEnergyResult, s::AbstractString, t::ZonedDateTime) + i_dr = findfirstunique(x.demandresponses, s) + i_t = findfirstunique(x.timestamps, t) + return x.energy_mean[i_dr, i_t], x.energy_regionperiod_std[i_dr, i_t] +end + +function finalize( + acc::DemandResponseEnergyAccumulator, + system::SystemModel{N,L,T,P,E}, +) where {N,L,T,P,E} + + _, period_std = mean_std(acc.energy_period) + demandresponseperiod_mean, demandresponseperiod_std = mean_std(acc.energy_demandresponseperiod) + + nsamples = first(first(acc.energy_period).stats).n + + return DemandResponseEnergyResult{N,L,T,E}( + nsamples, system.demandresponses.names, system.timestamps, + demandresponseperiod_mean, period_std, demandresponseperiod_std) + +end diff --git a/PRASCore.jl/src/Results/DemandResponseEnergySamples.jl b/PRASCore.jl/src/Results/DemandResponseEnergySamples.jl new file mode 100644 index 00000000..a9646ef5 --- /dev/null +++ b/PRASCore.jl/src/Results/DemandResponseEnergySamples.jl @@ -0,0 +1,87 @@ +""" + DemandResponseEnergySamples + +The `DemandResponseEnergySamples` result specification reports the sample-level state +of borrowed energy of `DemandResponses`, producing a `DemandResponseEnergySamplesResult`. + +A `DemandResponseEnergySamplesResult` can be indexed by demand response device name and +a timestamp to retrieve a vector of sample-level borrowed energy states for +the device in the given timestep. + +Example: + +```julia +demandresponseenergy, = + assess(sys, SequentialMonteCarlo(samples=10), DemandResponseEnergySamples()) + +samples = demandresponseenergy["MyDemandResponse123", ZonedDateTime(2020, 1, 1, 0, tz"UTC")] + +@assert samples isa Vector{Float64} +@assert length(samples) == 10 +``` + +Note that this result specification requires large amounts of memory for +larger sample sizes. See [`DemandResponseEnergy`](@ref) for estimated average demand response +borrowed energy when sample-level granularity isn't required. +""" +struct DemandResponseEnergySamples <: ResultSpec end + +struct DemandResponseEnergySamplesAccumulator <: ResultAccumulator{DemandResponseEnergySamples} + + energy::Array{Float64,3} + +end + +function accumulator( + sys::SystemModel{N}, nsamples::Int, ::DemandResponseEnergySamples +) where {N} + + ndrs = length(sys.demandresponses) + energy = zeros(Int, ndrs, N, nsamples) + + return DemandResponseEnergySamplesAccumulator(energy) + +end + +function merge!( + x::DemandResponseEnergySamplesAccumulator, y::DemandResponseEnergySamplesAccumulator +) + + x.energy .+= y.energy + return + +end + +accumulatortype(::DemandResponseEnergySamples) = DemandResponseEnergySamplesAccumulator + +struct DemandResponseEnergySamplesResult{N,L,T<:Period,E<:EnergyUnit} <: AbstractEnergyResult{N,L,T} + + demandresponses::Vector{String} + timestamps::StepRange{ZonedDateTime,T} + + energy::Array{Int,3} + +end + +names(x::DemandResponseEnergySamplesResult) = x.demandresponses + +function getindex(x::DemandResponseEnergySamplesResult, t::ZonedDateTime) + i_t = findfirstunique(x.timestamps, t) + return vec(sum(view(x.energy, :, i_t, :), dims=1)) +end + +function getindex(x::DemandResponseEnergySamplesResult, s::AbstractString, t::ZonedDateTime) + i_dr = findfirstunique(x.demandresponses, s) + i_t = findfirstunique(x.timestamps, t) + return vec(x.energy[i_dr, i_t, :]) +end + +function finalize( + acc::DemandResponseEnergySamplesAccumulator, + system::SystemModel{N,L,T,P,E}, +) where {N,L,T,P,E} + + return DemandResponseEnergySamplesResult{N,L,T,E}( + system.demandresponses.names, system.timestamps, acc.energy) + +end diff --git a/PRASCore.jl/src/Results/Results.jl b/PRASCore.jl/src/Results/Results.jl index ad6f4028..3e4e3300 100644 --- a/PRASCore.jl/src/Results/Results.jl +++ b/PRASCore.jl/src/Results/Results.jl @@ -16,12 +16,16 @@ export val, stderror, # Result specifications - Shortfall, ShortfallSamples, Surplus, SurplusSamples, + Shortfall, ShortfallSamples, + DemandResponseShortfall, DemandResponseShortfallSamples, + Surplus, SurplusSamples, Flow, FlowSamples, Utilization, UtilizationSamples, StorageEnergy, StorageEnergySamples, GeneratorStorageEnergy, GeneratorStorageEnergySamples, + DemandResponseEnergy, DemandResponseEnergySamples, GeneratorAvailability, StorageAvailability, - GeneratorStorageAvailability, LineAvailability + GeneratorStorageAvailability,DemandResponseAvailability, + LineAvailability include("utils.jl") include("metrics.jl") @@ -79,6 +83,7 @@ NEUE(x::AbstractShortfallResult, ::Colon, ::Colon) = include("Shortfall.jl") include("ShortfallSamples.jl") + abstract type AbstractSurplusResult{N,L,T} <: Result{N,L,T} end getindex(x::AbstractSurplusResult, ::Colon) = @@ -144,6 +149,7 @@ getindex(x::AbstractAvailabilityResult, ::Colon, ::Colon) = include("GeneratorAvailability.jl") include("StorageAvailability.jl") include("GeneratorStorageAvailability.jl") +include("DemandResponseAvailability.jl") include("LineAvailability.jl") abstract type AbstractEnergyResult{N,L,T} <: Result{N,L,T} end @@ -162,8 +168,10 @@ getindex(x::AbstractEnergyResult, ::Colon, ::Colon) = include("StorageEnergy.jl") include("GeneratorStorageEnergy.jl") +include("DemandResponseEnergy.jl") include("StorageEnergySamples.jl") include("GeneratorStorageEnergySamples.jl") +include("DemandResponseEnergySamples.jl") function resultchannel( results::T, threads::Int diff --git a/PRASCore.jl/src/Results/Shortfall.jl b/PRASCore.jl/src/Results/Shortfall.jl index e76c00c1..0070b39b 100644 --- a/PRASCore.jl/src/Results/Shortfall.jl +++ b/PRASCore.jl/src/Results/Shortfall.jl @@ -43,7 +43,9 @@ See [`ShortfallSamples`](@ref) for recording sample-level shortfall results. """ struct Shortfall <: ResultSpec end -mutable struct ShortfallAccumulator <: ResultAccumulator{Shortfall} +struct DemandResponseShortfall <: ResultSpec end + +mutable struct ShortfallAccumulator{S} <: ResultAccumulator{Shortfall} # Cross-simulation LOL period count mean/variances periodsdropped_total::MeanVariance @@ -68,8 +70,8 @@ mutable struct ShortfallAccumulator <: ResultAccumulator{Shortfall} end function accumulator( - sys::SystemModel{N}, nsamples::Int, ::Shortfall -) where {N} + sys::SystemModel{N}, nsamples::Int, ::S +) where {N,S<:Union{Shortfall,DemandResponseShortfall}} nregions = length(sys.regions) @@ -89,7 +91,7 @@ function accumulator( unservedload_total_currentsim = 0 unservedload_region_currentsim = zeros(Int, nregions) - return ShortfallAccumulator( + return ShortfallAccumulator{S}( periodsdropped_total, periodsdropped_region, periodsdropped_period, periodsdropped_regionperiod, periodsdropped_total_currentsim, periodsdropped_region_currentsim, @@ -117,9 +119,11 @@ function merge!( end -accumulatortype(::Shortfall) = ShortfallAccumulator +accumulatortype(::S) where { + S<:Union{Shortfall,DemandResponseShortfall} + } = ShortfallAccumulator{S} -struct ShortfallResult{N, L, T <: Period, E <: EnergyUnit} <: +struct ShortfallResult{N, L, T <: Period, E <: EnergyUnit, S} <: AbstractShortfallResult{N, L, T} nsamples::Union{Int, Nothing} regions::Regions @@ -145,7 +149,7 @@ struct ShortfallResult{N, L, T <: Period, E <: EnergyUnit} <: shortfall_period_std::Vector{Float64} shortfall_regionperiod_std::Matrix{Float64} - function ShortfallResult{N,L,T,E}( + function ShortfallResult{N,L,T,E,S}( nsamples::Union{Int,Nothing}, regions::Regions, timestamps::StepRange{ZonedDateTime,T}, @@ -162,7 +166,7 @@ struct ShortfallResult{N, L, T <: Period, E <: EnergyUnit} <: shortfall_region_std::Vector{Float64}, shortfall_period_std::Vector{Float64}, shortfall_regionperiod_std::Matrix{Float64} - ) where {N,L,T<:Period,E<:EnergyUnit} + ) where {N,L,T<:Period,E<:EnergyUnit,S <: Union{Shortfall,DemandResponseShortfall}} isnothing(nsamples) || nsamples > 0 || throw(DomainError("Sample count must be positive or `nothing`.")) @@ -184,7 +188,7 @@ struct ShortfallResult{N, L, T <: Period, E <: EnergyUnit} <: size(shortfall_regionperiod_std) == (nregions, N) || error("Inconsistent input data sizes") - new{N,L,T,E}(nsamples, regions, timestamps, + new{N,L,T,E,S}(nsamples, regions, timestamps, eventperiod_mean, eventperiod_std, eventperiod_region_mean, eventperiod_region_std, eventperiod_period_mean, eventperiod_period_std, @@ -206,12 +210,12 @@ function getindex(x::ShortfallResult, r::AbstractString) return sum(view(x.shortfall_mean, i_r, :)), x.shortfall_region_std[i_r] end -function getindex(x::ShortfallResult, t::ZonedDateTime) +function getindex(x::ShortfallResult, t::ZonedDateTime) i_t = findfirstunique(x.timestamps, t) return sum(view(x.shortfall_mean, :, i_t)), x.shortfall_period_std[i_t] end -function getindex(x::ShortfallResult, r::AbstractString, t::ZonedDateTime) +function getindex(x::ShortfallResult, r::AbstractString, t::ZonedDateTime) i_r = findfirstunique(x.regions.names, r) i_t = findfirstunique(x.timestamps, t) return x.shortfall_mean[i_r, i_t], x.shortfall_regionperiod_std[i_r, i_t] @@ -258,19 +262,19 @@ EUE(x::ShortfallResult{N,L,T,E}, t::ZonedDateTime) where {N,L,T,E} = EUE(x::ShortfallResult{N,L,T,E}, r::AbstractString, t::ZonedDateTime) where {N,L,T,E} = EUE{1,L,T,E}(MeanEstimate(x[r, t]..., x.nsamples)) -function NEUE(x::ShortfallResult{N,L,T,E}) where {N,L,T,E} +function NEUE(x::ShortfallResult) return NEUE(div(MeanEstimate(x[]..., x.nsamples),(sum(x.regions.load)/1e6))) end -function NEUE(x::ShortfallResult{N,L,T,E}, r::AbstractString) where {N,L,T,E} +function NEUE(x::ShortfallResult, r::AbstractString) i_r = findfirstunique(x.regions.names, r) return NEUE(div(MeanEstimate(x[r]..., x.nsamples),(sum(x.regions.load[i_r,:])/1e6))) end function finalize( - acc::ShortfallAccumulator, + acc::ShortfallAccumulator{S}, system::SystemModel{N,L,T,P,E}, -) where {N,L,T,P,E} +) where {N,L,T,P,E,S} ep_total_mean, ep_total_std = mean_std(acc.periodsdropped_total) ep_region_mean, ep_region_std = mean_std(acc.periodsdropped_region) @@ -293,7 +297,7 @@ function finalize( ue_period_std .*= p2e ue_regionperiod_std .*= p2e - return ShortfallResult{N,L,T,E}( + return ShortfallResult{N,L,T,E,S}( nsamples, system.regions, system.timestamps, ep_total_mean, ep_total_std, ep_region_mean, ep_region_std, ep_period_mean, ep_period_std, diff --git a/PRASCore.jl/src/Results/ShortfallSamples.jl b/PRASCore.jl/src/Results/ShortfallSamples.jl index 9d30b42b..bddca0f6 100644 --- a/PRASCore.jl/src/Results/ShortfallSamples.jl +++ b/PRASCore.jl/src/Results/ShortfallSamples.jl @@ -44,21 +44,22 @@ Note that this result specification requires large amounts of memory for larger sample sizes. See [`Shortfall`](@ref) for average shortfall outcomes when sample-level granularity isn't required. """ struct ShortfallSamples <: ResultSpec end +struct DemandResponseShortfallSamples <: ResultSpec end -struct ShortfallSamplesAccumulator <: ResultAccumulator{ShortfallSamples} +struct ShortfallSamplesAccumulator{S} <: ResultAccumulator{ShortfallSamples} shortfall::Array{Int,3} end function accumulator( - sys::SystemModel{N}, nsamples::Int, ::ShortfallSamples -) where {N} + sys::SystemModel{N}, nsamples::Int, ::S +) where {N,S<:Union{ShortfallSamples,DemandResponseShortfallSamples}} nregions = length(sys.regions) shortfall = zeros(Int, nregions, N, nsamples) - return ShortfallSamplesAccumulator(shortfall) + return ShortfallSamplesAccumulator{S}(shortfall) end @@ -71,9 +72,11 @@ function merge!( end -accumulatortype(::ShortfallSamples) = ShortfallSamplesAccumulator +accumulatortype(::S) where { + S<:Union{ShortfallSamples,DemandResponseShortfallSamples} + } = ShortfallSamplesAccumulator{S} -struct ShortfallSamplesResult{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} <: AbstractShortfallResult{N,L,T} +struct ShortfallSamplesResult{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit, S} <: AbstractShortfallResult{N,L,T} regions::Regions{N,P} timestamps::StepRange{ZonedDateTime,T} @@ -162,11 +165,11 @@ function NEUE(x::ShortfallSamplesResult{N,L,T,P,E}, r::AbstractString) where {N, end function finalize( - acc::ShortfallSamplesAccumulator, + acc::ShortfallSamplesAccumulator{S}, system::SystemModel{N,L,T,P,E}, -) where {N,L,T,P,E} +) where {N,L,T,P,E,S<:Union{ShortfallSamples,DemandResponseShortfallSamples}} - return ShortfallSamplesResult{N,L,T,P,E}( + return ShortfallSamplesResult{N,L,T,P,E,S}( system.regions, system.timestamps, acc.shortfall) end diff --git a/PRASCore.jl/src/Simulations/DispatchProblem.jl b/PRASCore.jl/src/Simulations/DispatchProblem.jl index 28240c14..63bac7bb 100644 --- a/PRASCore.jl/src/Simulations/DispatchProblem.jl +++ b/PRASCore.jl/src/Simulations/DispatchProblem.jl @@ -2,16 +2,18 @@ DispatchProblem(sys::SystemModel) Create a min-cost flow problem for the multi-region max power delivery problem -with generation and storage discharging in decreasing order of priority, and -storage charging with excess capacity. Storage and GeneratorStorage devices -within a region are represented individually on the network. +with generation and storage discharging in decreasing order of priority, +storage charging with excess capacity and demand response devices enabling load shed and shift functionality. +Storage, GeneratorStorage, and Demand Response devices within a region are represented individually on the network. + This involves injections/withdrawals at one node (regional capacity surplus/shortfall) for each modelled region, as well as two/three nodes associated with each Storage/GeneratorStorage device, and a supplementary "slack" node in the network that can absorb undispatched power or pass unserved energy or unused charging capability through to satisfy -power balance constraints. +power balance constraints. Demand Response devices are represented +in a structurally similar manner as storage charging and discharging. Flows from the generation nodes are free, while flows to charging and from discharging nodes are costed or rewarded according to the @@ -21,7 +23,10 @@ capacity is exhausted (implying an operational strategy that prioritizes resource adequacy over economic arbitrage). This is based on the storage dispatch strategy of Evans, Tindemans, and Angeli, as outlined in "Minimizing Unserved Energy Using Heterogenous Storage Units" (IEEE Transactions on Power -Systems, 2019). +Systems, 2019). Demand Response devices will borrow energy in devices with the +longest payback window first, and vice versa for payback energy. +Demand Response devices are utilized only after discharging all storage/genstor +and paid back before storage/genstor charging. Flows to the charging node have an attenuated negative cost (reward), incentivizing immediate storage charging if generation and transmission @@ -47,7 +52,9 @@ Nodes in the problem are ordered as: 5. GenerationStorage discharge capacity (GeneratorStorage order) 6. GenerationStorage grid injection (GeneratorStorage order) 7. GenerationStorage charge capacity (GeneratorStorage order) - 8. Slack node + 8. DemandResponse payback capacity (DemandResponse order) + 9. DemandResponse borrowing capacity (DemandResponse order) + 10. Slack node Edges are ordered as: @@ -67,6 +74,10 @@ Edges are ordered as: 14. GenerationStorage charge from inflow (GeneratorStorage order) 15. GenerationStorage charge unused (GeneratorStorage order) 16. GenerationStorage inflow unused (GeneratorStorage order) + 17. DemandResponse payback to grid (DemandResponse order) + 18. DemandResponse payback unused (DemandResponse order) + 19. DemandResponse borrowing from grid (DemandResponse order) + 20. DemandResponse borrowing unused (DemandResponse order) """ struct DispatchProblem @@ -80,6 +91,8 @@ struct DispatchProblem genstorage_discharge_nodes::UnitRange{Int} genstorage_togrid_nodes::UnitRange{Int} genstorage_charge_nodes::UnitRange{Int} + demandresponse_payback_nodes::UnitRange{Int} + demandresponse_borrow_nodes::UnitRange{Int} slack_node::Int # Edge labels @@ -99,10 +112,19 @@ struct DispatchProblem genstorage_inflowcharge_edges::UnitRange{Int} genstorage_chargeunused_edges::UnitRange{Int} genstorage_inflowunused_edges::UnitRange{Int} + demandresponse_payback_edges::UnitRange{Int} + demandresponse_paybackunused_edges::UnitRange{Int} + demandresponse_borrow_edges::UnitRange{Int} + demandresponse_borrowunused_edges::UnitRange{Int} + #stor/genstor costs min_chargecost::Int max_dischargecost::Int + #dr costs + max_borrowcost_dr::Int + min_paybackcost_dr::Int + function DispatchProblem( sys::SystemModel; unlimited::Int=999_999_999) @@ -110,14 +132,30 @@ struct DispatchProblem nifaces = length(sys.interfaces) nstors = length(sys.storages) ngenstors = length(sys.generatorstorages) + ndrs = length(sys.demandresponses) + + #Implied merit order: + # dr payback_cost < stor/genstor charge cost < stor/genstor discharge cost < borrow cost < unserved energy + + #for storage/genstor maxchargetime, maxdischargetime = maxtimetocharge_discharge(sys) min_chargecost = - maxchargetime - 1 max_dischargecost = - min_chargecost + maxdischargetime + 1 - shortagepenalty = 10 * (nifaces + max_dischargecost) + + #for demand response-we want to borrow energy in devices with the longest payback window, and payback energy from devices with the smallest payback window + maxpaybacktime_dr = max_payback_window_dr(sys) + cost_for_storage_wheeling_prevention = 75 #prevent storage discharging and wheeling across large number of regions which then results in it being cheaper to borrow than to discharge from storage + min_paybackcost_dr = - max_dischargecost - maxpaybacktime_dr - 1 - cost_for_storage_wheeling_prevention #add max_dischargecost (will mean greater than storage charge as well) to always have DR devices be paybacked first + max_borrowcost_dr = - min_paybackcost_dr + maxpaybacktime_dr + 1 #will always be greater than max_dischargecost and paybacktime, need to add 1 to not overlap with payback reward + + #for unserved energy + shortagepenalty = 10 * (nifaces + max_borrowcost_dr) + stor_regions = assetgrouplist(sys.region_stor_idxs) genstor_regions = assetgrouplist(sys.region_genstor_idxs) + dr_regions = assetgrouplist(sys.region_dr_idxs) region_nodes = 1:nregions stor_discharge_nodes = indices_after(region_nodes, nstors) @@ -126,7 +164,9 @@ struct DispatchProblem genstor_discharge_nodes = indices_after(genstor_inflow_nodes, ngenstors) genstor_togrid_nodes = indices_after(genstor_discharge_nodes, ngenstors) genstor_charge_nodes = indices_after(genstor_togrid_nodes, ngenstors) - slack_node = nnodes = last(genstor_charge_nodes) + 1 + dr_payback_nodes = indices_after(genstor_charge_nodes, ndrs) + dr_borrow_nodes = indices_after(dr_payback_nodes, ndrs) + slack_node = nnodes = last(dr_borrow_nodes) + 1 region_unservedenergy = 1:nregions region_unusedcapacity = indices_after(region_unservedenergy, nregions) @@ -144,7 +184,11 @@ struct DispatchProblem genstor_inflowcharge = indices_after(genstor_gridcharge, ngenstors) genstor_chargeunused = indices_after(genstor_inflowcharge, ngenstors) genstor_inflowunused = indices_after(genstor_chargeunused, ngenstors) - nedges = last(genstor_inflowunused) + dr_paybackused = indices_after(genstor_inflowunused, ndrs) + dr_paybackunused = indices_after(dr_paybackused, ndrs) + dr_borrowused = indices_after(dr_paybackunused, ndrs) + dr_borrowunused = indices_after(dr_borrowused, ndrs) + nedges = last(dr_borrowunused) nodesfrom = Vector{Int}(undef, nedges) nodesto = Vector{Int}(undef, nedges) @@ -196,16 +240,22 @@ struct DispatchProblem initedges(genstor_gridcharge, genstor_regions, genstor_charge_nodes) initedges(genstor_inflowcharge, genstor_inflow_nodes, genstor_charge_nodes) initedges(genstor_chargeunused, slack_node, genstor_charge_nodes) - initedges(genstor_inflowunused, genstor_inflow_nodes, slack_node) + # Demand Response borrowing / payback + initedges(dr_paybackused, dr_regions, dr_payback_nodes) + initedges(dr_paybackunused, slack_node,dr_payback_nodes) + initedges(dr_borrowused, dr_borrow_nodes, dr_regions) + initedges(dr_borrowunused, dr_borrow_nodes, slack_node) + return new( FlowProblem(nodesfrom, nodesto, limits, costs, injections), region_nodes, stor_discharge_nodes, stor_charge_nodes, genstor_inflow_nodes, genstor_discharge_nodes, - genstor_togrid_nodes, genstor_charge_nodes, slack_node, + genstor_togrid_nodes, genstor_charge_nodes, + dr_payback_nodes, dr_borrow_nodes, slack_node, region_unservedenergy, region_unusedcapacity, iface_forward, iface_reverse, @@ -214,7 +264,11 @@ struct DispatchProblem genstor_dischargegrid, genstor_dischargeunused, genstor_inflowgrid, genstor_totalgrid, genstor_gridcharge, genstor_inflowcharge, genstor_chargeunused, - genstor_inflowunused, min_chargecost, max_dischargecost + genstor_inflowunused, + dr_paybackused, dr_paybackunused, + dr_borrowused, dr_borrowunused, + min_chargecost, max_dischargecost, + max_borrowcost_dr, min_paybackcost_dr ) end @@ -242,7 +296,6 @@ function update_problem!( ) - system.regions.load[r, t] updateinjection!(region_node, slack_node, region_netgenavailable) - end # Update bidirectional interface limits (from lines) @@ -275,7 +328,6 @@ function update_problem!( maxenergy = system.storages.energy_capacity[i, t] # Update discharging - maxdischarge = stor_online * system.storages.discharge_capacity[i, t] dischargeefficiency = system.storages.discharge_efficiency[i, t] energydischargeable = stor_energy * dischargeefficiency @@ -380,6 +432,54 @@ function update_problem!( end + + + # Update Demand Response charge/discharge limits and priorities + for (i, (borrow_node, borrow_edge, payback_node, payback_edge)) in + enumerate(zip( + problem.demandresponse_borrow_nodes, problem.demandresponse_borrow_edges, + problem.demandresponse_payback_nodes, problem.demandresponse_payback_edges)) + + dr_online = state.drs_available[i] + dr_energy = state.drs_energy[i] + maxenergy = system.demandresponses.energy_capacity[i, t] + + # Update energy payback + + maxpayback = dr_online * system.demandresponses.payback_capacity[i, t] + paybackefficiency = system.demandresponses.payback_efficiency[i, t] + energy_payback_allowed = dr_energy * paybackefficiency + + allowablepayback = state.drs_paybackcounter[i] + + payback_capacity = min( + maxpayback, floor(Int, energytopower( + energy_payback_allowed, E, L, T, P)) + ) + updateinjection!( + fp.nodes[payback_node], slack_node, -payback_capacity) + + # smallest allowable payback window = highest priority (payback first) + paybackcost = problem.min_paybackcost_dr + allowablepayback + updateflowcost!(fp.edges[payback_edge], paybackcost) + + # Update borrowing-make sure no borrowing is allowed if allowable payback period is equal to zero + maxborrow = system.demandresponses.allowable_payback_period[i,t] != 0 ? dr_online * system.demandresponses.borrow_capacity[i, t] : 0 + borrowefficiency = system.demandresponses.borrow_efficiency[i, t] + energyborrowable = (maxenergy - dr_energy) * borrowefficiency + borrow_capacity = min( + maxborrow, floor(Int, energytopower( + energyborrowable, E, L, T, P)) + ) + updateinjection!( + fp.nodes[borrow_node], slack_node, borrow_capacity) + + # Longest allowable payback window = highest priority (borrow first) + borrowcost = problem.max_borrowcost_dr - allowablepayback # Positive cost + updateflowcost!(fp.edges[borrow_edge], borrowcost) + + end + end function update_state!( @@ -390,6 +490,9 @@ function update_state!( edges = problem.fp.edges p2e = conversionfactor(L, T, P, E) + # Storage Update + + #discharging for (i, e) in enumerate(problem.storage_discharge_edges) energy = state.stors_energy[i] energy_drop = ceil(Int, edges[e].flow * p2e / @@ -397,12 +500,16 @@ function update_state!( state.stors_energy[i] = max(0, energy - energy_drop) end - + + #charging for (i, e) in enumerate(problem.storage_charge_edges) state.stors_energy[i] += ceil(Int, edges[e].flow * p2e * system.storages.charge_efficiency[i, t]) end + # GeneratorStorage Update + + #discharging for (i, e) in enumerate(problem.genstorage_dischargegrid_edges) energy = state.genstors_energy[i] energy_drop = ceil(Int, edges[e].flow * p2e / @@ -410,6 +517,7 @@ function update_state!( state.genstors_energy[i] = max(0, energy - energy_drop) end + #charging for (i, (e1, e2)) in enumerate(zip(problem.genstorage_gridcharge_edges, problem.genstorage_inflowcharge_edges)) totalcharge = (edges[e1].flow + edges[e2].flow) * p2e @@ -417,6 +525,31 @@ function update_state!( ceil(Int, totalcharge * system.generatorstorages.charge_efficiency[i, t]) end + #Demand Response Update + for (i, e) in enumerate(problem.demandresponse_borrow_edges) + state.drs_energy[i] += + ceil(Int, edges[e].flow * p2e * system.demandresponses.borrow_efficiency[i, t]) + end + + #paybacking + for (i, e) in enumerate(problem.demandresponse_payback_edges) + energy = state.drs_energy[i] + energy_drop = ceil(Int, edges[e].flow * p2e / system.demandresponses.payback_efficiency[i, t]) + state.drs_energy[i] = max(0, energy - energy_drop) + + if (state.drs_paybackcounter[i] == 0) || (t == N) + #if payback window is over or you reach the end of the sim, count the unserved energy in drs_unservedenergy state and reset energy + #adding to exisiting unserved energy as we count any unserved energy above soc during borrowed_energy_interest + #we want to keep a running total (will be either zero or the unserved energy above energy capacity) + state.drs_unservedenergy[i] += state.drs_energy[i] + state.drs_energy[i] = 0 + end + end + + + + + end function assetgrouplist(idxss::Vector{UnitRange{Int}}) diff --git a/PRASCore.jl/src/Simulations/Simulations.jl b/PRASCore.jl/src/Simulations/Simulations.jl index b9f220b6..99b70e96 100644 --- a/PRASCore.jl/src/Simulations/Simulations.jl +++ b/PRASCore.jl/src/Simulations/Simulations.jl @@ -175,13 +175,19 @@ function initialize!( rng, state.genstors_available, state.genstors_nexttransition, system.generatorstorages, N) + initialize_availability!( + rng, state.drs_available, state.drs_nexttransition, + system.demandresponses, N) + initialize_availability!( rng, state.lines_available, state.lines_nexttransition, system.lines, N) fill!(state.stors_energy, 0) fill!(state.genstors_energy, 0) - + fill!(state.drs_energy, 0) + fill!(state.drs_unservedenergy, 0) + fill!(state.drs_paybackcounter, -1) return end @@ -204,15 +210,24 @@ function advance!( rng, state.genstors_available, state.genstors_nexttransition, system.generatorstorages, t, N) + update_availability!( + rng, state.drs_available, state.drs_nexttransition, + system.demandresponses, t, N) + update_availability!( rng, state.lines_available, state.lines_nexttransition, system.lines, t, N) update_energy!(state.stors_energy, system.storages, t) update_energy!(state.genstors_energy, system.generatorstorages, t) + update_dr_energy!(state.drs_energy, state.drs_unservedenergy, system.demandresponses, t) + + update_paybackcounter!(state.drs_paybackcounter,state.drs_energy, system.demandresponses,t) + update_problem!(dispatchproblem, state, system, t) + end function solve!( diff --git a/PRASCore.jl/src/Simulations/SystemState.jl b/PRASCore.jl/src/Simulations/SystemState.jl index a88296cb..2e249771 100644 --- a/PRASCore.jl/src/Simulations/SystemState.jl +++ b/PRASCore.jl/src/Simulations/SystemState.jl @@ -11,6 +11,11 @@ struct SystemState genstors_nexttransition::Vector{Int} genstors_energy::Vector{Int} + drs_available::Vector{Bool} + drs_nexttransition::Vector{Int} + drs_energy::Vector{Int} + drs_paybackcounter::Vector{Int} + drs_unservedenergy::Vector{Int} lines_available::Vector{Bool} lines_nexttransition::Vector{Int} @@ -30,6 +35,13 @@ struct SystemState genstors_nexttransition = Vector{Int}(undef, ngenstors) genstors_energy = Vector{Int}(undef, ngenstors) + ndrs = length(system.demandresponses) + drs_available = Vector{Bool}(undef, ndrs) + drs_nexttransition = Vector{Int}(undef, ndrs) + drs_energy = Vector{Int}(undef, ndrs) + drs_unservedenergy = Vector{Int}(undef, ndrs) + drs_paybackcounter = Vector{Int}(undef, ndrs) + nlines = length(system.lines) lines_available = Vector{Bool}(undef, nlines) lines_nexttransition = Vector{Int}(undef, nlines) @@ -38,6 +50,8 @@ struct SystemState gens_available, gens_nexttransition, stors_available, stors_nexttransition, stors_energy, genstors_available, genstors_nexttransition, genstors_energy, + drs_available, drs_nexttransition, drs_energy, + drs_paybackcounter,drs_unservedenergy, lines_available, lines_nexttransition) end diff --git a/PRASCore.jl/src/Simulations/recording.jl b/PRASCore.jl/src/Simulations/recording.jl index 613dda86..a156eb32 100644 --- a/PRASCore.jl/src/Simulations/recording.jl +++ b/PRASCore.jl/src/Simulations/recording.jl @@ -1,22 +1,33 @@ # Shortfall function record!( - acc::Results.ShortfallAccumulator, + acc::Results.ShortfallAccumulator{S}, system::SystemModel{N,L,T,P,E}, state::SystemState, problem::DispatchProblem, sampleid::Int, t::Int -) where {N,L,T,P,E} +) where {N,L,T,P,E,S} totalshortfall = 0 isshortfall = false edges = problem.fp.edges - for r in problem.region_unserved_edges + for (r, dr_idxs) in zip(problem.region_unserved_edges, system.region_dr_idxs) + + #count region shortfall and include dr shortfall + #if shortfall include region and dr shortfall + #if drshortfall don't include region shortfall - regionshortfall = edges[r].flow + 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 + fit!(acc.periodsdropped_regionperiod[r,t], isregionshortfall) fit!(acc.unservedload_regionperiod[r,t], regionshortfall) @@ -29,7 +40,6 @@ function record!( acc.unservedload_region_currentsim[r] += regionshortfall end - end if isshortfall @@ -68,14 +78,22 @@ end # ShortfallSamples function record!( - acc::Results.ShortfallSamplesAccumulator, + acc::Results.ShortfallSamplesAccumulator{S}, system::SystemModel{N,L,T,P,E}, state::SystemState, problem::DispatchProblem, sampleid::Int, t::Int -) where {N,L,T,P,E} +) where {N,L,T,P,E,S} + + for (r,dr_idxs) in zip(problem.region_unserved_edges,system.region_dr_idxs) + #getting dr shortfall + dr_shortfall = 0 + for i in dr_idxs + dr_shortfall += state.drs_unservedenergy[i] + end - for (r, e) in enumerate(problem.region_unserved_edges) - acc.shortfall[r, t, sampleid] = problem.fp.edges[e].flow + acc.shortfall[r, t, sampleid] = + init_regionshortfall(S,problem.fp.edges,r) + + dr_shortfall end return @@ -116,7 +134,6 @@ function record!( regionsurplus += min(grid_limit, total_unused) end - fit!(acc.surplus_regionperiod[r,t], regionsurplus) totalsurplus += regionsurplus @@ -161,7 +178,6 @@ function record!( regionsurplus += min(grid_limit, total_unused) end - acc.surplus[r, t, sampleid] = regionsurplus end @@ -326,6 +342,23 @@ end reset!(acc::Results.GenStorAvailabilityAccumulator, sampleid::Int) = nothing +# DemandResponseAvailability + +function record!( + acc::Results.DRAvailabilityAccumulator, + system::SystemModel{N,L,T,P,E}, + state::SystemState, problem::DispatchProblem, + sampleid::Int, t::Int +) where {N,L,T,P,E} + + acc.available[:, t, sampleid] .= state.drs_available + return + +end + +reset!(acc::Results.DRAvailabilityAccumulator, sampleid::Int) = nothing + + # LineAvailability function record!( @@ -398,6 +431,37 @@ end reset!(acc::Results.GenStorageEnergyAccumulator, sampleid::Int) = nothing + +# DemandResponseEnergy + +function record!( + acc::Results.DemandResponseEnergyAccumulator, + system::SystemModel{N,L,T,P,E}, + state::SystemState, problem::DispatchProblem, + sampleid::Int, t::Int +) where {N,L,T,P,E} + + totalenergy = 0 + ndemandresponses = length(system.demandresponses) + + for s in 1:ndemandresponses + + drenergy = state.drs_energy[s] + fit!(acc.energy_demandresponseperiod[s,t], drenergy) + totalenergy += drenergy + + end + + fit!(acc.energy_period[t], totalenergy) + + return + +end + +reset!(acc::Results.DemandResponseEnergyAccumulator, sampleid::Int) = nothing + + + # StorageEnergySamples function record!( @@ -429,3 +493,20 @@ function record!( end reset!(acc::Results.GenStorageEnergySamplesAccumulator, sampleid::Int) = nothing + + +# DemandResponseEnergySamples + +function record!( + acc::Results.DemandResponseEnergySamplesAccumulator, + system::SystemModel{N,L,T,P,E}, + state::SystemState, problem::DispatchProblem, + sampleid::Int, t::Int +) where {N,L,T,P,E} + + acc.energy[:, t, sampleid] .= state.drs_energy + return + +end + +reset!(acc::Results.DemandResponseEnergySamplesAccumulator, sampleid::Int) = nothing \ No newline at end of file diff --git a/PRASCore.jl/src/Simulations/utils.jl b/PRASCore.jl/src/Simulations/utils.jl index 533c9315..4c9e555a 100644 --- a/PRASCore.jl/src/Simulations/utils.jl +++ b/PRASCore.jl/src/Simulations/utils.jl @@ -120,6 +120,60 @@ function update_energy!( end +function update_dr_energy!( + drs_energy::Vector{Int}, + drs_unservedenergy::Vector{Int}, + drs::AbstractAssets, + t::Int +) + + for i in 1:length(drs_energy) + + soc = drs_energy[i] + borrowed_energy_interest = drs.borrowed_energy_interest[i,t] + 1.0 + maxenergy = drs.energy_capacity[i,t] + + # Decay SoC + soc = round(Int, soc * borrowed_energy_interest) + + unserved_energy_above_max_energy = max(0, soc - maxenergy) + + # Shed SoC above current energy limit + drs_energy[i] = min(soc, maxenergy) + + drs_unservedenergy[i] = unserved_energy_above_max_energy + end + +end + +function update_paybackcounter!( + payback_counter::Vector{Int}, + drs_energy::Vector{Int}, + drs::AbstractAssets, + t::Int +) + + for i in 1:length(payback_counter) + #if energy is zero or negative, set counter to -1 (to start counting new) + if drs_energy[i] <= 0 + if payback_counter[i] >= 0 + #if no energy borrowed and counter is positive, reset it to -1 + payback_counter[i] = -1 + end + elseif payback_counter[i] == -1 + #if energy is borrowed and counter is -1, set it to payback window-start of counting + payback_counter[i] = drs.allowable_payback_period[i,t]-1 + elseif payback_counter[i] >= 0 + #if counter is positive, decrement by one + payback_counter[i] -= 1 + end + + + end + +end + + function maxtimetocharge_discharge(system::SystemModel) if length(system.storages) > 0 @@ -177,6 +231,37 @@ function maxtimetocharge_discharge(system::SystemModel) end +function max_payback_window_dr(system::SystemModel) + + if length(system.demandresponses) > 0 + #no need to check for zero since allowable_payback_period is force to positive + maxpaybacktime_dr = maximum(system.demandresponses.allowable_payback_period) + else + maxpaybacktime_dr = 0 + end + + return maxpaybacktime_dr + +end + +function init_regionshortfall( + ::Type{S}, + edges, + region) where {S <: Union{ + Results.Shortfall, + Results.ShortfallSamples}} + return edges[region].flow +end + +function init_regionshortfall( + ::Type{S}, + edges, + region) where {S <: Union{ + Results.DemandResponseShortfall, + Results.DemandResponseShortfallSamples}} + return 0 +end + function utilization(f::MinCostFlows.Edge, b::MinCostFlows.Edge) flow_forward = f.flow diff --git a/PRASCore.jl/src/Systems/SystemModel.jl b/PRASCore.jl/src/Systems/SystemModel.jl index c9091e9e..65f7c926 100644 --- a/PRASCore.jl/src/Systems/SystemModel.jl +++ b/PRASCore.jl/src/Systems/SystemModel.jl @@ -56,6 +56,9 @@ struct SystemModel{N, L, T <: Period, P <: PowerUnit, E <: EnergyUnit} generatorstorages::GeneratorStorages{N,L,T,P,E} region_genstor_idxs::Vector{UnitRange{Int}} + demandresponses::DemandResponses{N,L,T,P,E} + region_dr_idxs::Vector{UnitRange{Int}} + lines::Lines{N,L,T,P} interface_line_idxs::Vector{UnitRange{Int}} @@ -63,12 +66,14 @@ struct SystemModel{N, L, T <: Period, P <: PowerUnit, E <: EnergyUnit} attrs::Dict{String, String} + #base system constructor-demand response devices function SystemModel{}( regions::Regions{N,P}, interfaces::Interfaces{N,P}, generators::Generators{N,L,T,P}, region_gen_idxs::Vector{UnitRange{Int}}, storages::Storages{N,L,T,P,E}, region_stor_idxs::Vector{UnitRange{Int}}, generatorstorages::GeneratorStorages{N,L,T,P,E}, region_genstor_idxs::Vector{UnitRange{Int}}, + demandresponses::DemandResponses{N,L,T,P,E}, region_dr_idxs::Vector{UnitRange{Int}}, lines::Lines{N,L,T,P}, interface_line_idxs::Vector{UnitRange{Int}}, timestamps::StepRange{ZonedDateTime,T}, attrs::Dict{String, String}=Dict{String, String}() @@ -78,6 +83,7 @@ struct SystemModel{N, L, T <: Period, P <: PowerUnit, E <: EnergyUnit} n_gens = length(generators) n_stors = length(storages) n_genstors = length(generatorstorages) + n_drs = length(demandresponses) n_interfaces = length(interfaces) n_lines = length(lines) @@ -85,6 +91,7 @@ struct SystemModel{N, L, T <: Period, P <: PowerUnit, E <: EnergyUnit} @assert consistent_idxs(region_gen_idxs, n_gens, n_regions) @assert consistent_idxs(region_stor_idxs, n_stors, n_regions) @assert consistent_idxs(region_genstor_idxs, n_genstors, n_regions) + @assert consistent_idxs(region_dr_idxs, n_drs, n_regions) @assert consistent_idxs(interface_line_idxs, n_lines, n_interfaces) @assert all( @@ -97,19 +104,44 @@ struct SystemModel{N, L, T <: Period, P <: PowerUnit, E <: EnergyUnit} new{N,L,T,P,E}( regions, interfaces, generators, region_gen_idxs, storages, region_stor_idxs, - generatorstorages, region_genstor_idxs, lines, interface_line_idxs, + generatorstorages, region_genstor_idxs, + demandresponses, region_dr_idxs, + lines, interface_line_idxs, timestamps, attrs) end end -# No time zone constructor +#base system constructor- no demand response devices +function SystemModel{}( + regions::Regions{N,P}, interfaces::Interfaces{N,P}, + generators::Generators{N,L,T,P}, region_gen_idxs::Vector{UnitRange{Int}}, + storages::Storages{N,L,T,P,E}, region_stor_idxs::Vector{UnitRange{Int}}, + generatorstorages::GeneratorStorages{N,L,T,P,E}, + region_genstor_idxs::Vector{UnitRange{Int}}, + lines::Lines{N,L,T,P}, interface_line_idxs::Vector{UnitRange{Int}}, + timestamps::StepRange{ZonedDateTime,T}, + attrs::Dict{String, String}=Dict{String, String}() +) where {N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} + + return SystemModel( + regions, interfaces, + generators, region_gen_idxs, + storages, region_stor_idxs, + generatorstorages, region_genstor_idxs, + DemandResponses{N,L,T,P,E}(), repeat([1:0],length(regions)), + lines, interface_line_idxs, + timestamps, attrs) +end + +# No time zone constructor - demand responses included function SystemModel( regions::Regions{N,P}, interfaces::Interfaces{N,P}, generators::Generators{N,L,T,P}, region_gen_idxs::Vector{UnitRange{Int}}, storages::Storages{N,L,T,P,E}, region_stor_idxs::Vector{UnitRange{Int}}, generatorstorages::GeneratorStorages{N,L,T,P,E}, region_genstor_idxs::Vector{UnitRange{Int}}, + demandresponses::DemandResponses{N,L,T,P,E}, region_dr_idxs::Vector{UnitRange{Int}}, lines, interface_line_idxs::Vector{UnitRange{Int}}, timestamps::StepRange{DateTime,T}, attrs::Dict{String, String}=Dict{String, String}() @@ -129,37 +161,78 @@ function SystemModel( generators, region_gen_idxs, storages, region_stor_idxs, generatorstorages, region_genstor_idxs, + demandresponses, region_dr_idxs, lines, interface_line_idxs, timestamps_tz,attrs) end -# Single-node constructor +# No time zone constructor - demand responses not included +function SystemModel( + regions::Regions{N,P}, interfaces::Interfaces{N,P}, + generators::Generators{N,L,T,P}, region_gen_idxs::Vector{UnitRange{Int}}, + storages::Storages{N,L,T,P,E}, region_stor_idxs::Vector{UnitRange{Int}}, + generatorstorages::GeneratorStorages{N,L,T,P,E}, region_genstor_idxs::Vector{UnitRange{Int}}, + lines, interface_line_idxs::Vector{UnitRange{Int}}, + timestamps::StepRange{DateTime,T}, + attrs::Dict{String, String}=Dict{String, String}() +) where {N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} + + return SystemModel( + regions, interfaces, + generators, region_gen_idxs, + storages, region_stor_idxs, + generatorstorages, region_genstor_idxs, + DemandResponses{N,L,T,P,E}(), repeat([1:0],length(regions)), + lines, interface_line_idxs, + timestamps,attrs) + +end + +# Single-node constructor - demand responses included function SystemModel( generators::Generators{N,L,T,P}, storages::Storages{N,L,T,P,E}, generatorstorages::GeneratorStorages{N,L,T,P,E}, + demandresponses::DemandResponses{N,L,T,P,E}, timestamps::StepRange{<:AbstractDateTime,T}, load::Vector{Int}, attrs::Dict{String, String}=Dict{String, String}() ) where {N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} return SystemModel( - Regions{N,P}(["Region"], reshape(load, 1, :)), - Interfaces{N,P}( - Int[], Int[], - Matrix{Int}(undef, 0, N), Matrix{Int}(undef, 0, N)), + Regions{N,P}(load), + Interfaces{N,P}(), generators, [1:length(generators)], storages, [1:length(storages)], generatorstorages, [1:length(generatorstorages)], - Lines{N,L,T,P}( - String[], String[], - Matrix{Int}(undef, 0, N), Matrix{Int}(undef, 0, N), - Matrix{Float64}(undef, 0, N), Matrix{Float64}(undef, 0, N)), + demandresponses, [1:length(demandresponses)], + Lines{N,L,T,P}(), UnitRange{Int}[], timestamps, attrs) end +# Single-node constructor - demand responses not included +function SystemModel( + generators::Generators{N,L,T,P}, + storages::Storages{N,L,T,P,E}, + generatorstorages::GeneratorStorages{N,L,T,P,E}, + timestamps::StepRange{<:AbstractDateTime,T}, + load::Vector{Int}, + attrs::Dict{String, String}=Dict{String, String}() +) where {N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} + return SystemModel( + Regions{N,P}(load), + Interfaces{N,P}(), + generators, [1:length(generators)], + storages, [1:length(storages)], + generatorstorages, [1:length(generatorstorages)], + DemandResponses{N,L,T,P,E}(), [1:0], + Lines{N,L,T,P}(), + UnitRange{Int}[], timestamps, attrs) +end + + Base.:(==)(x::T, y::T) where {T <: SystemModel} = x.regions == y.regions && x.interfaces == y.interfaces && @@ -169,6 +242,8 @@ Base.:(==)(x::T, y::T) where {T <: SystemModel} = x.region_stor_idxs == y.region_stor_idxs && x.generatorstorages == y.generatorstorages && x.region_genstor_idxs == y.region_genstor_idxs && + x.demandresponses == y.demandresponses && + x.region_dr_idxs == y.region_dr_idxs && x.lines == y.lines && x.interface_line_idxs == y.interface_line_idxs && x.timestamps == y.timestamps && @@ -181,6 +256,7 @@ unitsymbol(::SystemModel{N,L,T,P,E}) where { unitsymbol(T), unitsymbol(P), unitsymbol(E) isnonnegative(x::Real) = x >= 0 +ispositive(x::Real) = x > 0 isfractional(x::Real) = 0 <= x <= 1 get_params(::SystemModel{N,L,T,P,E}) where {N,L,T,P,E} = @@ -206,6 +282,7 @@ function Base.show(io::IO, sys::SystemModel{N,L,T,P,E}) where {N,L,T<:Period,P<: print(io, "SystemModel($(length(sys.regions)) regions, $(length(sys.interfaces)) interfaces, ", "$(length(sys.generators)) generators, $(length(sys.storages)) storages, ", "$(length(sys.generatorstorages)) generator-storages,", + "$(length(sys.demandresponses)) demand-responses,", " $(N) $(time_unit)s)") end @@ -217,6 +294,7 @@ function Base.show(io::IO, ::MIME"text/plain", sys::SystemModel{N,L,T,P,E}) wher println(io, " Generators: $(length(sys.generators)) units") println(io, " Storages: $(length(sys.storages)) units") println(io, " GeneratorStorages: $(length(sys.generatorstorages)) units") + println(io, " DemandResponses: $(length(sys.demandresponses)) units") println(io, " Lines: $(length(sys.lines))") println(io, "\nTime series:") println(io, " Start time: $(first(sys.timestamps))") @@ -241,6 +319,7 @@ struct RegionInfo generators::NamedTuple storages::NamedTuple generatorstorages::NamedTuple + demandresponses::NamedTuple peak_load::Int power_unit::String end @@ -261,7 +340,8 @@ function Base.getindex(sys::SystemModel, region_idx::Int) gen_range = sys.region_gen_idxs[region_idx] stor_range = sys.region_stor_idxs[region_idx] genstor_range = sys.region_genstor_idxs[region_idx] - + dr_range = sys.region_dr_idxs[region_idx] + # Get peak load for this region peak_load = maximum(sys.regions.load[region_idx, :]) @@ -281,6 +361,10 @@ function Base.getindex(sys::SystemModel, region_idx::Int) indices = genstor_range, count = length(genstor_range), ), + ( + indices = dr_range, + count = length(dr_range), + ), peak_load, power_unit, ) @@ -304,7 +388,8 @@ function Base.show(io::IO, info::RegionInfo) println(io, " Peak load: $(info.peak_load) $(info.power_unit)") println(io, " Generators: $(info.generators.count) units [indices - $(info.generators.indices)]") println(io, " Storages: $(info.storages.count) units [indices - $(info.storages.indices)]") - print(io, " GeneratorStorages: $(info.generatorstorages.count) units [indices - $(info.generatorstorages.indices)]") + println(io, " GeneratorStorages: $(info.generatorstorages.count) units [indices - $(info.generatorstorages.indices)]") + println(io, " DemandResponses: $(info.demandresponses.count) units [indices - $(info.demandresponses.indices)]") end """ @@ -349,7 +434,12 @@ function _get_asset_by_type(sys::SystemModel, region_idx::Int, ::Type{GeneratorS return sys.generatorstorages[genstor_range] end +function _get_asset_by_type(sys::SystemModel, region_idx::Int, ::Type{DemandResponses}) + dr_range = sys.region_dr_idxs[region_idx] + return sys.demandresponses[dr_range] +end + # Fallback for unsupported types function _get_asset_by_type(sys::SystemModel, region_idx::Int, ::Type{T}) where {T} - error("Asset type $T is not supported. Supported types are: Generators, Storages, GeneratorStorages") + error("Asset type $T is not supported. Supported types are: Generators, Storages, GeneratorStorages, DemandResponses") end diff --git a/PRASCore.jl/src/Systems/Systems.jl b/PRASCore.jl/src/Systems/Systems.jl index 6b4cd8b1..405e9a80 100644 --- a/PRASCore.jl/src/Systems/Systems.jl +++ b/PRASCore.jl/src/Systems/Systems.jl @@ -12,7 +12,7 @@ export # System assets Regions, Interfaces, - AbstractAssets, Generators, Storages, GeneratorStorages, Lines, + AbstractAssets, Generators, Storages, GeneratorStorages,DemandResponses, Lines, # Units Period, Minute, Hour, Day, Year, diff --git a/PRASCore.jl/src/Systems/TestData.jl b/PRASCore.jl/src/Systems/TestData.jl index bf9d004c..6b5ba264 100644 --- a/PRASCore.jl/src/Systems/TestData.jl +++ b/PRASCore.jl/src/Systems/TestData.jl @@ -26,6 +26,7 @@ emptygenstors1 = GeneratorStorages{4,1,Hour,MW,MWh}( (empty_int(4) for _ in 1:3)..., (empty_float(4) for _ in 1:3)..., (empty_int(4) for _ in 1:3)..., (empty_float(4) for _ in 1:2)...) + singlenode_a = SystemModel( gens1, emptystors1, emptygenstors1, ZonedDateTime(2010,1,1,0,tz):Hour(1):ZonedDateTime(2010,1,1,3,tz), @@ -53,6 +54,7 @@ emptygenstors1_5min = GeneratorStorages{4,5,Minute,MW,MWh}( (empty_int(4) for _ in 1:3)..., (empty_float(4) for _ in 1:3)..., (empty_int(4) for _ in 1:3)..., (empty_float(4) for _ in 1:2)...) + singlenode_a_5min = SystemModel( gens1_5min, emptystors1_5min, emptygenstors1_5min, ZonedDateTime(2010,1,1,0,0,tz):Minute(5):ZonedDateTime(2010,1,1,0,15,tz), @@ -80,6 +82,7 @@ emptygenstors2 = GeneratorStorages{6,1,Hour,MW,MWh}( (empty_int(6) for _ in 1:3)..., (empty_float(6) for _ in 1:3)..., (empty_int(6) for _ in 1:3)..., (empty_float(6) for _ in 1:2)...) + genstors2 = GeneratorStorages{6,1,Hour,MW,MWh}( ["Genstor1", "Genstor2"], ["Genstorage", "Genstorage"], fill(0, 2, 6), fill(0, 2, 6), fill(4, 2, 6), @@ -141,7 +144,7 @@ lines = Lines{4,1,Hour,MW}( threenode = SystemModel( regions, interfaces, generators, [1:2, 3:5, 6:8], - emptystors1, fill(1:0, 3), emptygenstors1, fill(1:0, 3), + emptystors1, fill(1:0, 3), emptygenstors1, fill(1:0, 3), lines, [1:1, 2:2, 3:3], ZonedDateTime(2018,10,30,0,tz):Hour(1):ZonedDateTime(2018,10,30,3,tz)) @@ -173,6 +176,7 @@ emptygenstors = GeneratorStorages{1,1,Hour,MW,MWh}( (zeros(Int, 0, 1) for _ in 1:3)..., (zeros(Float64, 0, 1) for _ in 1:3)..., (zeros(Int, 0, 1) for _ in 1:3)..., (zeros(Float64, 0, 1) for _ in 1:2)...) + interfaces = Interfaces{1,MW}([1], [2], fill(8, 1, 1), fill(8, 1, 1)) lines = Lines{1,1,Hour,MW}( @@ -216,6 +220,7 @@ emptygenstors = GeneratorStorages{2,1,Hour,MW,MWh}( (zeros(Int, 0, 2) for _ in 1:3)..., (zeros(Float64, 0, 2) for _ in 1:3)..., (zeros(Int, 0, 2) for _ in 1:3)..., (zeros(Float64, 0, 2) for _ in 1:2)...) + test2 = SystemModel(gen, stor, emptygenstors, timestamps, [8, 9]) test2_lole = 0.2 @@ -267,4 +272,172 @@ test3_util_t = [0.8614, 0.626674] test3_eenergy = [6.561, 7.682202] -end + +# Test System 4 (Gen + DR, 1 Region for 6 hours) +# Copper plage, testing DR interactions in one system, that borrowing is occurring and payback is happening +timestamps = ZonedDateTime(2020,1,1,1, tz):Hour(1):ZonedDateTime(2020,1,1,6, tz) + +gen = Generators{6,1,Hour,MW}( + ["Gen 1"], ["Generators"], + fill(57, 1, 6), fill(0.1, 1, 6), fill(0.9, 1, 6)) + +emptystors = Storages{6,1,Hour,MW,MWh}( + (String[] for _ in 1:2)..., + (zeros(Int, 0, 6) for _ in 1:3)..., + (zeros(Float64, 0, 6) for _ in 1:5)...) + +emptygenstors = GeneratorStorages{6,1,Hour,MW,MWh}( + (String[] for _ in 1:2)..., + (zeros(Int, 0, 6) for _ in 1:3)..., (zeros(Float64, 0, 6) for _ in 1:3)..., + (zeros(Int, 0, 6) for _ in 1:3)..., (zeros(Float64, 0, 6) for _ in 1:2)...) + +dr = DemandResponses{6,1,Hour,MW,MWh}( + ["DR1"], ["DemandResponse Category"], + fill(10, 1, 6), fill(10, 1, 6), fill(10, 1, 6), + fill(0.0, 1, 6), + fill(2, 1, 6), fill(0.1, 1, 6), fill(0.9, 1, 6)) + + +full_day_load_profile = [56,58,60,61,59,53] + + +test4 = SystemModel(gen, emptystors, emptygenstors, dr, timestamps, full_day_load_profile) + +test4_lole = 2.118 +test4_lolps = [0.09979000000000159, 0.26288399999999845, 0.3300120000000098, 0.8603499999999678, 0.26384700000001193, 0.30136599999999447] + + +test4_eue = 42.14 +test4_eues = [4.689609999999829, 5.1521070000000115, 6.940174999999926, 12.988998999999957, 6.031562999999943, 6.333204999999944] + + +test4_esurplus = [0.9,0,0,0,0,1.9818] + +test4_eenergy = [0.89863, 2.45616, 4.2544, 0.997662, 2.673, 0.0] + + +# Test System 5 (Gen + DR + Stor, 1 Region for 6 hours) +#Copper plate, testing DR and storage interactions in one system +timestamps = ZonedDateTime(2020,1,1,1, tz):Hour(1):ZonedDateTime(2020,1,1,6, tz) + +gen = Generators{6,1,Hour,MW}( + ["Gen 1"], ["Generators"], + fill(57, 1, 6), fill(0.1, 1, 6), fill(0.9, 1, 6)) + +stor = Storages{6,1,Hour,MW,MWh}( + ["Stor 1"], ["Storages"], + fill(5, 1, 6), fill(5, 1, 6), fill(5, 1, 6), + fill(1., 1, 6), fill(1., 1, 6), fill(1., 1, 6), fill(0.1, 1, 6), fill(0.9, 1, 6)) + +emptygenstors = GeneratorStorages{6,1,Hour,MW,MWh}( + (String[] for _ in 1:2)..., + (zeros(Int, 0, 6) for _ in 1:3)..., (zeros(Float64, 0, 6) for _ in 1:3)..., + (zeros(Int, 0, 6) for _ in 1:3)..., (zeros(Float64, 0, 6) for _ in 1:2)...) + +dr = DemandResponses{6,1,Hour,MW,MWh}( + ["DR1"], ["DemandResponse Category"], + fill(10, 1, 6), fill(10, 1, 6), fill(10, 1, 6), + fill(0.0, 1, 6), + fill(2, 1, 6), fill(0.1, 1, 6), fill(0.9, 1, 6)) + + +full_day_load_profile = [56,58,60,61,59,53] + + +test5 = SystemModel(gen, stor, emptygenstors, dr, timestamps, full_day_load_profile) + +test5_lole = 2.007 +test5_lolps = [0.09979000000000159, 0.19739399999999457, 0.33007500000000567, 0.4260809999999895, 0.698582999999988, 0.25501299999998867] + + + +test5_eue = 42.11 +test5_eues = [4.6888800000001805, 5.013817000000116, 6.877069999999737, 8.325742999999783, 11.226304000000445, 5.983083000000084] + + +test5_esurplus = [0.0901899999999984,0,0,0,0,0.27479499999999374] + +test5_eenergy = [0.89936, 1.86623, 3.66255, 5.05573, 1.54738, 0.0] + + +# Multiregion with DR +#testing multiple DR interactions and that intra dispatch priority is working properly +regions = Regions{6,MW}(["Region 1", "Region 2", "Region 3"], + [19 20 25 26 24 25; 20 21 30 27 23 24; 22 26 27 25 23 24]) + +generators = Generators{6,1,Hour,MW}( + ["Gen1", "VG A", "Gen 2", "Gen 3", "VG B", "Gen 4", "Gen 5", "VG C"], + ["Gens", "VG", "Gens", "Gens", "VG", "Gens", "Gens", "VG"], + [10 10 10 10 10 10; 4 3 2 3 4 3; # A + 10 10 10 10 10 10; 10 10 10 10 10 10; 6 5 3 4 3 2; # B + 10 10 15 10 10 10; 20 20 25 20 22 24; 2 1 2 1 2 2], # C + fill(0.1, 8, 6), + fill(0.9, 8, 6) +) + +drs = DemandResponses{6,1,Hour,MW,MWh}( + ["DR1", "DR2", "DR3"], + ["DR_TYPE1", "DR_TYPE1", "DR_TYPE1"], + [fill(5, 1, 6); # A borrow capacity + fill(4, 1, 6); # B borrow capacity + fill(3, 1, 6);], # C borrow capacity + [fill(5, 1, 6); # A payback capacity + fill(4, 1, 6); # B payback capacity + fill(3, 1, 6);], # C payback capacity + [fill(10, 1, 6); # A energy capacity + fill(8, 1, 6); # B energy capacity + fill(6, 1, 6);], # C energy capacity + fill(0.0, 3, 6), # All regions 0% borrowed energy interest + fill(4, 3, 6), # All regions 3 allowable payback time periods + [fill(0.1, 1, 6); # A + fill(0.1, 1, 6); # B + fill(0.1, 1, 6)], # C + [fill(0.9, 1, 6); # A + fill(0.9, 1, 6); # B + fill(0.9, 1, 6)]) # C) + +interfaces = Interfaces{6,MW}( + [1,1,2], [2,3,3], fill(100, 3, 6), fill(100, 3, 6)) + +lines = Lines{6,1,Hour,MW}( + ["L1", "L2", "L3"], ["Lines", "Lines", "Lines"], + fill(100, 3, 6), fill(100, 3, 6), fill(0.1, 3, 6), fill(0.9, 3, 6)) + +threenode_dr = + SystemModel( + regions, interfaces, generators, [1:2, 3:5, 6:8], + emptystors, fill(1:0, 3), emptygenstors, fill(1:0, 3), + drs, [1:1, 2:2, 3:3], + lines, [1:1, 2:2, 3:3], + ZonedDateTime(2018,10,30,0,tz):Hour(1):ZonedDateTime(2018,10,30,5,tz)) + +threenode_dr_lole = 3.204734 +threenode_dr_lole_r = [2.748671; 2.078730; 1.581805] +threenode_dr_lole_t = [0.0817; 0.2169519; 0.47371299; 0.843717; 0.588652; 1.0] +threenode_dr_lole_rt = [0.01415000; 0.082619; 0.39146799; 0.360588; 0.279216; 0.9506899] + +threenode_dr_eue = 53.0449649 +threenode_dr_eue_r = [26.40881199; 15.280649; 11.35548099] +threenode_dr_eue_t = [0.566655; 1.82072; 5.47627; 9.7399670; 8.7754709; 26.66588199] +threenode_dr_eue_rt = [0.147045; 0.5105039; 2.39299799; 6.236285; 5.0194399; 12.102539] + +threenode_dr_esurplus_t = [6.066344; 0.805918; 0.148378; 0.03628699; 0.0952889; 0.10477099] +threenode_dr_esurplus_rt = [2.67869699; 0.5570409; 0.0; 0.0; 0.0; 0.0] + +threenode_dr_flow = -1.0203119 +threenode_dr_flow_t = [-1.4839049; -2.3727269; -0.2560769; -0.6444219; -0.7271149; -0.6376259] +threenode_dr_util = 0.23285 +threenode_dr_util_t = [0.115079209;0.12469228;0.1083558;0.106332479;0.10810597; 0.10776672] + +threenode_dr_energy = 57.750022 + +threenode_dr_energy_samples = 58.403 + +threenode_dr_shortfall_specific_eue = 23.524975 +threenode_dr_shortfall_specific_eue_r = [10.196; 7.944; 5.385] +threenode_dr_shortfall_specific_eue_t = [0.00; 0.00; 0.00; 0.00; 3.582; 19.943] +threenode_dr_shortfall_specific_eue_rt = [0.0000; 0.0000; 0.0000; 0.0000; 1.895; 8.300] + +threenode_dr_shortfall_samples = 1895487 + +end \ No newline at end of file diff --git a/PRASCore.jl/src/Systems/assets.jl b/PRASCore.jl/src/Systems/assets.jl index 04ee3469..a794a473 100644 --- a/PRASCore.jl/src/Systems/assets.jl +++ b/PRASCore.jl/src/Systems/assets.jl @@ -91,6 +91,15 @@ struct Generators{N,L,T<:Period,P<:PowerUnit} <: AbstractAssets{N,L,T,P} end +# Empty Generators constructor +function Generators{N,L,T,P}() where {N,L,T,P} + + return Generators{N,L,T,P}( + String[], String[], zeros(Int, 0, N), + zeros(Float64, 0, N), zeros(Float64, 0, N)) + +end + Base.:(==)(x::T, y::T) where {T <: Generators} = x.names == y.names && x.categories == y.categories && @@ -226,6 +235,16 @@ struct Storages{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} <: AbstractAssets{N,L, end +# Empty Storages constructor +function Storages{N,L,T,P,E}() where {N,L,T,P,E} + + return Storages{N,L,T,P,E}( + String[], String[], + zeros(Int, 0, N), zeros(Int, 0, N), zeros(Int, 0, N), + zeros(Float64, 0, N), zeros(Float64, 0, N), zeros(Float64, 0, N), + zeros(Float64, 0, N), zeros(Float64, 0, N)) +end + Base.:(==)(x::T, y::T) where {T <: Storages} = x.names == y.names && x.categories == y.categories && @@ -408,6 +427,18 @@ struct GeneratorStorages{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} <: AbstractAs end +# Empty GeneratorStorages constructor +function GeneratorStorages{N,L,T,P,E}() where {N,L,T,P,E} + + return GeneratorStorages{N,L,T,P,E}( + String[], String[], + zeros(Int, 0, N), zeros(Int, 0, N), zeros(Int, 0, N), + zeros(Float64, 0, N), zeros(Float64, 0, N), zeros(Float64, 0, N), + zeros(Int, 0, N), zeros(Int, 0, N), zeros(Int, 0, N), + zeros(Float64, 0, N), zeros(Float64, 0, N)) + +end + Base.:(==)(x::T, y::T) where {T <: GeneratorStorages} = x.names == y.names && x.categories == y.categories && @@ -487,6 +518,206 @@ function Base.vcat(gen_stors::GeneratorStorages{N,L,T,P,E}...) where {N, L, T, P end +""" + DemandResponses{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} + +A struct representing demand response devices in the system. + +# Type Parameters +- `N`: Number of timesteps in the system model +- `L`: Length of each timestep in T units +- `T`: The time period type used for temporal representation, subtype of `Period` +- `P`: The power unit used for capacity measurements, subtype of `PowerUnit` +- `E`: The energy unit used for energy storage, subtype of `EnergyUnit` + +# Fields + - `names`: Name of demand response device + - `categories`: Category of demand response device + - `borrow_capacity`: Maximum available borrowing capacity for each demand response unit in each + timeperiod, expressed in units given by the `power_units` (`P`) type parameter + - `payback_capacity`: Maximum available payback capacity for each demand response unit in + each timeperiod, expressed in units given by the `power_units` (`P`) type parameter + - `energy_capacity`: Maximum available energy capable of being held for each demand response unit in + each timeperiod, expressed in units given by the `energy_units` (`E`) type parameter + - `borrowed_energy_interest`: Growth or decay rate of borrowed energy in the demand response device + from the beginning of one period to energy retained at the end of the previous period, for + each demand response unit in each timeperiod. A `borrowed_energy_interest` of 0.0 has no growth or decay. Unitless. + - `allowable_payback_period`: Maximum number of time steps a demand response device can hold borrowed load. + Any energy still contained at the end of the period will be counted as unserved load for that hour. + If borrowed load is paid back before the end of the `allowable_payback_period`, counter is reset upn further use. (`T`) type parameter + - `λ` (failure probability): Probability the unit transitions from operational to forced + outage during a given simulation timestep, for each storage unit in each timeperiod. + Unitless. + - `μ` (repair probability): Probability the unit transitions from forced outage to + operational during a given simulation timestep, for each storage unit in each + timeperiod. Unitless. +""" +struct DemandResponses{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} <: AbstractAssets{N,L,T,P} + + names::Vector{String} + categories::Vector{String} + + borrow_capacity::Matrix{Int} # power + payback_capacity::Matrix{Int} # power + energy_capacity::Matrix{Int} # energy + + borrowed_energy_interest::Matrix{Float64} + + allowable_payback_period::Matrix{Int} + + λ::Matrix{Float64} + μ::Matrix{Float64} + + borrow_efficiency::Matrix{Float64} + payback_efficiency::Matrix{Float64} + + function DemandResponses{N,L,T,P,E}( + names::Vector{<:AbstractString}, categories::Vector{<:AbstractString}, + borrowcapacity::Matrix{Int}, paybackcapacity::Matrix{Int}, + energycapacity::Matrix{Int}, borrowedenergyinterest::Matrix{Float64}, + allowablepaybackperiod::Matrix{Int}, + λ::Matrix{Float64}, μ::Matrix{Float64}, + borrowefficiency::Matrix{Float64},paybackefficiency::Matrix{Float64} + ) where {N,L,T,P,E} + + n_drs = length(names) + @assert length(categories) == n_drs + @assert allunique(names) + + + @assert size(borrowcapacity) == (n_drs, N) + @assert size(paybackcapacity) == (n_drs, N) + @assert size(energycapacity) == (n_drs, N) + @assert all(isnonnegative, borrowcapacity) + @assert all(isnonnegative, paybackcapacity) + @assert all(isnonnegative, energycapacity) + + @assert size(borrowefficiency) == (n_drs, N) + @assert size(paybackefficiency) == (n_drs, N) + @assert size(borrowedenergyinterest) == (n_drs, N) + @assert all(isfractional, borrowefficiency) + @assert all(isfractional, paybackefficiency) + @assert all(borrowedenergyinterest .<= 1.0) + @assert all(borrowedenergyinterest .>= -1.0) + + @assert size(allowablepaybackperiod) == (n_drs, N) + @assert all(ispositive, allowablepaybackperiod) + + + @assert size(λ) == (n_drs, N) + @assert size(μ) == (n_drs, N) + @assert all(isfractional, λ) + @assert all(isfractional, μ) + + new{N,L,T,P,E}(string.(names), string.(categories), + borrowcapacity, paybackcapacity, energycapacity, + borrowedenergyinterest, + allowablepaybackperiod, + λ, μ,borrowefficiency, paybackefficiency,) + end +end + +# second constructor if borrow and payback efficiencies are not provided +function DemandResponses{N,L,T,P,E}( + names::Vector{<:AbstractString}, categories::Vector{<:AbstractString}, + borrowcapacity::Matrix{Int}, paybackcapacity::Matrix{Int}, + energycapacity::Matrix{Int}, borrowedenergyinterest::Matrix{Float64}, + allowablepaybackperiod::Matrix{Int}, + λ::Matrix{Float64}, μ::Matrix{Float64} +) where {N,L,T,P,E} + return DemandResponses{N,L,T,P,E}( + names, categories, + borrowcapacity, paybackcapacity, energycapacity, + borrowedenergyinterest, allowablepaybackperiod, + λ, μ, + ones(Float64, size(borrowcapacity)), ones(Float64, size(paybackcapacity)) + ) +end + + +# Empty DemandResponses constructor +function DemandResponses{N,L,T,P,E}() where {N,L,T,P,E} + + return DemandResponses{N,L,T,P,E}( + String[], String[], + Matrix{Int}(undef, 0, N),Matrix{Int}(undef, 0, N),Matrix{Int}(undef, 0, N), + Matrix{Float64}(undef, 0, N), + Matrix{Int}(undef, 0, N),Matrix{Float64}(undef, 0, N),Matrix{Float64}(undef, 0, N), + Matrix{Float64}(undef, 0, N),Matrix{Float64}(undef, 0, N)) +end + +Base.:(==)(x::T, y::T) where {T <: DemandResponses} = + x.names == y.names && + x.categories == y.categories && + x.borrow_capacity == y.borrow_capacity && + x.payback_capacity == y.payback_capacity && + x.energy_capacity == y.energy_capacity && + x.borrow_efficiency == y.borrow_efficiency && + x.payback_efficiency == y.payback_efficiency && + x.borrowed_energy_interest == y.borrowed_energy_interest && + x.allowable_payback_period == y.allowable_payback_period && + x.λ == y.λ && + x.μ == y.μ + +Base.getindex(dr::DR, idxs::AbstractVector{Int}) where {DR <: DemandResponses} = + DR(dr.names[idxs], dr.categories[idxs],dr.borrow_capacity[idxs,:], + dr.payback_capacity[idxs, :],dr.energy_capacity[idxs, :], + dr.borrowed_energy_interest[idxs, :],dr.allowable_payback_period[idxs, :],dr.λ[idxs, :], dr.μ[idxs, :], + dr.borrow_efficiency[idxs, :], dr.payback_efficiency[idxs, :]) + +function Base.vcat(drs::DemandResponses{N,L,T,P,E}...) where {N, L, T, P, E} + + n_drs = sum(length(dr) for dr in drs) + + names = Vector{String}(undef, n_drs) + categories = Vector{String}(undef, n_drs) + + borrow_capacity = Matrix{Int}(undef, n_drs, N) + payback_capacity = Matrix{Int}(undef, n_drs, N) + energy_capacity = Matrix{Int}(undef, n_drs, N) + + borrow_efficiency = Matrix{Float64}(undef, n_drs, N) + payback_efficiency = Matrix{Float64}(undef, n_drs, N) + borrowed_energy_interest = Matrix{Float64}(undef, n_drs, N) + + allowable_payback_period = Matrix{Int}(undef, n_drs, N) + + + λ = Matrix{Float64}(undef, n_drs, N) + μ = Matrix{Float64}(undef, n_drs, N) + + last_idx = 0 + + for dr in drs + + n = length(dr) + rows = last_idx .+ (1:n) + + names[rows] = dr.names + categories[rows] = dr.categories + + borrow_capacity[rows, :] = dr.borrow_capacity + payback_capacity[rows, :] = dr.payback_capacity + energy_capacity[rows, :] = dr.energy_capacity + + borrow_efficiency[rows, :] = dr.borrow_efficiency + payback_efficiency[rows, :] = dr.payback_efficiency + borrowed_energy_interest[rows, :] = dr.borrowed_energy_interest + + allowable_payback_period[rows, :] = dr.allowable_payback_period + + λ[rows, :] = dr.λ + μ[rows, :] = dr.μ + + last_idx += n + + end + + return DemandResponses{N,L,T,P,E}(names, categories, borrow_capacity, payback_capacity, energy_capacity, + borrowed_energy_interest,allowable_payback_period, λ, μ, borrow_efficiency, payback_efficiency) + +end + """ Lines{N,L,T<:Period,P<:PowerUnit} @@ -549,6 +780,13 @@ struct Lines{N,L,T<:Period,P<:PowerUnit} <: AbstractAssets{N,L,T,P} end +# Empty Lines constructor +function Lines{N,L,T,P}() where {N,L,T,P} + return Lines{N,L,T,P}( + String[], String[], zeros(Int, 0, N), zeros(Int, 0, N), + zeros(Float64, 0, N), zeros(Float64, 0, N)) +end + Base.:(==)(x::T, y::T) where {T <: Lines} = x.names == y.names && x.categories == y.categories && diff --git a/PRASCore.jl/src/Systems/collections.jl b/PRASCore.jl/src/Systems/collections.jl index a265ae94..60fe4939 100644 --- a/PRASCore.jl/src/Systems/collections.jl +++ b/PRASCore.jl/src/Systems/collections.jl @@ -32,6 +32,11 @@ struct Regions{N,P<:PowerUnit} end +# Single Regions constructor +function Regions{N,P}(load::Vector{Int}) where {N,P} + return Regions{N,P}(["Region"], reshape(load, 1, :)) +end + Base.:(==)(x::T, y::T) where {T <: Regions} = x.names == y.names && x.load == y.load @@ -81,6 +86,12 @@ struct Interfaces{N,P<:PowerUnit} end +# Empty Interfaces constructor +function Interfaces{N,P}() where {N,P} + return Interfaces{N,P}( + Int[], Int[], zeros(Int, 0, N), zeros(Int, 0, N)) +end + Base.:(==)(x::T, y::T) where {T <: Interfaces} = x.regions_from == y.regions_from && x.regions_to == y.regions_to && diff --git a/PRASCore.jl/test/Results/availability.jl b/PRASCore.jl/test/Results/availability.jl index fd06020e..df75424c 100644 --- a/PRASCore.jl/test/Results/availability.jl +++ b/PRASCore.jl/test/Results/availability.jl @@ -41,6 +41,19 @@ @test_throws BoundsError result[r_bad, t] @test_throws BoundsError result[r_bad, t_bad] + # DemandResponses + + result = PRASCore.Results.DemandResponseAvailabilityResult{N,1,Hour}( + DD.resourcenames, DD.periods, available) + + @test length(result[r, t]) == DD.nsamples + @test result[r, t] ≈ vec(available[r_idx, t_idx, :]) + + @test_throws BoundsError result[r, t_bad] + @test_throws BoundsError result[r_bad, t] + @test_throws BoundsError result[r_bad, t_bad] + + # Lines result = PRASCore.Results.LineAvailabilityResult{N,1,Hour}( diff --git a/PRASCore.jl/test/Results/energy.jl b/PRASCore.jl/test/Results/energy.jl index 15f620f1..98d6ab89 100644 --- a/PRASCore.jl/test/Results/energy.jl +++ b/PRASCore.jl/test/Results/energy.jl @@ -34,6 +34,22 @@ @test_throws BoundsError result[r_bad, t] @test_throws BoundsError result[r_bad, t_bad] + # DemandResponses + + result = PRASCore.Results.DemandResponseEnergyResult{N,1,Hour,MWh}( + DD.nsamples, DD.resourcenames, DD.periods, + DD.d1_resourceperiod, DD.d2_period, DD.d2_resourceperiod) + + @test result[t] ≈ (sum(DD.d1_resourceperiod[:, t_idx]), DD.d2_period[t_idx]) + @test result[r, t] ≈ + (DD.d1_resourceperiod[r_idx, t_idx], DD.d2_resourceperiod[r_idx, t_idx]) + + @test_throws BoundsError result[t_bad] + @test_throws BoundsError result[r, t_bad] + @test_throws BoundsError result[r_bad, t] + @test_throws BoundsError result[r_bad, t_bad] + + end @testset "EnergySamplesResult" begin @@ -74,4 +90,20 @@ end @test_throws BoundsError result[r_bad, t] @test_throws BoundsError result[r_bad, t_bad] + # DemandResponses + + result = PRASCore.Results.DemandResponseEnergySamplesResult{N,1,Hour,MWh}( + DD.resourcenames, DD.periods, DD.d) + + @test length(result[t]) == DD.nsamples + @test result[t] ≈ vec(sum(view(DD.d, :, t_idx, :), dims=1)) + + @test length(result[r, t]) == DD.nsamples + @test result[r, t] ≈ vec(DD.d[r_idx, t_idx, :]) + + @test_throws BoundsError result[t_bad] + @test_throws BoundsError result[r, t_bad] + @test_throws BoundsError result[r_bad, t] + @test_throws BoundsError result[r_bad, t_bad] + end diff --git a/PRASCore.jl/test/Results/shortfall.jl b/PRASCore.jl/test/Results/shortfall.jl index 9c6d73d1..9f5d4a01 100644 --- a/PRASCore.jl/test/Results/shortfall.jl +++ b/PRASCore.jl/test/Results/shortfall.jl @@ -5,7 +5,7 @@ r, r_idx, r_bad = DD.testresource, DD.testresource_idx, DD.notaresource t, t_idx, t_bad = DD.testperiod, DD.testperiod_idx, DD.notaperiod - result = PRASCore.Results.ShortfallResult{N,1,Hour,MWh}( + result = PRASCore.Results.ShortfallResult{N,1,Hour,MWh,Shortfall}( DD.nsamples, Regions{N,MW}(DD.resourcenames, DD.resource_vals), DD.periods, DD.d1, DD.d2, DD.d1_resource, DD.d2_resource, DD.d1_period, DD.d2_period, DD.d1_resourceperiod, DD.d2_resourceperiod, @@ -101,7 +101,7 @@ end r, r_idx, r_bad = DD.testresource, DD.testresource_idx, DD.notaresource t, t_idx, t_bad = DD.testperiod, DD.testperiod_idx, DD.notaperiod - result = PRASCore.Results.ShortfallSamplesResult{N,1,Hour,MW,MWh}( + result = PRASCore.Results.ShortfallSamplesResult{N,1,Hour,MW,MWh,ShortfallSamples}( Regions{N,MW}(DD.resourcenames, DD.resource_vals), DD.periods, DD.d) # Overall diff --git a/PRASCore.jl/test/Simulations/runtests.jl b/PRASCore.jl/test/Simulations/runtests.jl index 014f11b3..3c518634 100644 --- a/PRASCore.jl/test/Simulations/runtests.jl +++ b/PRASCore.jl/test/Simulations/runtests.jl @@ -9,7 +9,8 @@ simspec = SequentialMonteCarlo(samples=100_000, seed=1, threaded=false) smallsample = SequentialMonteCarlo(samples=10, seed=123, threaded=false) - resultspecs = (Shortfall(), Surplus(), Flow(), Utilization(), + resultspecs = (Shortfall(), Surplus(), + Flow(), Utilization(), ShortfallSamples(), SurplusSamples(), FlowSamples(), UtilizationSamples(), GeneratorAvailability()) @@ -46,11 +47,11 @@ shortfall2_3, _, flow2_3, util2_3, _ = assess(TestData.threenode, simspec, resultspecs...) - assess(TestData.threenode, smallsample, + assess(TestData.threenode, smallsample, DemandResponseShortfall(), GeneratorAvailability(), LineAvailability(), - StorageAvailability(), GeneratorStorageAvailability(), - StorageEnergy(), GeneratorStorageEnergy(), - StorageEnergySamples(), GeneratorStorageEnergySamples()) + StorageAvailability(), GeneratorStorageAvailability(),DemandResponseAvailability(), + StorageEnergy(), GeneratorStorageEnergy(),DemandResponseEnergy(), + StorageEnergySamples(), GeneratorStorageEnergySamples(),DemandResponseEnergySamples()) @testset "Shortfall Results" begin @@ -409,6 +410,185 @@ TestData.test3_eenergy, simspec.nsamples, nstderr_tol)) + + end + + @testset "Test System 4: Gen + DR, 1 Region" begin + + simspec = SequentialMonteCarlo(samples=1_000_000, seed=112) + region = first(TestData.test4.regions.names) + dr = first(TestData.test4.demandresponses.names) + dts = TestData.test4.timestamps + + shortfall, surplus, energy = + assess(TestData.test4, simspec, + Shortfall(), Surplus(), DemandResponseEnergy()) + + # Shortfall - LOLE + @test withinrange(LOLE(shortfall), + TestData.test4_lole, nstderr_tol) + @test withinrange(LOLE(shortfall, region), + TestData.test4_lole, nstderr_tol) + @test all(withinrange.(LOLE.(shortfall, dts), + TestData.test4_lolps, nstderr_tol)) + @test all(withinrange.(LOLE.(shortfall, region, dts), + TestData.test4_lolps, nstderr_tol)) + + # Shortfall - EUE + @test withinrange(EUE(shortfall), + TestData.test4_eue, nstderr_tol) + @test withinrange(EUE(shortfall, region), + TestData.test4_eue, nstderr_tol) + @test all(withinrange.(EUE.(shortfall, dts), + TestData.test4_eues, nstderr_tol)) + @test all(withinrange.(EUE.(shortfall, region, dts), + TestData.test4_eues, nstderr_tol)) + + # Surplus + @test all(withinrange.(getindex.(surplus, dts), + TestData.test4_esurplus, + simspec.nsamples, nstderr_tol)) + @test all(withinrange.(getindex.(surplus, region, dts), + TestData.test4_esurplus, + simspec.nsamples, nstderr_tol)) + # Energy + @test all(withinrange.(getindex.(energy, dts), + TestData.test4_eenergy, + simspec.nsamples, nstderr_tol)) + @test all(withinrange.(getindex.(energy, dr, dts), + TestData.test4_eenergy, + simspec.nsamples, nstderr_tol)) + end + @testset "Test System 5: Gen + DR + Stor, 1 Region" begin + + simspec = SequentialMonteCarlo(samples=1_000_000, seed=112) + region = first(TestData.test5.regions.names) + dr = first(TestData.test5.demandresponses.names) + dts = TestData.test5.timestamps + + shortfall, surplus, energy = + assess(TestData.test5, simspec, + Shortfall(), Surplus(), DemandResponseEnergy()) + + # Shortfall - LOLE + @test withinrange(LOLE(shortfall), + TestData.test5_lole, nstderr_tol) + @test withinrange(LOLE(shortfall, region), + TestData.test5_lole, nstderr_tol) + @test all(withinrange.(LOLE.(shortfall, dts), + TestData.test5_lolps, nstderr_tol)) + @test all(withinrange.(LOLE.(shortfall, region, dts), + TestData.test5_lolps, nstderr_tol)) + + # Shortfall - EUE + @test withinrange(EUE(shortfall), + TestData.test5_eue, nstderr_tol) + @test withinrange(EUE(shortfall, region), + TestData.test5_eue, nstderr_tol) + @test all(withinrange.(EUE.(shortfall, dts), + TestData.test5_eues, nstderr_tol)) + @test all(withinrange.(EUE.(shortfall, region, dts), + TestData.test5_eues, nstderr_tol)) + + # Surplus + @test all(withinrange.(getindex.(surplus, dts), + TestData.test5_esurplus, + simspec.nsamples, nstderr_tol)) + @test all(withinrange.(getindex.(surplus, region, dts), + TestData.test5_esurplus, + simspec.nsamples, nstderr_tol)) + # Energy + @test all(withinrange.(getindex.(energy, dts), + TestData.test5_eenergy, + simspec.nsamples, nstderr_tol)) + @test all(withinrange.(getindex.(energy, dr, dts), + TestData.test5_eenergy, + simspec.nsamples, nstderr_tol)) + + end + + @testset "Test System 6: Gen + DR, 3 Regions" begin + + simspec = SequentialMonteCarlo(samples=1_000_000, seed=112) + regions = TestData.threenode_dr.regions.names + dr = first(TestData.threenode_dr.demandresponses.names) + dts = TestData.threenode_dr.timestamps + + shortfall, surplus, dr_energy, dr_energy_samples, dr_shortfall, dr_shortfall_samples, flow, utilization = + assess(TestData.threenode_dr, simspec, + Shortfall(), Surplus(), + DemandResponseEnergy(),DemandResponseEnergySamples(), + DemandResponseShortfall(), DemandResponseShortfallSamples(), + Flow(), Utilization()) + + + # Shortfall - LOLE + @test withinrange(LOLE(shortfall), + TestData.threenode_dr_lole, nstderr_tol) + @test all(withinrange.(LOLE.(shortfall, regions), + TestData.threenode_dr_lole_r, nstderr_tol)) + @test all(withinrange.(LOLE.(shortfall, dts), + TestData.threenode_dr_lole_t, nstderr_tol)) + @test all(withinrange.(LOLE.(shortfall, "Region 2", dts), + TestData.threenode_dr_lole_rt, nstderr_tol)) + + # Shortfall - EUE + @test withinrange(EUE(shortfall), + TestData.threenode_dr_eue, nstderr_tol) + @test all(withinrange.(EUE.(shortfall, regions), + TestData.threenode_dr_eue_r, nstderr_tol)) + @test all(withinrange.(EUE.(shortfall, dts), + TestData.threenode_dr_eue_t, nstderr_tol)) + @test all(withinrange.(EUE.(shortfall, "Region 1", dts), + TestData.threenode_dr_eue_rt, nstderr_tol)) + + # Surplus + @test all(withinrange.(getindex.(surplus, dts), + TestData.threenode_dr_esurplus_t, + simspec.nsamples, nstderr_tol)) + @test all(withinrange.(getindex.(surplus, "Region 2", dts), + TestData.threenode_dr_esurplus_rt, + simspec.nsamples, nstderr_tol)) + + # Flow + @test withinrange(flow["Region 1" => "Region 2"], + TestData.threenode_dr_flow, + simspec.nsamples, nstderr_tol) + @test all(withinrange.(getindex.(flow, "Region 1"=>"Region 2", dts), + TestData.threenode_dr_flow_t, + simspec.nsamples, nstderr_tol)) + + # Utilization + @test withinrange((sum(x[1] for x in getindex.(utilization, "Region 1"=>"Region 2")), sum(x[end] for x in getindex.(utilization, "Region 1"=>"Region 2"))), + TestData.threenode_dr_util, + simspec.nsamples, nstderr_tol) + + @test all(withinrange.(getindex.(utilization, "Region 1"=>"Region 2",dts), + TestData.threenode_dr_util_t, + simspec.nsamples, nstderr_tol)) + # DR Energy + @test withinrange((sum(x[1] for x in getindex.(dr_energy, dts)), sum(x[end] for x in getindex.(dr_energy, dts))), + TestData.threenode_dr_energy, + simspec.nsamples, nstderr_tol) + + # DR Energy Samples + @test isapprox(mean(sum(dr_energy_samples.energy[:, :, i]) for i in 1:1_000),TestData.threenode_dr_energy_samples, rtol=0.01) + + # DR Shortfall + @test withinrange(EUE(dr_shortfall), + TestData.threenode_dr_shortfall_specific_eue, nstderr_tol) + @test all(withinrange.(EUE.(dr_shortfall, regions), + TestData.threenode_dr_shortfall_specific_eue_r, nstderr_tol)) + @test all(withinrange.(EUE.(dr_shortfall, dts), + TestData.threenode_dr_shortfall_specific_eue_t, nstderr_tol)) + @test all(withinrange.(EUE.(dr_shortfall, "Region 1", dts), + TestData.threenode_dr_shortfall_specific_eue_rt, nstderr_tol)) + + # DR Shortfall Samples + @test isapprox(sum(dr_shortfall_samples["Region 1",dts[5]])/1e4,TestData.threenode_dr_shortfall_samples/1e4, rtol=0.01) + end + + end diff --git a/PRASCore.jl/test/Systems/SystemModel.jl b/PRASCore.jl/test/Systems/SystemModel.jl index 741b114b..0bf5f577 100644 --- a/PRASCore.jl/test/Systems/SystemModel.jl +++ b/PRASCore.jl/test/Systems/SystemModel.jl @@ -17,11 +17,18 @@ rand(1:10, 1, 10), rand(1:10, 1, 10), rand(1:10, 1, 10), fill(0.1, 1, 10), fill(0.5, 1, 10)) + demandresponses = DemandResponses{10,1,Hour,MW,MWh}( + ["S1", "S2"], ["HVAC", "Industrial"], + rand(1:10, 2, 10), rand(1:10, 2, 10), rand(1:10, 2, 10), + fill(0.99, 2, 10), + fill(4, 2, 10),fill(0.1, 2, 10), fill(0.5, 2, 10), + fill(0.9, 2, 10), fill(1.0, 2, 10)) + tz = tz"UTC" timestamps = ZonedDateTime(2020, 1, 1, 0, tz):Hour(1):ZonedDateTime(2020,1,1,9, tz) attrs = Dict("type" => "Single-Region System") single_reg_sys_wo_attrs = SystemModel( - generators, storages, generatorstorages, timestamps, rand(1:20, 10)) + generators, storages, generatorstorages, demandresponses, timestamps, rand(1:20, 10)) single_reg_sys_with_attrs = SystemModel( generators, storages, generatorstorages, timestamps, rand(1:20, 10), attrs) @@ -56,6 +63,7 @@ gen_regions = [1:1, 2:2] stor_regions = [1:0, 1:2] genstor_regions = [1:1, 2:1] + dr_regions = [1:0, 1:2] line_interfaces = [1:2] # Multi-region constructor @@ -63,6 +71,7 @@ regions, interfaces, generators, gen_regions, storages, stor_regions, generatorstorages, genstor_regions, + demandresponses, dr_regions, lines, line_interfaces, timestamps) diff --git a/PRASCore.jl/test/Systems/assets.jl b/PRASCore.jl/test/Systems/assets.jl index 46a256f4..33b5df29 100644 --- a/PRASCore.jl/test/Systems/assets.jl +++ b/PRASCore.jl/test/Systems/assets.jl @@ -131,6 +131,25 @@ end end + @testset "DemandResponses" begin + + DemandResponses{10,1,Hour,MW,MWh}( + names, categories, vals_int, vals_int, vals_int, + vals_float, vals_int, vals_float, vals_float, vals_float, vals_float) + + @test_throws AssertionError DemandResponses{5,1,Hour,MW,MWh}( + names, categories, vals_int, vals_int, vals_int, + vals_float, vals_int, vals_float, vals_float, vals_float, vals_float,) + + @test_throws AssertionError DemandResponses{10,1,Hour,MW,MWh}( + names, categories[1:2], vals_int, vals_int, vals_int, + vals_float, vals_int, vals_float, vals_float,vals_float, vals_float) + + @test_throws AssertionError DemandResponses{10,1,Hour,MW,MWh}( + names[1:2], categories[1:2], vals_int, vals_int, vals_int, + vals_float, vals_int, vals_float, vals_float,vals_float, vals_float) + end + @testset "Lines" begin lines = Lines{10,1,Hour,MW}( diff --git a/PRASFiles.jl/src/PRASFiles.jl b/PRASFiles.jl/src/PRASFiles.jl index 0106e2e1..1baed2da 100644 --- a/PRASFiles.jl/src/PRASFiles.jl +++ b/PRASFiles.jl/src/PRASFiles.jl @@ -1,8 +1,8 @@ module PRASFiles import PRASCore.Systems: SystemModel, Regions, Interfaces, - Generators, Storages, GeneratorStorages, Lines, - timeunits, powerunits, energyunits, unitsymbol + Generators, Storages, GeneratorStorages, DemandResponses, Lines, + timeunits, powerunits, energyunits, unitsymbol import PRASCore.Results: EUE, LOLE, NEUE, ShortfallResult, ShortfallSamplesResult, AbstractShortfallResult, Result import StatsBase: mean diff --git a/PRASFiles.jl/src/Systems/read.jl b/PRASFiles.jl/src/Systems/read.jl index fd37b210..d49f26cc 100644 --- a/PRASFiles.jl/src/Systems/read.jl +++ b/PRASFiles.jl/src/Systems/read.jl @@ -50,6 +50,8 @@ function _systemmodel_core(f::File) P = powerunits[read(metadata["power_unit"])] E = energyunits[read(metadata["energy_unit"])] + type_params = (N,L,T,P,E) + timestep = T(L) end_timestamp = start_timestamp + (N-1)*timestep timestamps = StepRange(start_timestamp, timestep, end_timestamp) @@ -98,9 +100,7 @@ function _systemmodel_core(f::File) else - generators = Generators{N,L,T,P}( - String[], String[], zeros(Int, 0, N), - zeros(Float64, 0, N), zeros(Float64, 0, N)) + generators = Generators{N,L,T,P}() region_gen_idxs = fill(1:0, n_regions) @@ -131,11 +131,7 @@ function _systemmodel_core(f::File) else - storages = Storages{N,L,T,P,E}( - String[], String[], - zeros(Int, 0, N), zeros(Int, 0, N), zeros(Int, 0, N), - zeros(Float64, 0, N), zeros(Float64, 0, N), zeros(Float64, 0, N), - zeros(Float64, 0, N), zeros(Float64, 0, N)) + storages = Storages{N,L,T,P,E}() region_stor_idxs = fill(1:0, n_regions) @@ -170,12 +166,7 @@ function _systemmodel_core(f::File) else - generatorstorages = GeneratorStorages{N,L,T,P,E}( - String[], String[], - zeros(Int, 0, N), zeros(Int, 0, N), zeros(Int, 0, N), - zeros(Float64, 0, N), zeros(Float64, 0, N), zeros(Float64, 0, N), - zeros(Int, 0, N), zeros(Int, 0, N), zeros(Int, 0, N), - zeros(Float64, 0, N), zeros(Float64, 0, N)) + generatorstorages = GeneratorStorages{N,L,T,P,E}() region_genstor_idxs = fill(1:0, n_regions) @@ -257,12 +248,9 @@ function _systemmodel_core(f::File) else - interfaces = Interfaces{N,P}( - Int[], Int[], zeros(Int, 0, N), zeros(Int, 0, N)) + interfaces = Interfaces{N,P}() - lines = Lines{N,L,T,P}( - String[], String[], zeros(Int, 0, N), zeros(Int, 0, N), - zeros(Float64, 0, N), zeros(Float64, 0, N)) + lines = Lines{N,L,T,P}() interface_line_idxs = UnitRange{Int}[] @@ -273,7 +261,7 @@ function _systemmodel_core(f::File) storages, region_stor_idxs, generatorstorages, region_genstor_idxs, lines, interface_line_idxs, - timestamps) + timestamps),type_params end """ @@ -281,7 +269,9 @@ Read a SystemModel from a PRAS file in version 0.5.x - 0.7.x format. """ function systemmodel_0_5(f::File) - return SystemModel(_systemmodel_core(f)...) + systemmodel_0_5_objs, _ = _systemmodel_core(f) + + return SystemModel(systemmodel_0_5_objs...) end @@ -290,12 +280,47 @@ Read a SystemModel from a PRAS file in version 0.8.x format. """ function systemmodel_0_8(f::File) - regions, interfaces, + (regions, interfaces, generators, region_gen_idxs, storages, region_stor_idxs, generatorstorages, region_genstor_idxs, lines, interface_line_idxs, - timestamps = _systemmodel_core(f) + timestamps), (N,L,T,P,E) = _systemmodel_core(f) + + n_regions = length(regions) + regionlookup = Dict(n=>i for (i, n) in enumerate(regions.names)) + + if haskey(f, "demandresponses") + + + dr_core = read(f["demandresponses/_core"]) + dr_names, dr_categories, dr_regionnames = readvector.( + Ref(dr_core), [:name, :category, :region]) + + dr_regions = getindex.(Ref(regionlookup), dr_regionnames) + region_order = sortperm(dr_regions) + + demandresponses = DemandResponses{N,L,T,P,E}( + dr_names[region_order], dr_categories[region_order], + load_matrix(f["demandresponses/borrowcapacity"], region_order, Int), + load_matrix(f["demandresponses/paybackcapacity"], region_order, Int), + load_matrix(f["demandresponses/energycapacity"], region_order, Int), + load_matrix(f["demandresponses/borrowedenergyinterest"], region_order, Float64), + load_matrix(f["demandresponses/allowablepaybackperiod"], region_order, Int), + load_matrix(f["demandresponses/failureprobability"], region_order, Float64), + load_matrix(f["demandresponses/repairprobability"], region_order, Float64), + load_matrix(f["demandresponses/borrowefficiency"], region_order, Float64), + load_matrix(f["demandresponses/paybackefficiency"], region_order, Float64), + ) + + region_dr_idxs = makeidxlist(dr_regions[region_order], n_regions) + + else + demandresponses = DemandResponses{N,L,T,P,E}() + + region_dr_idxs = fill(1:0, n_regions) + + end attrs = read_attrs(f) @@ -304,6 +329,7 @@ function systemmodel_0_8(f::File) generators, region_gen_idxs, storages, region_stor_idxs, generatorstorages, region_genstor_idxs, + demandresponses, region_dr_idxs, lines, interface_line_idxs, timestamps, attrs) diff --git a/PRASFiles.jl/src/Systems/write.jl b/PRASFiles.jl/src/Systems/write.jl index 4bda464f..1b3cb8c0 100644 --- a/PRASFiles.jl/src/Systems/write.jl +++ b/PRASFiles.jl/src/Systems/write.jl @@ -39,6 +39,14 @@ function savemodel( process_generatorstorages!(f, sys, string_length, compression_level) end + verbose && @info "$(length(sys.demandresponses)) DemandResponses " * + "found in the SystemModel." + + if length(sys.demandresponses) > 0 + verbose && @info "Processing DemandResponses for .pras file ..." + process_demandresponses!(f, sys, string_length, compression_level) + end + if length(sys.regions) > 1 verbose && @info "Processing Lines and Interfaces for .pras file ..." process_lines_interfaces!(f, sys, string_length, compression_level) @@ -213,6 +221,51 @@ function process_generatorstorages!( end +function process_demandresponses!( + f::File, sys::SystemModel, strlen::Int, compression::Int) + + demandresponses = create_group(f, "demandresponses") + + drs_core = Matrix{String}(undef, length(sys.demandresponses), 3) + drs_core_colnames = ["name", "category", "region"] + + drs_core[:, 1] = sys.demandresponses.names + drs_core[:, 2] = sys.demandresponses.categories + drs_core[:, 3] = regionnames( + length(sys.demandresponses), sys.regions.names, sys.region_dr_idxs) + + string_table!(demandresponses, "_core", drs_core_colnames, drs_core, strlen) + + demandresponses["borrowcapacity", deflate = compression] = + sys.demandresponses.borrow_capacity + + demandresponses["paybackcapacity", deflate = compression] = + sys.demandresponses.payback_capacity + + demandresponses["energycapacity", deflate = compression] = + sys.demandresponses.energy_capacity + + demandresponses["borrowefficiency", deflate = compression] = + sys.demandresponses.borrow_efficiency + + demandresponses["paybackefficiency", deflate = compression] = + sys.demandresponses.payback_efficiency + + demandresponses["borrowedenergyinterest", deflate = compression] = + sys.demandresponses.borrowed_energy_interest + + demandresponses["allowablepaybackperiod", deflate = compression] = + sys.demandresponses.allowable_payback_period + + demandresponses["failureprobability", deflate = compression] = sys.demandresponses.λ + + demandresponses["repairprobability", deflate = compression] = sys.demandresponses.μ + + return + +end + + function process_lines_interfaces!( f::File, sys::SystemModel, strlen::Int, compression::Int) diff --git a/docs/Project.toml b/docs/Project.toml index 4f3c2baa..c39c6e12 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -3,6 +3,7 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589" +Measures = "442fdcdd-2543-5da2-b0f3-8c86c306513e" PRASCapacityCredits = "2e1a2ed5-e89d-4cd3-bc86-c0e88a73d3a3" PRASCore = "c5c32b99-e7c3-4530-a685-6f76e19f7fe2" PRASFiles = "a2806276-6d43-4ef5-91c0-491704cd7cf1" diff --git a/docs/make.jl b/docs/make.jl index 801de9cd..85beb1c9 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -18,7 +18,6 @@ for file in files documenter = true, credit = true) end -examples_nav = fix_suffix.("./examples/" .* files) # Generate the unified documentation makedocs( @@ -33,7 +32,7 @@ makedocs( "Resource Adequacy" => "resourceadequacy.md", "Getting Started" => [ "Installation" => "installation.md", - "Quick start" => "quickstart.md", + "Quick Start" => "quickstart.md", ], "PRAS Components " => [ "System Model Specification" => "PRAS/sysmodelspec.md", @@ -42,7 +41,10 @@ makedocs( "Capacity Credit Calculation" => "PRAS/capacitycredit.md", ], ".pras File Format" => "SystemModel_HDF5_spec.md", - "Tutorials" => examples_nav, + "Tutorials" => [ + "PRAS 101 Walkthrough" => "examples/pras_walkthrough.md", + "Demand Response Walkthrough" => "examples/demand_response_walkthrough.md", + ], "Extending PRAS" => "extending.md", # "Contributing" => "contributing.md", "Changelog" => "changelog.md", diff --git a/docs/src/PRAS/sysmodelspec.md b/docs/src/PRAS/sysmodelspec.md index f92d1a8f..3bcafb77 100644 --- a/docs/src/PRAS/sysmodelspec.md +++ b/docs/src/PRAS/sysmodelspec.md @@ -44,33 +44,37 @@ methods -- while this is also possible with a single multi-year run, it requires some additional post-processing work. PRAS represents a power system as one or more **regions**, each containing -zero or more **generators**, **storages**, and -**generator-storages**. **Interfaces** contain **lines** and +zero or more **generators**, **storages**,**generator-storages**, and +**demand responses**. **Interfaces** contain **lines** and allow power transfer between two regions. The table below summarizes the characteristics of the different resource types (generators, storages, generator-storages, and lines), and the remainder of this section provides more details about each resource type and their associated resource collections (regions or interfaces). -| Parameter | Generator | Storage | Generator-Storage | Line | -|-----------|-----------|---------|-------------------|------| -| *Associated with a(n)...* | *Region* | *Region* | *Region* | *Interface* | -| *Name* | • | • | • | • | -| *Category* | • | • | • | • | +| Parameter | Generator | Storage | Generator-Storage | Demand Response | Line | +|-----------|-----------|---------|-------------------|-----------------|------| +| *Associated with a(n)...* | *Region* | *Region* | *Region* | *Region* | *Interface* | +| *Name* | • | • | • | • | • | +| *Category* | • | • | • | • | • | | Generation Capacity | • | | | | | Inflow Capacity | | | • | | | Charge Capacity | | • | • | | | Discharge Capacity | | • | • | | +| Load Borrow Capacity | | | | • | | +| Load Payback Capacity | | | | • | | | Energy Capacity | | • | • | | +| Borrowed Load Capacity | | | | • | | | Charge Efficiency | | • | • | | | Discharge Efficiency | | • | • | | | Carryover Efficiency | | • | • | | +| Borrowed Energy Interest | | | | • | | | Grid Injection Capacity | | | • | | | Grid Withdrawal Capacity | | | • | | -| Forward Transfer Capacity | | | | • | -| Backward Transfer Capacity | | | | • | -| Available→Unavailable Transition Probability | • | • | • | • | -| Unavailable→Available Transition Probability | • | • | • | • | +| Forward Transfer Capacity | | | | | • | +| Backward Transfer Capacity | | | | | • | +| Available→Unavailable Transition Probability | • | • | • | • | • | +| Unavailable→Available Transition Probability | • | • | • | • | • | *Table: PRAS resource parameters. Parameters in italic are fixed values; all others are provided as a time series.* @@ -113,8 +117,8 @@ time series, and so may be different during different time periods. ```@raw html
-Relations between power and energy parameters for generator, storage, and generator-storage resources. -
Relations between power and energy parameters for generator, storage, and generator-storage resources
+Relations between power and energy parameters for generator, storage, and generator-storage resources. +
Relations between power and energy parameters for generator, storage, generator-storage, and demand response resources
``` @@ -157,7 +161,7 @@ Just as with generators, storages may be in available or unavailable states, and move between these states randomly over time, according to provided state transition probabilities. Unavailable storages cannot inject power into or withdraw power from the grid, but they do maintain their energy state of charge -during an outage (minus any carryover losses occuring over time). +during an outage (minus any carryover losses occurring over time). ## Generator-Storages @@ -190,6 +194,61 @@ its state of charge during outages (subject to carryover losses). ``` +## Demand Responses + +Resources that can shift or shed electrical load forward in time but do +not provide an overall net addition of energy into the system (e.g., dr event programs) +are referred to as **demand responses** in PRAS. Like storages, demand response components are +associated with descriptive name and category metadata. +Each demand response unit has both a load borrow and load payback capacity time series, representing the +device's maximum ability to borrow load from or payback load into the grid +at a given point in time (as with storage capacity, these values may remain +constant over the simulation or may vary to reflect external constraints). + +Demand response units also have a maximum load capacity time series, reflecting the +maximum amount of borrowed load the device can hold at a given point in +time (increasing or decreasing this value will change the duration of time for +which the device could borrow or payback at maximum power). The demand response's +state of load increases with borrowing and decreases with +payback, and must always remain between zero and the maximum load +capacity in that time period. The energy flow relationships between +these capacities are depicted visually in the energy relation diagram. +If the maximum load capacity is less than the previous hour's state, the load borrowed will +be counted as unserved energy. + + +| **Demand Response Type** | **Description** | **Relevance to PRAS?** | +|:----------|:---------------|:----------------------| +| **Shift** | Demand response that encourages the movement of energy consumption from times of high demand to times of day when there is a surplus of generation.​ | Yes: can be dispatched as part of the PRAS optimization, to move load out of (net) peak periods.​ | +| **Shed** | Loads that can be curtailed to provide peak capacity reduction and support the system in emergency or contingency events with a range in advance notice times.​ | Yes: but caution must be used to avoid too frequent load shed and to capture only load that can truly be shed (i.e., not able to be paid back at a later time). Setting the `borrowed_energy_interest` to -0.99 will drop all effective borrowed load.​ | +| **Shape** | Demand response that reshapes customer load profiles through price response or behavioral campaigns — ‘load-modifying DR’ — with advance notice of months to days.​ | Only through pre-processing to load; no original implementation option within the PRAS optimization.| +| **Shimmy** | Loads that dynamically adjust to alleviate short-run ramps and disturbances on the system at timescales ranging from seconds up to an hour.​ | No: the system conditions and ramping are more granular than what is possible in PRAS. | +*Table: Common demand response categories and what is possible to model in PRAS. Demand Response type categories are taken from Piette, Mary Ann, et al. "2025 California Demand Response Potential Study.".* + + +Demand response units may incur losses or gains forward in time (borrowed energy interest). +The available load borrowed in the next time period is calculated by multiplying the borrowed energy interest +by the current load borrowed and adding to the current load borrowed. +The borrowed energy interest may be positive or negative, +indicating a growing or shrinking respectively borrowed load hour to hour. Borrowed energy interest +may lead to cases where the maximum load capacity of the demand response device is passed. If this occurs, +any load above the maximum capacity will be tracked as unserved energy. + + +Just as with storage, demand responses may be in available or unavailable states, and +move between these states randomly over time, according to provided state +transition probabilities. Unavailable demand responses cannot pay back load into or +withdraw load from the grid, but they do maintain their current borrowed load +during an outage (plus or minus any borrowed energy interest occurring over time). + +Demand response devices also have a maximum amount of time that they can hold borrowed energy. +This cutoff, where borrowed load is unable to be repaid and transitions over to unserved_energy, +is refereed to as the allowable payback period. This parameter can be time varying, and therefore +enable unique tailoring to the real world device being modeled. The allowable payback window is a integer +and follows the timestep units set for the system. If any surplus exists in the region, the demand response +device will attempt to payback any borrowed load, before charging storage. If the demand response device +pays back all borrowed load before the end of the period, the counter is reset upon further use. + ## Interfaces **Interfaces** define a potential capability to directly exchange power diff --git a/docs/src/PRASCore/api.md b/docs/src/PRASCore/api.md index 4cb5426c..a555af07 100644 --- a/docs/src/PRASCore/api.md +++ b/docs/src/PRASCore/api.md @@ -7,6 +7,7 @@ PRASCore.Systems.Regions PRASCore.Systems.Generators PRASCore.Systems.Storages PRASCore.Systems.GeneratorStorages +PRASCore.Systems.DemandResponses PRASCore.Systems.Lines PRASCore.Systems.Interfaces ``` @@ -37,5 +38,8 @@ PRASCore.Results.GeneratorStorageEnergySamples PRASCore.Results.StorageAvailability PRASCore.Results.StorageEnergy PRASCore.Results.StorageEnergySamples +PRASCore.Results.DemandResponseAvailability +PRASCore.Results.DemandResponseEnergy +PRASCore.Results.DemandResponseEnergySamples PRASCore.Results.LineAvailability ``` diff --git a/docs/src/changelog.md b/docs/src/changelog.md index acc22738..9282e19f 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -1,9 +1,11 @@ # Changelog -## [Unreleased] -- Initial changelog created. +## [0.8.0], 2025 - October -## [0.8.0], 2025 - September +- Add a demand response component which can model shift and shed type DR devices +- Add SystemModel attributes +- Enable pretty printing of different PRAS component information +- Enable auto-generated documentation website with detailed tutorials and walk-throughs ## [0.7], 2024 - December diff --git a/docs/src/images/resourceparameters.png b/docs/src/images/resourceparameters.png new file mode 100644 index 00000000..84fb9775 Binary files /dev/null and b/docs/src/images/resourceparameters.png differ