PRAS walkthrough
This is a complete example of running a PRAS assessment, using the RTS-GMLC system
Load the PRAS package and other tools necessary for analyses
using PRAS
using Plots
using DataFrames
using PrintfRead and explore a SystemModel
You can load in a system model from a .pras file if you have one like so:
sys = SystemModel("mysystem.pras")For the purposes of this example, we'll just use the built-in RTS-GMLC model.
sys = PRAS.rts_gmlc();We see some information about the system by just typing its name (or rather the variable that holds it):
sysPRAS system with 3 regions, and 3 interfaces between these regions.
Region names: 1, 2, 3
Assets:
Generators: 153 units
Storages: 1 units
GeneratorStorages: 1 units
Lines: 6
Time series:
Start time: 2020-01-01T00:00:00+00:00
Resolution: 1 Hour
Number of time steps: 8784
End time: 2020-12-31T23:00:00+00:00
Time zone: UTC
We retrieve the parameters of the system using the get_params function and use this for the plots below to ensure we have correct units:
(timesteps,periodlen,periodunit,powerunit,energyunit) = get_params(rts_gmlc())(8784, 1, Hour, MW, MWh)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:
sys["2"]Region: 2 (index - 2)
Peak load: 2850 MW
Generators: 35 units [indices - 52:86]
Storages: 0 units [indices - 1:0]
GeneratorStorages: 1 units [indices - 1:1]We can visualize a time series of the regional load in region "2":
region_2_load = sys.regions.load[sys["2"].index,:]
plot(sys.timestamps, region_2_load,
xlabel="Time", ylabel="Region 2 load ($(powerunit))",
legend=false)
We can find more information about all the Generators in the system by retriving the generators in the SystemModel:
system_generators = sys.generators153 Generators:
Category | Count
--------------------
Gas CC | 10
Gas CT | 27
Oil CT | 12
Coal | 16
Wind | 4
Hydro | 20
Nuclear | 1
Oil ST | 7
Solar PV | 25
Solar RTPV | 31
This returns an object of the asset type 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):
system_generators.capacity153×8784 Matrix{Int64}:
76 76 76 76 76 76 76 76 … 76 76 76 76 76 76 76
76 76 76 76 76 76 76 76 76 76 76 76 76 76 76
76 76 76 76 76 76 76 76 76 76 76 76 76 76 76
76 76 76 76 76 76 76 76 76 76 76 76 76 76 76
155 155 155 155 155 155 155 155 155 155 155 155 155 155 155
155 155 155 155 155 155 155 155 … 155 155 155 155 155 155 155
155 155 155 155 155 155 155 155 155 155 155 155 155 155 155
350 350 350 350 350 350 350 350 350 350 350 350 350 350 350
355 355 355 355 355 355 355 355 355 355 355 355 355 355 355
355 355 355 355 355 355 355 355 355 355 355 355 355 355 355
⋮ ⋮ ⋱ ⋮
0 0 0 0 0 0 0 4 0 0 0 0 0 0 0
0 0 0 0 0 0 0 4 … 0 0 0 0 0 0 0
0 0 0 0 0 0 0 4 0 0 0 0 0 0 0
0 0 0 0 0 0 0 4 0 0 0 0 0 0 0
0 0 0 0 0 0 0 5 0 0 0 0 0 0 0
0 0 0 0 0 0 0 2 0 0 0 0 0 0 0
481 635 487 433 408 440 377 199 … 222 281 203 167 221 251 220
143 139 145 145 137 99 62 47 0 0 0 0 0 0 0
795 794 774 767 752 719 655 595 11 33 74 90 60 21 16We can visualize a time series of the total system capacity (sum over individual generators' capacity at each time step)
plot(sys.timestamps, sum(system_generators.capacity, dims=1)',
xlabel="Time", ylabel="Total system capacity (MW)", legend=false)
Or, by category of generators:
category_indices = Dict([cat => findall(==(cat), system_generators.categories)
for cat in unique(system_generators.categories)]);
capacity_matrix = Vector{Vector{Int}}();
for (category,indices) in category_indices
push!(capacity_matrix, sum(system_generators.capacity[indices, :], dims=1)[1,:])
end
areaplot(sys.timestamps, hcat(capacity_matrix...),
label=permutedims(collect(keys(category_indices))),
xlabel="Time", ylabel="Total system capacity (MW)")
Similarly we can also retrieve all the Storages in the system and GenerationStorages in the system using sys.storages and sys.generatorstorages, respectively.
To retrieve the assets in a particular region, we can index by the region name and asset type (Generators here):
region_2_generators = sys["2", Generators]35 Generators:
Category | Count
--------------------
Gas CC | 3
Gas CT | 9
Oil CT | 4
Coal | 7
Hydro | 10
Solar RTPV | 1
Solar PV | 1
We get the storage device in region "3" like so:
region_3_storage = sys["3", Storages]1 Storages:
Category | Count
--------------------
Storage | 1
and the generation-storage device in region "2" like so:
region_2_genstorage = sys["2", GeneratorStorages]1 GeneratorStorages:
Category | Count
--------------------
CSP | 1
Run a Sequential Monte Carlo simulation
We can run a Sequential Monte Carlo simulation on this system using the assess function. Here we will also use four different result specifications:
shortfall, surplus, utilization, storage = assess(
sys, SequentialMonteCarlo(samples=100, seed=1),
Shortfall(), Surplus(), Utilization(), StorageEnergy());Start by checking the overall system adequacy:
lole = LOLE(shortfall); # event-hours per year
eue = EUE(shortfall); # unserved energy per year
println("System $(lole), $(eue)")System LOLE = 0.00000 event-h/8784h, EUE = 0.00000 MWh/8784hGiven we use only 100 samples and the RTS-GMLC system is quite reliable, we see a system which is reliable, with LOLE and EUE both near zero. For the purposes of this example, let's increase the system load homogenously by 700MW in every hour and region, and re-run the assessment.
sys.regions.load .+= 700.0
shortfall, surplus, utilization, storage = assess(
sys, SequentialMonteCarlo(samples=100, seed=1),
Shortfall(), Surplus(), Utilization(), StorageEnergy());
lole = LOLE(shortfall); # event-hours per year
eue = EUE(shortfall); # unserved energy per year
neue = NEUE(shortfall); # unserved energy per year
println("System $(lole), $(eue), $(neue)")System LOLE = 85±2 event-h/8784h, EUE = 26400±1100 MWh/8784h, NEUE = 470±20 ppmNow we see a system which is slightly unreliable with a normalized expected unserved energy (NEUE) of close to 470 parts per million of total load.
We can now look at the hourly loss-of-load expectation (LOLE) to see when when shortfalls are occurring. LOLE.(shortfall, many_hours) is Julia shorthand for calling LOLE on every timestep in the collection many_hours
lolps = LOLE.(shortfall, sys.timestamps)8784-element Vector{LOLE{1, 1, Hour}}:
LOLE = 0.00000 event-h/h
LOLE = 0.00000 event-h/h
LOLE = 0.00000 event-h/h
LOLE = 0.00000 event-h/h
LOLE = 0.00000 event-h/h
LOLE = 0.00000 event-h/h
LOLE = 0.00000 event-h/h
LOLE = 0.00000 event-h/h
LOLE = 0.00000 event-h/h
LOLE = 0.00000 event-h/h
⋮
LOLE = 0.00000 event-h/h
LOLE = 0.00000 event-h/h
LOLE = 0.00000 event-h/h
LOLE = 0.00000 event-h/h
LOLE = 0.00000 event-h/h
LOLE = 0.00000 event-h/h
LOLE = 0.00000 event-h/h
LOLE = 0.00000 event-h/h
LOLE = 0.00000 event-h/hHere results are in terms of event-hours per hour, which is equivalent to the loss-of-load probability (LOLP) for each hour. The LOLE object is shown as mean ± standard error. We are mostly interested in the mean here, we can retrieve this using val.(lolps) and visualize this:
plot(sys.timestamps, val.(lolps),
xlabel="Time", ylabel="Hourly LOLE (event-hours/hour)", legend=false)
We see the shortfall is concentrated in a few hours and there are many hours with LOLE = 1, which means that hour had a shortfall in every Monte Carlo sample.
We can find the regional NEUE for the entire simulation period, and obtain it in as a DataFrame for easier viewing:
regional_eue = DataFrame([(Region=reg_name, NEUE=val(NEUE(shortfall, reg_name)))
for reg_name in sys.regions.names],
)| Row | Region | NEUE |
|---|---|---|
| String | Float64 | |
| 1 | 1 | 684.134 |
| 2 | 2 | 569.856 |
| 3 | 3 | 173.376 |
So, region "1" has the highest overall NEUE, and has a higher load normalized shortfall
We may be interested in the EUE in the hour with highest LOLE
max_lole_ts = sys.timestamps[findfirst(val.(lolps).==1)];
println("Hour with first LOLE of 1.0: ", max_lole_ts)Hour with first LOLE of 1.0: 2020-07-24T18:00:00+00:00And we can find the unserved energy by region in that hour:
unserved_by_region = EUE.(shortfall, sys.regions.names, max_lole_ts)3-element Vector{EUE{1, 1, Hour, MWh}}:
EUE = 190±20 MWh/h
EUE = 220±20 MWh/h
EUE = 70±10 MWh/hwhich returns a Vector of EUE values for each region.
Region 2 has highest EUE in that hour, and we can look at the utilization of interfaces into that region in that period:
utilization_str = join([@sprintf("Interface between regions 1 and 2 utilization: %0.2f",
utilization["1" => "2", max_lole_ts][1]),
@sprintf("Interface between regions 3 and 2 utilization: %0.2f",
utilization["3" => "2", max_lole_ts][1]),
@sprintf("Interface between regions 1 and 3 utilization: %0.2f",
utilization["1" => "3", max_lole_ts][1])], "\n");
println(utilization_str)Interface between regions 1 and 2 utilization: 0.00
Interface between regions 3 and 2 utilization: 0.13
Interface between regions 1 and 3 utilization: 0.01We see that the interfaces are not fully utilized, meaning there is no excess generation in the system that could be wheeled 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]))Surplus in
region 1: 0.00
region 2: 0.00
region 3: 0.00Is 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]))Storage energy T-1: 63.39
Storage energy T: 24.40It 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 dropped. The device's ability to mitigate the shortfall must then be limited only by its discharge capacity, so increasing the regions storage capacity by adding more storage devices may help mitigate some shortfall.
Note that if the event was less consistent, this analysis could also have been performed on the subset of samples in which the event was observed, using the ShortfallSamples, UtilizationSamples, and StorageEnergySamples result specifications instead.
This page was generated using Literate.jl.