diff --git a/.gitignore b/.gitignore
index c725c05e..9b0f6caa 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,6 +4,7 @@
 Manifest.toml
 *.DS_Store
 *.pras
+!PRASFiles.jl/test/versioned_toy/*
 PRASFiles.jl/test/PRAS_Results_Export/
 
 docs/build/
diff --git a/PRASCore.jl/src/Simulations/Simulations.jl b/PRASCore.jl/src/Simulations/Simulations.jl
index 99b70e96..b3ed302f 100644
--- a/PRASCore.jl/src/Simulations/Simulations.jl
+++ b/PRASCore.jl/src/Simulations/Simulations.jl
@@ -182,10 +182,24 @@ function initialize!(
         initialize_availability!(
             rng, state.lines_available, state.lines_nexttransition,
             system.lines, N)
+        if size(system.storages.energy_capacity, 1) > 0
+            state.stors_energy .= (system.storages.initial_soc .* system.storages.energy_capacity[:,1])
+        else
+            fill!(state.stors_energy, 0.0)
+        end
+
+        if size(system.generatorstorages.energy_capacity, 1) > 0
+            state.genstors_energy .= (system.generatorstorages.initial_soc .* system.generatorstorages.energy_capacity[:,1])
+        else
+            fill!(state.genstors_energy, 0.0)
+        end
+        
+        if size(system.demandresponses.energy_capacity, 1) > 0
+            state.drs_energy .= (system.demandresponses.initial_borrowed_load .* system.demandresponses.energy_capacity[:,1])
+        else
+            fill!(state.drs_energy, 0.0)
+        end
 
-        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
diff --git a/PRASCore.jl/src/Systems/assets.jl b/PRASCore.jl/src/Systems/assets.jl
index a794a473..ba1dc6ad 100644
--- a/PRASCore.jl/src/Systems/assets.jl
+++ b/PRASCore.jl/src/Systems/assets.jl
@@ -178,6 +178,8 @@ A struct representing storage devices in the system.
  - `μ` (repair probability): Probability the unit transitions from forced outage to
    operational during a given simulation timestep, for each storage unit in each
    timeperiod. Unitless.
+ - `initial_soc`: Optional keyword for initial state of charge as a fraction [0,1.0] of `energy_capacity` at first timestep for 
+   each storage unit at the beginning of the simulation. Default is zero.
 """
 struct Storages{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} <: AbstractAssets{N,L,T,P}
 
@@ -195,12 +197,15 @@ struct Storages{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} <: AbstractAssets{N,L,
     λ::Matrix{Float64}
     μ::Matrix{Float64}
 
+    initial_soc::Vector{Float64} # energy
+
     function Storages{N,L,T,P,E}(
         names::Vector{<:AbstractString}, categories::Vector{<:AbstractString},
         chargecapacity::Matrix{Int}, dischargecapacity::Matrix{Int},
         energycapacity::Matrix{Int}, chargeefficiency::Matrix{Float64},
         dischargeefficiency::Matrix{Float64}, carryoverefficiency::Matrix{Float64},
-        λ::Matrix{Float64}, μ::Matrix{Float64}
+        λ::Matrix{Float64}, μ::Matrix{Float64}; 
+        initial_soc::Vector{Float64} = zeros(Float64, length(names))
     ) where {N,L,T,P,E}
 
         n_stors = length(names)
@@ -226,10 +231,14 @@ struct Storages{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} <: AbstractAssets{N,L,
         @assert all(isfractional, λ)
         @assert all(isfractional, μ)
 
+        @assert length(initial_soc) == n_stors
+        @assert all(isfractional, initial_soc)
+
+        #order of this is tied to how struct order is defined
         new{N,L,T,P,E}(string.(names), string.(categories),
                        chargecapacity, dischargecapacity, energycapacity,
                        chargeefficiency, dischargeefficiency, carryoverefficiency,
-                       λ, μ)
+                       λ, μ, initial_soc)
 
     end
 
@@ -242,7 +251,8 @@ function Storages{N,L,T,P,E}() where {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))
+                zeros(Float64, 0, N), zeros(Float64, 0, N);
+                initial_soc =  zeros(Float64, 0))
 end
 
 Base.:(==)(x::T, y::T) where {T <: Storages} =
@@ -255,13 +265,15 @@ Base.:(==)(x::T, y::T) where {T <: Storages} =
     x.discharge_efficiency == y.discharge_efficiency &&
     x.carryover_efficiency == y.carryover_efficiency &&
     x.λ == y.λ &&
-    x.μ == y.μ
+    x.μ == y.μ &&
+    x.initial_soc == y.initial_soc
 
 Base.getindex(s::S, idxs::AbstractVector{Int}) where {S <: Storages} =
     S(s.names[idxs], s.categories[idxs],s.charge_capacity[idxs,:],
       s.discharge_capacity[idxs, :],s.energy_capacity[idxs, :],
       s.charge_efficiency[idxs, :], s.discharge_efficiency[idxs, :], 
-      s.carryover_efficiency[idxs, :],s.λ[idxs, :], s.μ[idxs, :])
+      s.carryover_efficiency[idxs, :],s.λ[idxs, :], s.μ[idxs, :];
+      initial_soc = s.initial_soc[idxs])
 
 function Base.vcat(stors::Storages{N,L,T,P,E}...) where {N, L, T, P, E}
 
@@ -281,6 +293,8 @@ function Base.vcat(stors::Storages{N,L,T,P,E}...) where {N, L, T, P, E}
     λ = Matrix{Float64}(undef, n_stors, N)
     μ = Matrix{Float64}(undef, n_stors, N)
 
+    initial_soc = Vector{Float64}(undef, n_stors)
+
     last_idx = 0
 
     for s in stors
@@ -302,12 +316,15 @@ function Base.vcat(stors::Storages{N,L,T,P,E}...) where {N, L, T, P, E}
         λ[rows, :] = s.λ
         μ[rows, :] = s.μ
 
+        initial_soc[rows] = s.initial_soc
+
         last_idx += n
 
     end
 
     return Storages{N,L,T,P,E}(names, categories, charge_capacity, discharge_capacity, energy_capacity, charge_efficiency, discharge_efficiency, 
-                               carryover_efficiency, λ, μ)
+                               carryover_efficiency, λ, μ;
+                               initial_soc = initial_soc)
 
 end
 
@@ -351,6 +368,8 @@ A struct representing generator-storage hybrid devices within a power system.
  - `μ` (repair probability): Probability the unit transitions from forced outage to
    operational during a given simulation timestep, for each generator-storage unit in each
    timeperiod. Unitless.
+ - `initial_soc`: Optional keyword for initial state of charge as a fraction [0.0, 1.0] of 
+   `energy_capacity` at the first timestep for each storage unit. Default is zero.
 """
 struct GeneratorStorages{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} <: AbstractAssets{N,L,T,P}
 
@@ -372,6 +391,8 @@ struct GeneratorStorages{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} <: AbstractAs
     λ::Matrix{Float64}
     μ::Matrix{Float64}
 
+    initial_soc::Vector{Float64} # energy
+
     function GeneratorStorages{N,L,T,P,E}(
         names::Vector{<:AbstractString}, categories::Vector{<:AbstractString},
         charge_capacity::Matrix{Int}, discharge_capacity::Matrix{Int},
@@ -380,7 +401,8 @@ struct GeneratorStorages{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} <: AbstractAs
         carryover_efficiency::Matrix{Float64},
         inflow::Matrix{Int},
         gridwithdrawal_capacity::Matrix{Int}, gridinjection_capacity::Matrix{Int},
-        λ::Matrix{Float64}, μ::Matrix{Float64}
+        λ::Matrix{Float64}, μ::Matrix{Float64}; 
+        initial_soc::Vector{Float64} = zeros(Float64, length(names)),
     ) where {N,L,T,P,E}
 
         n_stors = length(names)
@@ -416,12 +438,15 @@ struct GeneratorStorages{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} <: AbstractAs
         @assert all(isfractional, λ)
         @assert all(isfractional, μ)
 
+        @assert length(initial_soc) == n_stors
+        @assert all(isfractional, initial_soc)
+
         new{N,L,T,P,E}(
             string.(names), string.(categories),
             charge_capacity, discharge_capacity, energy_capacity,
             charge_efficiency, discharge_efficiency, carryover_efficiency,
             inflow, gridwithdrawal_capacity, gridinjection_capacity,
-            λ, μ)
+            λ, μ,initial_soc)
 
     end
 
@@ -435,7 +460,7 @@ function GeneratorStorages{N,L,T,P,E}() where {N,L,T,P,E}
                 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))
+                zeros(Float64, 0, N), zeros(Float64, 0, N); initial_soc =  zeros(Float64, 0))
     
 end
 
@@ -452,7 +477,8 @@ Base.:(==)(x::T, y::T) where {T <: GeneratorStorages} =
     x.gridwithdrawal_capacity == y.gridwithdrawal_capacity &&
     x.gridinjection_capacity == y.gridinjection_capacity &&
     x.λ == y.λ &&
-    x.μ == y.μ
+    x.μ == y.μ &&
+    x.initial_soc == y.initial_soc
 
 Base.getindex(g_s::G, idxs::AbstractVector{Int}) where {G <: GeneratorStorages} =
     G(g_s.names[idxs], g_s.categories[idxs], g_s.charge_capacity[idxs,:],
@@ -460,7 +486,7 @@ Base.getindex(g_s::G, idxs::AbstractVector{Int}) where {G <: GeneratorStorages}
       g_s.charge_efficiency[idxs, :], g_s.discharge_efficiency[idxs, :], 
       g_s.carryover_efficiency[idxs, :],g_s.inflow[idxs, :],
       g_s.gridwithdrawal_capacity[idxs, :],g_s.gridinjection_capacity[idxs, :],
-      g_s.λ[idxs, :], g_s.μ[idxs, :])
+      g_s.λ[idxs, :], g_s.μ[idxs, :]; initial_soc =  g_s.initial_soc[idxs])
 
 function Base.vcat(gen_stors::GeneratorStorages{N,L,T,P,E}...) where {N, L, T, P, E}
 
@@ -484,6 +510,8 @@ function Base.vcat(gen_stors::GeneratorStorages{N,L,T,P,E}...) where {N, L, T, P
     λ = Matrix{Float64}(undef, n_gen_stors, N)
     μ = Matrix{Float64}(undef, n_gen_stors, N)
 
+    initial_soc = Vector{Float64}(undef, n_gen_stors)
+
     last_idx = 0
 
     for g_s in gen_stors
@@ -509,12 +537,14 @@ function Base.vcat(gen_stors::GeneratorStorages{N,L,T,P,E}...) where {N, L, T, P
         λ[rows, :] = g_s.λ
         μ[rows, :] = g_s.μ
 
+        initial_soc[rows] = g_s.initial_soc
+
         last_idx += n
 
     end
 
     return GeneratorStorages{N,L,T,P,E}(names, categories, charge_capacity, discharge_capacity, energy_capacity, charge_efficiency, discharge_efficiency, 
-                               carryover_efficiency,inflow, gridwithdrawal_capacity, gridinjection_capacity, λ, μ)
+                               carryover_efficiency,inflow, gridwithdrawal_capacity, gridinjection_capacity, λ, μ; initial_soc =  initial_soc)
 
 end
 
@@ -551,6 +581,8 @@ A struct representing demand response devices in the system.
  - `μ` (repair probability): Probability the unit transitions from forced outage to
    operational during a given simulation timestep, for each storage unit in each
    timeperiod. Unitless.
+ - `initial_borrowed_load`: Optional initial state of borrowed load as a fraction [0,1.0] of `energy_capacity` at first timestep for 
+   each demand response unit at the beginning of the simulation. Default is zero.
 """
 struct DemandResponses{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} <: AbstractAssets{N,L,T,P}
 
@@ -568,6 +600,8 @@ struct DemandResponses{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} <: AbstractAsse
     λ::Matrix{Float64}
     μ::Matrix{Float64}
 
+    initial_borrowed_load::Vector{Float64} # energy
+
     borrow_efficiency::Matrix{Float64}
     payback_efficiency::Matrix{Float64}
 
@@ -576,8 +610,10 @@ struct DemandResponses{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} <: AbstractAsse
         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}
+        λ::Matrix{Float64}, μ::Matrix{Float64};
+        initial_borrowed_load::Vector{Float64} = zeros(Float64, length(names)),
+        borrow_efficiency::Matrix{Float64} = ones(Float64, size(borrowcapacity)),
+        payback_efficiency::Matrix{Float64} = ones(Float64, size(paybackcapacity)) 
     ) where {N,L,T,P,E}
 
         n_drs = length(names)
@@ -592,11 +628,11 @@ struct DemandResponses{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} <: AbstractAsse
         @assert all(isnonnegative, paybackcapacity)
         @assert all(isnonnegative, energycapacity)
 
-        @assert size(borrowefficiency) == (n_drs, N)
-        @assert size(paybackefficiency) == (n_drs, N)
+        @assert size(borrow_efficiency) == (n_drs, N)
+        @assert size(payback_efficiency) == (n_drs, N)
         @assert size(borrowedenergyinterest) == (n_drs, N)
-        @assert all(isfractional, borrowefficiency)
-        @assert all(isfractional, paybackefficiency)
+        @assert all(isfractional, borrow_efficiency)
+        @assert all(isfractional, payback_efficiency)
         @assert all(borrowedenergyinterest .<= 1.0)
         @assert all(borrowedenergyinterest .>= -1.0)
 
@@ -609,41 +645,31 @@ struct DemandResponses{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit} <: AbstractAsse
         @assert all(isfractional, λ)
         @assert all(isfractional, μ)
 
+        @assert length(initial_borrowed_load) == n_drs
+        @assert all(isfractional, initial_borrowed_load)
+
         new{N,L,T,P,E}(string.(names), string.(categories),
                        borrowcapacity, paybackcapacity, energycapacity,
                       borrowedenergyinterest,
                        allowablepaybackperiod,
-                       λ, μ,borrowefficiency, paybackefficiency,)
+                       λ, μ,
+                       initial_borrowed_load,
+                       borrow_efficiency,
+                       payback_efficiency)
     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))
+  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);
+              initial_borrowed_load = zeros(Float64, 0),
+              borrow_efficiency = Matrix{Float64}(undef, 0, N),
+              payback_efficiency = Matrix{Float64}(undef, 0, N))
 end
 
 Base.:(==)(x::T, y::T) where {T <: DemandResponses} =
@@ -657,13 +683,16 @@ Base.:(==)(x::T, y::T) where {T <: DemandResponses} =
     x.borrowed_energy_interest == y.borrowed_energy_interest &&
     x.allowable_payback_period == y.allowable_payback_period &&
     x.λ == y.λ &&
-    x.μ == y.μ
+    x.μ == y.μ &&
+    x.initial_borrowed_load == y.initial_borrowed_load
 
 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, :])
+      dr.borrowed_energy_interest[idxs, :],dr.allowable_payback_period[idxs, :],dr.λ[idxs, :], dr.μ[idxs, :];
+      initial_borrowed_load = dr.initial_borrowed_load[idxs],
+      borrow_efficiency = dr.borrow_efficiency[idxs, :], 
+      payback_efficiency = dr.payback_efficiency[idxs, :])
 
 function Base.vcat(drs::DemandResponses{N,L,T,P,E}...) where {N, L, T, P, E}
 
@@ -686,6 +715,8 @@ function Base.vcat(drs::DemandResponses{N,L,T,P,E}...) where {N, L, T, P, E}
     λ = Matrix{Float64}(undef, n_drs, N)
     μ = Matrix{Float64}(undef, n_drs, N)
 
+    initial_borrowed_load = Vector{Float64}(undef, n_drs)
+
     last_idx = 0
 
     for dr in drs
@@ -709,12 +740,17 @@ function Base.vcat(drs::DemandResponses{N,L,T,P,E}...) where {N, L, T, P, E}
         λ[rows, :] = dr.λ
         μ[rows, :] = dr.μ
 
+        initial_borrowed_load[rows] = dr.initial_borrowed_load
+
         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)
+    borrowed_energy_interest,allowable_payback_period, λ, μ;
+    initial_borrowed_load = initial_borrowed_load,
+    borrow_efficiency = borrow_efficiency, 
+    payback_efficiency = payback_efficiency)
 
 end
 
diff --git a/PRASCore.jl/test/Systems/SystemModel.jl b/PRASCore.jl/test/Systems/SystemModel.jl
index 0bf5f577..d3a8b80b 100644
--- a/PRASCore.jl/test/Systems/SystemModel.jl
+++ b/PRASCore.jl/test/Systems/SystemModel.jl
@@ -21,8 +21,10 @@
         ["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))
+        fill(4, 2, 10),fill(0.1, 2, 10), fill(0.5, 2, 10);
+        initial_borrowed_load = zeros(Float64, 2),
+        borrow_efficiency = fill(0.9, 2, 10),
+        payback_efficiency = fill(1.0, 2, 10))
 
     tz = tz"UTC"
     timestamps = ZonedDateTime(2020, 1, 1, 0, tz):Hour(1):ZonedDateTime(2020,1,1,9, tz)
diff --git a/PRASCore.jl/test/Systems/assets.jl b/PRASCore.jl/test/Systems/assets.jl
index 33b5df29..e788ffa6 100644
--- a/PRASCore.jl/test/Systems/assets.jl
+++ b/PRASCore.jl/test/Systems/assets.jl
@@ -5,7 +5,7 @@
 
     vals_int   = rand(1:100, 3, 10)
     vals_float = rand(3, 10)
-
+    vals_float_soc = rand(0.0:0.1:1.0, length(names))
     @testset "Generators" begin
 
         gens = Generators{10,1,Hour,MW}(
@@ -46,24 +46,24 @@
 
         stors = Storages{10,1,Hour,MW,MWh}(
             names, categories, vals_int, vals_int, vals_int,
-            vals_float, vals_float, vals_float, vals_float, vals_float)
+            vals_float, vals_float, vals_float, vals_float, vals_float;initial_soc = vals_float_soc)
         @test stors isa Storages
 
         @test_throws AssertionError Storages{5,1,Hour,MW,MWh}(
             names, categories, vals_int, vals_int, vals_int,
-            vals_float, vals_float, vals_float, vals_float, vals_float)
+            vals_float, vals_float, vals_float, vals_float, vals_float;initial_soc =  vals_float_soc)
 
         @test_throws AssertionError Storages{10,1,Hour,MW,MWh}(
             names, categories[1:2], vals_int, vals_int, vals_int,
-            vals_float, vals_float, vals_float, vals_float, vals_float)
+            vals_float, vals_float, vals_float, vals_float, vals_float;initial_soc =  vals_float_soc)
 
         @test_throws AssertionError Storages{10,1,Hour,MW,MWh}(
             names[1:2], categories[1:2], vals_int, vals_int, vals_int,
-            vals_float, vals_float, vals_float, vals_float, vals_float)
+            vals_float, vals_float, vals_float, vals_float, vals_float;initial_soc =  vals_float_soc)
 
         @test_throws AssertionError Storages{10,1,Hour,MW,MWh}(
             names, categories, vals_int, vals_int, vals_int,
-            vals_float, vals_float, -vals_float, vals_float, vals_float)
+            vals_float, vals_float, -vals_float, vals_float, vals_float;initial_soc =  vals_float_soc)
 
         @testset "Printing" begin
             io = IOBuffer()
@@ -89,29 +89,29 @@
         gen_stors = GeneratorStorages{10,1,Hour,MW,MWh}(
             names, categories,
             vals_int, vals_int, vals_int, vals_float, vals_float, vals_float,
-            vals_int, vals_int, vals_int, vals_float, vals_float)
+            vals_int, vals_int, vals_int, vals_float, vals_float;initial_soc =  vals_float_soc)
         @test gen_stors isa GeneratorStorages
 
         @test_throws AssertionError GeneratorStorages{5,1,Hour,MW,MWh}(
             names, categories,
             vals_int, vals_int, vals_int, vals_float, vals_float, vals_float,
-            vals_int, vals_int, vals_int, vals_float, vals_float)
+            vals_int, vals_int, vals_int, vals_float, vals_float;initial_soc =  vals_float_soc)
 
 
         @test_throws AssertionError GeneratorStorages{10,1,Hour,MW,MWh}(
             names, categories[1:2],
             vals_int, vals_int, vals_int, vals_float, vals_float, vals_float,
-            vals_int, vals_int, vals_int, vals_float, vals_float)
+            vals_int, vals_int, vals_int, vals_float, vals_float;initial_soc =  vals_float_soc)
 
         @test_throws AssertionError GeneratorStorages{10,1,Hour,MW,MWh}(
             names[1:2], categories[1:2],
             vals_int, vals_int, vals_int, vals_float, vals_float, vals_float,
-            vals_int, vals_int, vals_int, vals_float, vals_float)
+            vals_int, vals_int, vals_int, vals_float, vals_float;initial_soc =  vals_float_soc)
 
         @test_throws AssertionError GeneratorStorages{10,1,Hour,MW,MWh}(
             names, categories,
             vals_int, vals_int, vals_int, vals_float, vals_float, -vals_float,
-            vals_int, vals_int, vals_int, vals_float, vals_float)
+            vals_int, vals_int, vals_int, vals_float, vals_float;initial_soc =  vals_float_soc)
 
         @testset "Printing" begin
             io = IOBuffer()
@@ -135,19 +135,19 @@
 
         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)
+            vals_float, vals_int, vals_float, vals_float;initial_borrowed_load =  vals_float_soc, borrow_efficiency = vals_float, payback_efficiency = 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,)
+            vals_float, vals_int, vals_float, vals_float;initial_borrowed_load =  vals_float_soc, borrow_efficiency = vals_float, payback_efficiency = 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)
+            vals_float, vals_int, vals_float, vals_float;initial_borrowed_load =  vals_float_soc, borrow_efficiency = vals_float, payback_efficiency = 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)
+            vals_float, vals_int, vals_float, vals_float;initial_borrowed_load =  vals_float_soc, borrow_efficiency = vals_float, payback_efficiency = vals_float)
     end
 
     @testset "Lines" begin
diff --git a/PRASFiles.jl/src/Systems/read.jl b/PRASFiles.jl/src/Systems/read.jl
index d49f26cc..8b2af244 100644
--- a/PRASFiles.jl/src/Systems/read.jl
+++ b/PRASFiles.jl/src/Systems/read.jl
@@ -20,8 +20,10 @@ function SystemModel(inputfile::String)
         # Determine the appropriate version of the importer to use
         return if (0,5,0) <= version < (0,8,0)
             systemmodel_0_5(f)
-        elseif (0,8,0) <= version < (0,9,0)
-            systemmodel_0_8(f)
+        elseif version == (0,8,0)
+            systemmodel_0_8_0(f)
+        elseif version >= (0,8,1)
+            systemmodel_0_8_1(f)
         else
             error("PRAS file format $versionstring not supported by this version of PRASBase.")
         end
@@ -58,7 +60,6 @@ function _systemmodel_core(f::File)
 
     has_regions = haskey(f, "regions")
     has_generators = haskey(f, "generators")
-    has_storages = haskey(f, "storages")
     has_generatorstorages = haskey(f, "generatorstorages")
     has_interfaces = haskey(f, "interfaces")
     has_lines = haskey(f, "lines")
@@ -106,72 +107,6 @@ function _systemmodel_core(f::File)
 
     end
 
-    if has_storages
-
-        stor_core = read(f["storages/_core"])
-        stor_names, stor_categories, stor_regionnames = readvector.(
-            Ref(stor_core), [:name, :category, :region])
-
-        stor_regions = getindex.(Ref(regionlookup), stor_regionnames)
-        region_order = sortperm(stor_regions)
-
-        storages = Storages{N,L,T,P,E}(
-            stor_names[region_order], stor_categories[region_order],
-            load_matrix(f["storages/chargecapacity"], region_order, Int),
-            load_matrix(f["storages/dischargecapacity"], region_order, Int),
-            load_matrix(f["storages/energycapacity"], region_order, Int),
-            load_matrix(f["storages/chargeefficiency"], region_order, Float64),
-            load_matrix(f["storages/dischargeefficiency"], region_order, Float64),
-            load_matrix(f["storages/carryoverefficiency"], region_order, Float64),
-            load_matrix(f["storages/failureprobability"], region_order, Float64),
-            load_matrix(f["storages/repairprobability"], region_order, Float64)
-        )
-
-        region_stor_idxs = makeidxlist(stor_regions[region_order], n_regions)
-
-    else
-
-        storages = Storages{N,L,T,P,E}()
-
-        region_stor_idxs = fill(1:0, n_regions)
-
-    end
-
-
-    if has_generatorstorages
-
-        genstor_core = read(f["generatorstorages/_core"])
-        genstor_names, genstor_categories, genstor_regionnames = readvector.(
-            Ref(genstor_core), [:name, :category, :region])
-
-        genstor_regions = getindex.(Ref(regionlookup), genstor_regionnames)
-        region_order = sortperm(genstor_regions)
-
-        generatorstorages = GeneratorStorages{N,L,T,P,E}(
-            genstor_names[region_order], genstor_categories[region_order],
-            load_matrix(f["generatorstorages/chargecapacity"], region_order, Int),
-            load_matrix(f["generatorstorages/dischargecapacity"], region_order, Int),
-            load_matrix(f["generatorstorages/energycapacity"], region_order, Int),
-            load_matrix(f["generatorstorages/chargeefficiency"], region_order, Float64),
-            load_matrix(f["generatorstorages/dischargeefficiency"], region_order, Float64),
-            load_matrix(f["generatorstorages/carryoverefficiency"], region_order, Float64),
-            load_matrix(f["generatorstorages/inflow"], region_order, Int),
-            load_matrix(f["generatorstorages/gridwithdrawalcapacity"], region_order, Int),
-            load_matrix(f["generatorstorages/gridinjectioncapacity"], region_order, Int),
-            load_matrix(f["generatorstorages/failureprobability"], region_order, Float64),
-            load_matrix(f["generatorstorages/repairprobability"], region_order, Float64)
-        )
-
-        region_genstor_idxs = makeidxlist(genstor_regions[region_order], n_regions)
-
-    else
-
-        generatorstorages = GeneratorStorages{N,L,T,P,E}()
-
-        region_genstor_idxs = fill(1:0, n_regions)
-
-    end
-
     if has_interfaces
 
         interfaces_core = read(f["interfaces/_core"])
@@ -258,8 +193,6 @@ function _systemmodel_core(f::File)
 
     return (regions, interfaces,
             generators, region_gen_idxs,
-            storages, region_stor_idxs,
-            generatorstorages, region_genstor_idxs,
             lines, interface_line_idxs,
             timestamps),type_params
 end
@@ -269,29 +202,179 @@ Read a SystemModel from a PRAS file in version 0.5.x - 0.7.x format.
 """
 function systemmodel_0_5(f::File)
 
-    systemmodel_0_5_objs, _ = _systemmodel_core(f)
+    (regions, interfaces,
+    generators, region_gen_idxs,
+    lines, interface_line_idxs,
+    timestamps), (N,L,T,P,E) = _systemmodel_core(f)
 
-    return SystemModel(systemmodel_0_5_objs...)
+    n_regions = length(regions)
+    regionlookup = Dict(n=>i for (i, n) in enumerate(regions.names))
+    attrs = read_attrs(f)
+
+    has_storages = haskey(f, "storages")
+    has_generatorstorages = haskey(f, "generatorstorages")
+
+    if has_storages
+
+        stor_core = read(f["storages/_core"])
+        stor_names, stor_categories, stor_regionnames = readvector.(
+            Ref(stor_core), [:name, :category, :region])
+
+        stor_regions = getindex.(Ref(regionlookup), stor_regionnames)
+        region_order = sortperm(stor_regions)
+
+        storages = Storages{N,L,T,P,E}(
+            stor_names[region_order], stor_categories[region_order],
+            load_matrix(f["storages/chargecapacity"], region_order, Int),
+            load_matrix(f["storages/dischargecapacity"], region_order, Int),
+            load_matrix(f["storages/energycapacity"], region_order, Int),
+            load_matrix(f["storages/chargeefficiency"], region_order, Float64),
+            load_matrix(f["storages/dischargeefficiency"], region_order, Float64),
+            load_matrix(f["storages/carryoverefficiency"], region_order, Float64),
+            load_matrix(f["storages/failureprobability"], region_order, Float64),
+            load_matrix(f["storages/repairprobability"], region_order, Float64)
+        )
+
+        region_stor_idxs = makeidxlist(stor_regions[region_order], n_regions)
+
+    else
+
+        storages = Storages{N,L,T,P,E}()
+
+        region_stor_idxs = fill(1:0, n_regions)
+
+    end
+
+
+    if has_generatorstorages
+
+        genstor_core = read(f["generatorstorages/_core"])
+        genstor_names, genstor_categories, genstor_regionnames = readvector.(
+            Ref(genstor_core), [:name, :category, :region])
+
+        genstor_regions = getindex.(Ref(regionlookup), genstor_regionnames)
+        region_order = sortperm(genstor_regions)
+
+        generatorstorages = GeneratorStorages{N,L,T,P,E}(
+            genstor_names[region_order], genstor_categories[region_order],
+            load_matrix(f["generatorstorages/chargecapacity"], region_order, Int),
+            load_matrix(f["generatorstorages/dischargecapacity"], region_order, Int),
+            load_matrix(f["generatorstorages/energycapacity"], region_order, Int),
+            load_matrix(f["generatorstorages/chargeefficiency"], region_order, Float64),
+            load_matrix(f["generatorstorages/dischargeefficiency"], region_order, Float64),
+            load_matrix(f["generatorstorages/carryoverefficiency"], region_order, Float64),
+            load_matrix(f["generatorstorages/inflow"], region_order, Int),
+            load_matrix(f["generatorstorages/gridwithdrawalcapacity"], region_order, Int),
+            load_matrix(f["generatorstorages/gridinjectioncapacity"], region_order, Int),
+            load_matrix(f["generatorstorages/failureprobability"], region_order, Float64),
+            load_matrix(f["generatorstorages/repairprobability"], region_order, Float64)
+        )
+
+        region_genstor_idxs = makeidxlist(genstor_regions[region_order], n_regions)
+
+    else
+
+        generatorstorages = GeneratorStorages{N,L,T,P,E}()
+
+        region_genstor_idxs = fill(1:0, n_regions)
+
+    end
+
+    return SystemModel(
+        regions, interfaces,
+        generators, region_gen_idxs,
+        storages, region_stor_idxs,
+        generatorstorages, region_genstor_idxs,
+        lines, interface_line_idxs,
+        timestamps, attrs)
 
 end
 
 """
-Read a SystemModel from a PRAS file in version 0.8.x format.
+Read a SystemModel from a PRAS file in version 0.8.0 format. Adds Demand Response component info
 """
-function systemmodel_0_8(f::File)
+function systemmodel_0_8_0(f::File)
     
     (regions, interfaces,
     generators, region_gen_idxs,
-    storages, region_stor_idxs,
-    generatorstorages, region_genstor_idxs,
     lines, interface_line_idxs,
     timestamps), (N,L,T,P,E) = _systemmodel_core(f)
 
     n_regions = length(regions)
     regionlookup = Dict(n=>i for (i, n) in enumerate(regions.names))
+    attrs = read_attrs(f)
+
+    has_storages = haskey(f, "storages")
+    has_generatorstorages = haskey(f, "generatorstorages")
+    has_demandresponses = haskey(f, "demandresponses")
+
+    if has_storages
+
+        stor_core = read(f["storages/_core"])
+        stor_names, stor_categories, stor_regionnames = readvector.(
+            Ref(stor_core), [:name, :category, :region])
+
+        stor_regions = getindex.(Ref(regionlookup), stor_regionnames)
+        region_order = sortperm(stor_regions)
+
+        storages = Storages{N,L,T,P,E}(
+            stor_names[region_order], stor_categories[region_order],
+            load_matrix(f["storages/chargecapacity"], region_order, Int),
+            load_matrix(f["storages/dischargecapacity"], region_order, Int),
+            load_matrix(f["storages/energycapacity"], region_order, Int),
+            load_matrix(f["storages/chargeefficiency"], region_order, Float64),
+            load_matrix(f["storages/dischargeefficiency"], region_order, Float64),
+            load_matrix(f["storages/carryoverefficiency"], region_order, Float64),
+            load_matrix(f["storages/failureprobability"], region_order, Float64),
+            load_matrix(f["storages/repairprobability"], region_order, Float64)
+        )
+
+        region_stor_idxs = makeidxlist(stor_regions[region_order], n_regions)
+
+    else
+
+        storages = Storages{N,L,T,P,E}()
+
+        region_stor_idxs = fill(1:0, n_regions)
+
+    end
 
-    if haskey(f, "demandresponses")
 
+    if has_generatorstorages
+
+        genstor_core = read(f["generatorstorages/_core"])
+        genstor_names, genstor_categories, genstor_regionnames = readvector.(
+            Ref(genstor_core), [:name, :category, :region])
+
+        genstor_regions = getindex.(Ref(regionlookup), genstor_regionnames)
+        region_order = sortperm(genstor_regions)
+
+        generatorstorages = GeneratorStorages{N,L,T,P,E}(
+            genstor_names[region_order], genstor_categories[region_order],
+            load_matrix(f["generatorstorages/chargecapacity"], region_order, Int),
+            load_matrix(f["generatorstorages/dischargecapacity"], region_order, Int),
+            load_matrix(f["generatorstorages/energycapacity"], region_order, Int),
+            load_matrix(f["generatorstorages/chargeefficiency"], region_order, Float64),
+            load_matrix(f["generatorstorages/dischargeefficiency"], region_order, Float64),
+            load_matrix(f["generatorstorages/carryoverefficiency"], region_order, Float64),
+            load_matrix(f["generatorstorages/inflow"], region_order, Int),
+            load_matrix(f["generatorstorages/gridwithdrawalcapacity"], region_order, Int),
+            load_matrix(f["generatorstorages/gridinjectioncapacity"], region_order, Int),
+            load_matrix(f["generatorstorages/failureprobability"], region_order, Float64),
+            load_matrix(f["generatorstorages/repairprobability"], region_order, Float64)
+        )
+
+        region_genstor_idxs = makeidxlist(genstor_regions[region_order], n_regions)
+
+    else
+
+        generatorstorages = GeneratorStorages{N,L,T,P,E}()
+
+        region_genstor_idxs = fill(1:0, n_regions)
+
+    end
+
+    if has_demandresponses
 
         dr_core = read(f["demandresponses/_core"])
         dr_names, dr_categories, dr_regionnames = readvector.(
@@ -322,8 +405,135 @@ function systemmodel_0_8(f::File)
 
     end
 
+    return SystemModel(
+        regions, interfaces,
+        generators, region_gen_idxs,
+        storages, region_stor_idxs,
+        generatorstorages, region_genstor_idxs,
+        demandresponses, region_dr_idxs,
+        lines, interface_line_idxs,
+        timestamps, attrs)
+    
+end
+
+"""
+Read a SystemModel from a PRAS file in version 0.8.1+ format. Requires initial SOC levels for storages/generator-storages and initial borrowed load/efficiency parameters for demand responses.
+"""
+function systemmodel_0_8_1(f::File)
+    
+    (regions, interfaces,
+    generators, region_gen_idxs,
+    lines, interface_line_idxs,
+    timestamps), (N,L,T,P,E) = _systemmodel_core(f)
+
+    n_regions = length(regions)
+    regionlookup = Dict(n=>i for (i, n) in enumerate(regions.names))
     attrs = read_attrs(f)
 
+    has_storages = haskey(f, "storages")
+    has_generatorstorages = haskey(f, "generatorstorages")
+    has_demandresponses = haskey(f, "demandresponses")
+
+    if has_storages
+
+        stor_core = read(f["storages/_core"])
+        stor_names, stor_categories, stor_regionnames = readvector.(
+            Ref(stor_core), [:name, :category, :region])
+
+        stor_regions = getindex.(Ref(regionlookup), stor_regionnames)
+        region_order = sortperm(stor_regions)
+
+        storages = Storages{N,L,T,P,E}(
+            stor_names[region_order], stor_categories[region_order],
+            load_matrix(f["storages/chargecapacity"], region_order, Int),
+            load_matrix(f["storages/dischargecapacity"], region_order, Int),
+            load_matrix(f["storages/energycapacity"], region_order, Int),
+            load_matrix(f["storages/chargeefficiency"], region_order, Float64),
+            load_matrix(f["storages/dischargeefficiency"], region_order, Float64),
+            load_matrix(f["storages/carryoverefficiency"], region_order, Float64),
+            load_matrix(f["storages/failureprobability"], region_order, Float64),
+            load_matrix(f["storages/repairprobability"], region_order, Float64);
+            initial_soc = load_matrix(f["storages/initialsoc"], region_order, Float64)
+        )
+
+        region_stor_idxs = makeidxlist(stor_regions[region_order], n_regions)
+
+    else
+
+        storages = Storages{N,L,T,P,E}()
+
+        region_stor_idxs = fill(1:0, n_regions)
+
+    end
+
+
+    if has_generatorstorages
+
+        genstor_core = read(f["generatorstorages/_core"])
+        genstor_names, genstor_categories, genstor_regionnames = readvector.(
+            Ref(genstor_core), [:name, :category, :region])
+
+        genstor_regions = getindex.(Ref(regionlookup), genstor_regionnames)
+        region_order = sortperm(genstor_regions)
+
+        generatorstorages = GeneratorStorages{N,L,T,P,E}(
+            genstor_names[region_order], genstor_categories[region_order],
+            load_matrix(f["generatorstorages/chargecapacity"], region_order, Int),
+            load_matrix(f["generatorstorages/dischargecapacity"], region_order, Int),
+            load_matrix(f["generatorstorages/energycapacity"], region_order, Int),
+            load_matrix(f["generatorstorages/chargeefficiency"], region_order, Float64),
+            load_matrix(f["generatorstorages/dischargeefficiency"], region_order, Float64),
+            load_matrix(f["generatorstorages/carryoverefficiency"], region_order, Float64),
+            load_matrix(f["generatorstorages/inflow"], region_order, Int),
+            load_matrix(f["generatorstorages/gridwithdrawalcapacity"], region_order, Int),
+            load_matrix(f["generatorstorages/gridinjectioncapacity"], region_order, Int),
+            load_matrix(f["generatorstorages/failureprobability"], region_order, Float64),
+            load_matrix(f["generatorstorages/repairprobability"], region_order, Float64);
+            initial_soc = load_matrix(f["generatorstorages/initialsoc"], region_order, Float64)
+        )
+
+        region_genstor_idxs = makeidxlist(genstor_regions[region_order], n_regions)
+
+    else
+
+        generatorstorages = GeneratorStorages{N,L,T,P,E}()
+
+        region_genstor_idxs = fill(1:0, n_regions)
+
+    end
+
+
+    if has_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);
+            initial_borrowed_load = load_matrix(f["demandresponses/initialborrowedload"], region_order, Float64),
+            borrow_efficiency = load_matrix(f["demandresponses/borrowefficiency"], region_order, Float64),
+            payback_efficiency = 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
+
     return SystemModel(
         regions, interfaces,
         generators, region_gen_idxs,
@@ -332,9 +542,10 @@ function systemmodel_0_8(f::File)
         demandresponses, region_dr_idxs,
         lines, interface_line_idxs,
         timestamps, attrs)
-    
+
 end
 
+
 """
 Attempts to parse the file's "vX.Y.Z" version label into (x::Int, y::Int, z::Int).
 Errors if the label cannot be found or parsed as expected.
diff --git a/PRASFiles.jl/src/Systems/write.jl b/PRASFiles.jl/src/Systems/write.jl
index 1b3cb8c0..93b0e57c 100644
--- a/PRASFiles.jl/src/Systems/write.jl
+++ b/PRASFiles.jl/src/Systems/write.jl
@@ -164,6 +164,9 @@ function process_storages!(
 
     storages["repairprobability", deflate = compression] = sys.storages.μ
 
+    storages["initialsoc", deflate = compression] =
+         sys.storages.initial_soc
+
     return
 
 end
@@ -217,6 +220,8 @@ function process_generatorstorages!(
     generatorstorages["repairprobability", deflate = compression] =
         sys.generatorstorages.μ
 
+    generatorstorages["initialsoc", deflate = compression] =
+        sys.generatorstorages.initial_soc
     return
 
 end
@@ -261,6 +266,8 @@ function process_demandresponses!(
 
     demandresponses["repairprobability", deflate = compression] = sys.demandresponses.μ
 
+    demandresponses["initialborrowedload", deflate = compression] =
+        sys.demandresponses.initial_borrowed_load
     return
 
 end
diff --git a/PRASFiles.jl/test/runtests.jl b/PRASFiles.jl/test/runtests.jl
index 2249c6dd..626f6ccf 100644
--- a/PRASFiles.jl/test/runtests.jl
+++ b/PRASFiles.jl/test/runtests.jl
@@ -30,6 +30,25 @@ using JSON3
 
     end
 
+
+    @testset "Read version 0.7.0 to 0.8.1 compatibility" begin
+
+        path = dirname(@__FILE__)
+        version_0_7_0 = SystemModel(path * "/versioned_toy/toymodel_v0_7_0.pras")
+        version_0_8_1 = SystemModel(path * "/versioned_toy/toymodel_v0_8_1.pras")
+
+        @test version_0_7_0 == version_0_8_1
+        end
+
+    @testset "Read version 0.8.0 to 0.8.1 compatibility" begin
+
+        path = dirname(@__FILE__)
+        version_0_8_0 = SystemModel(path * "/versioned_toy/toymodel_v0_8_0.pras")
+        version_0_8_1 = SystemModel(path * "/versioned_toy/toymodel_v0_8_1.pras")
+
+        @test version_0_8_0 == version_0_8_1
+        end
+
     @testset "Run RTS-GMLC" begin
 
         assess(PRASFiles.rts_gmlc(), SequentialMonteCarlo(samples=100), Shortfall())
diff --git a/PRASFiles.jl/test/versioned_toy/toymodel_v0_7_0.pras b/PRASFiles.jl/test/versioned_toy/toymodel_v0_7_0.pras
new file mode 100644
index 00000000..080f1e69
Binary files /dev/null and b/PRASFiles.jl/test/versioned_toy/toymodel_v0_7_0.pras differ
diff --git a/PRASFiles.jl/test/versioned_toy/toymodel_v0_8_0.pras b/PRASFiles.jl/test/versioned_toy/toymodel_v0_8_0.pras
new file mode 100644
index 00000000..410bfd09
Binary files /dev/null and b/PRASFiles.jl/test/versioned_toy/toymodel_v0_8_0.pras differ
diff --git a/PRASFiles.jl/test/versioned_toy/toymodel_v0_8_1.pras b/PRASFiles.jl/test/versioned_toy/toymodel_v0_8_1.pras
new file mode 100644
index 00000000..b809daad
Binary files /dev/null and b/PRASFiles.jl/test/versioned_toy/toymodel_v0_8_1.pras differ
diff --git a/docs/src/SystemModel_HDF5_spec.md b/docs/src/SystemModel_HDF5_spec.md
index 898818c5..69e64d19 100644
--- a/docs/src/SystemModel_HDF5_spec.md
+++ b/docs/src/SystemModel_HDF5_spec.md
@@ -298,6 +298,11 @@ The `storages` group should also contain the following datasets describing
    transitions from forced outage to operational during a given simulation
    timestep, for each storage unit in each timeperiod. Unitless.
 
+An optional parameter is available to set the initial state of charge for the unit:
+
+ - `initialsoc`, as 64-bit floats representing the initial state of charge of
+   the unit as a fraction [0.0, 1.0] of energy capacity at timestep 1, for each storage unit.
+
 #### `generatorstorages` group
 
 Information relating to the combination generation-storage resources in the
@@ -362,6 +367,11 @@ generator-storage devices:
    transitions from forced outage to operational during a given simulation
    timestep, for each generator-storage unit in each timeperiod. Unitless.
 
+An optional parameter is available to set the initial state of charge for the unit:
+
+ - `initialsoc`, as 64-bit floats representing the initial state of charge of
+   the unit as a fraction [0.0, 1.0] of energy capacity at timestep 1, for each generator-storage unit.
+
 #### `demandresponses` group
 
 Information relating to the demand response only devices of the represented system is
@@ -406,6 +416,11 @@ The `demandresponse` group should also contain the following datasets describing
    transitions from forced outage to operational during a given simulation
    timestep, for each demand response unit in each timeperiod. Unitless.
 
+An optional parameter is available to set the initial state of borrowed load for the unit:
+
+ - `initialborrowedload`, as 64-bit floats representing the initial borrowed load
+   for the unit as a fraction [0.0, 1.0] of energy capacity at timestep 1, for each demand response unit.
+
 #### `interfaces` group
 
 Information relating to transmission interfaces between regions of the