diff --git a/src/PowerFlowData.jl b/src/PowerFlowData.jl
index 43470a12..f064a02c 100644
--- a/src/PowerFlowData.jl
+++ b/src/PowerFlowData.jl
@@ -1,3 +1,20 @@
+abstract type PowerFlowContainer end
+
+"""
+Trait signifying whether the `PowerFlowContainer` can represent multi-period data. Must be
+implemented for all concrete subtypes.
+"""
+supports_multi_period(x::PowerFlowContainer) =
+    throw(
+        IS.NotImplementedError(
+            "supports_multi_period must be implemented for $(typeof(x))"),
+    )
+
+"A `PowerFlowContainer` that represents its data as a `PSY.System`"
+abstract type SystemPowerFlowContainer <: PowerFlowContainer end
+
+get_system(container::SystemPowerFlowContainer) = container.system
+
 """
 Structure containing all the data required for the evaluation of the power
 flows and angles, as well as these ones.
@@ -51,7 +68,7 @@ flows and angles, as well as these ones.
 struct PowerFlowData{
     M <: PNM.PowerNetworkMatrix,
     N <: Union{PNM.PowerNetworkMatrix, Nothing},
-}
+} <: PowerFlowContainer
     bus_lookup::Dict{Int, Int}
     branch_lookup::Dict{String, Int}
     bus_activepower_injection::Matrix{Float64}
@@ -60,6 +77,7 @@ struct PowerFlowData{
     bus_reactivepower_withdrawals::Matrix{Float64}
     bus_reactivepower_bounds::Vector{Vector{Float64}}
     bus_type::Vector{PSY.ACBusTypes}
+    branch_type::Vector{DataType}
     bus_magnitude::Matrix{Float64}
     bus_angles::Matrix{Float64}
     branch_flow_values::Matrix{Float64}
@@ -70,6 +88,33 @@ struct PowerFlowData{
     neighbors::Vector{Set{Int}}
 end
 
+get_bus_lookup(pfd::PowerFlowData) = pfd.bus_lookup
+get_branch_lookup(pfd::PowerFlowData) = pfd.branch_lookup
+get_bus_activepower_injection(pfd::PowerFlowData) = pfd.bus_activepower_injection
+get_bus_reactivepower_injection(pfd::PowerFlowData) = pfd.bus_reactivepower_injection
+get_bus_activepower_withdrawals(pfd::PowerFlowData) = pfd.bus_activepower_withdrawals
+get_bus_reactivepower_withdrawals(pfd::PowerFlowData) = pfd.bus_reactivepower_withdrawals
+get_bus_reactivepower_bounds(pfd::PowerFlowData) = pfd.bus_reactivepower_bounds
+get_bus_type(pfd::PowerFlowData) = pfd.bus_type
+get_branch_type(pfd::PowerFlowData) = pfd.branch_type
+get_bus_magnitude(pfd::PowerFlowData) = pfd.bus_magnitude
+get_bus_angles(pfd::PowerFlowData) = pfd.bus_angles
+get_branch_flow_values(pfd::PowerFlowData) = pfd.branch_flow_values
+get_timestep_map(pfd::PowerFlowData) = pfd.timestep_map
+get_valid_ix(pfd::PowerFlowData) = pfd.valid_ix
+get_power_network_matrix(pfd::PowerFlowData) = pfd.power_network_matrix
+get_aux_network_matrix(pfd::PowerFlowData) = pfd.aux_network_matrix
+get_neighbor(pfd::PowerFlowData) = pfd.neighbors
+supports_multi_period(::PowerFlowData) = true
+
+function clear_injection_data!(pfd::PowerFlowData)
+    pfd.bus_activepower_injection .= 0.0
+    pfd.bus_reactivepower_injection .= 0.0
+    pfd.bus_activepower_withdrawals .= 0.0
+    pfd.bus_reactivepower_withdrawals .= 0.0
+    return
+end
+
 # AC Power Flow Data
 # TODO -> MULTI PERIOD: AC Power Flow Data
 function _calculate_neighbors(
@@ -99,12 +144,12 @@ NOTE: use it for AC power flow computations.
 - `sys::PSY.System`:
         container storing the system data to consider in the PowerFlowData
         structure.
-- `timesteps::Int`:
+- `time_steps::Int`:
         number of time periods to consider in the PowerFlowData structure. It
         defines the number of columns of the matrices used to store data.
         Default value = 1.
 - `timestep_names::Vector{String}`:
-        names of the time periods defines by the argmunet "timesteps". Default
+        names of the time periods defines by the argmunet "time_steps". Default
         value = String[].
 - `check_connectivity::Bool`:
         Perform connectivity check on the network matrix. Default value = true.
@@ -114,17 +159,17 @@ WARNING: functions for the evaluation of the multi-period AC PF still to be impl
 function PowerFlowData(
     ::ACPowerFlow,
     sys::PSY.System;
-    timesteps::Int = 1,
+    time_steps::Int = 1,
     timestep_names::Vector{String} = String[],
     check_connectivity::Bool = true)
 
     # assign timestep_names
     # timestep names are then allocated in a dictionary to map matrix columns
-    if timesteps != 0
+    if time_steps != 0
         if length(timestep_names) == 0
-            timestep_names = [string(i) for i in 1:timesteps]
-        elseif length(timestep_names) != timesteps
-            error("timestep_names field must have same length as timesteps")
+            timestep_names = [string(i) for i in 1:time_steps]
+        elseif length(timestep_names) != time_steps
+            error("timestep_names field must have same length as time_steps")
         end
     end
 
@@ -140,85 +185,41 @@ function PowerFlowData(
     n_branches = length(branches)
 
     bus_lookup = power_network_matrix.lookup[2]
-    branch_lookup =
-        Dict{String, Int}(PSY.get_name(b) => ix for (ix, b) in enumerate(branches))
-
-    # TODO: bus_type might need to also be a Matrix since the type can change for a particular scenario
-    bus_type = Vector{PSY.ACBusTypes}(undef, n_buses)
-    bus_angles = zeros(Float64, n_buses)
-    bus_magnitude = zeros(Float64, n_buses)
+    branch_lookup = Dict{String, Int}()
     temp_bus_map = Dict{Int, String}(
         PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.Bus, sys)
     )
-
-    _initialize_bus_data!(
-        bus_type,
-        bus_angles,
-        bus_magnitude,
-        temp_bus_map,
-        bus_lookup,
-        sys,
-    )
-
-    bus_activepower_injection = zeros(Float64, n_buses)
-    bus_reactivepower_injection = zeros(Float64, n_buses)
-    _get_injections!(
-        bus_activepower_injection,
-        bus_reactivepower_injection,
-        bus_lookup,
-        sys,
-    )
-
-    bus_activepower_withdrawals = zeros(Float64, n_buses)
-    bus_reactivepower_withdrawals = zeros(Float64, n_buses)
-    _get_withdrawals!(
-        bus_activepower_withdrawals,
-        bus_reactivepower_withdrawals,
-        bus_lookup,
-        sys,
-    )
-
-    # define fields as matrices whose number of columns is eqault to the number of timesteps
-    bus_activepower_injection_1 = zeros(n_buses, timesteps)
-    bus_reactivepower_injection_1 = zeros(n_buses, timesteps)
-    bus_activepower_withdrawals_1 = zeros(n_buses, timesteps)
-    bus_reactivepower_withdrawals_1 = zeros(n_buses, timesteps)
-    bus_magnitude_1 = zeros(n_buses, timesteps)
-    bus_angles_1 = zeros(n_buses, timesteps)
-    branch_flow_values_1 = zeros(n_branches, timesteps)
+    branch_types = Vector{DataType}(undef, n_branches)
+    for (ix, b) in enumerate(branches)
+        branch_lookup[PSY.get_name(b)] = ix
+        branch_types[ix] = typeof(b)
+    end
 
     bus_reactivepower_bounds = Vector{Vector{Float64}}(undef, n_buses)
     for i in 1:n_buses
         bus_reactivepower_bounds[i] = [0.0, 0.0]
     end
     _get_reactive_power_bound!(bus_reactivepower_bounds, bus_lookup, sys)
+    timestep_map = Dict(1 => "1")
+    valid_ix = setdiff(1:n_buses, ref_bus_positions)
+    neighbors = _calculate_neighbors(power_network_matrix)
+    aux_network_matrix = nothing
 
-    # initial values related to first timestep allocated in the first column
-    bus_activepower_injection_1[:, 1] .= bus_activepower_injection
-    bus_reactivepower_injection_1[:, 1] .= bus_reactivepower_injection
-    bus_activepower_withdrawals_1[:, 1] .= bus_activepower_withdrawals
-    bus_reactivepower_withdrawals_1[:, 1] .= bus_reactivepower_withdrawals
-    bus_magnitude_1[:, 1] .= bus_magnitude
-    bus_angles_1[:, 1] .= bus_angles
-    branch_flow_values_1[:, 1] .= zeros(n_branches)
-
-    return PowerFlowData(
+    return make_powerflowdata(
+        sys,
+        time_steps,
+        power_network_matrix,
+        aux_network_matrix,
+        n_buses,
+        n_branches,
         bus_lookup,
         branch_lookup,
-        bus_activepower_injection_1,
-        bus_reactivepower_injection_1,
-        bus_activepower_withdrawals_1,
-        bus_reactivepower_withdrawals_1,
+        temp_bus_map,
+        branch_types,
         bus_reactivepower_bounds,
-        bus_type,
-        bus_magnitude_1,
-        bus_angles_1,
-        branch_flow_values_1,
-        Dict(1 => "1"),
-        setdiff(1:n_buses, ref_bus_positions),
-        power_network_matrix,
-        nothing,
-        _calculate_neighbors(power_network_matrix),
+        timestep_map,
+        valid_ix,
+        neighbors,
     )
 end
 
@@ -236,12 +237,12 @@ NOTE: use it for DC power flow computations.
 - `sys::PSY.System`:
         container storing the system data to consider in the PowerFlowData
         structure.
-- `timesteps::Int`:
+- `time_steps::Int`:
         number of time periods to consider in the PowerFlowData structure. It
         defines the number of columns of the matrices used to store data.
         Default value = 1.
 - `timestep_names::Vector{String}`:
-        names of the time periods defines by the argmunet "timesteps". Default
+        names of the time periods defines by the argmunet "time_steps". Default
         value = String[].
 - `check_connectivity::Bool`:
         Perform connectivity check on the network matrix. Default value = true.
@@ -249,16 +250,16 @@ NOTE: use it for DC power flow computations.
 function PowerFlowData(
     ::DCPowerFlow,
     sys::PSY.System;
-    timesteps::Int = 1,
+    time_steps::Int = 1,
     timestep_names::Vector{String} = String[],
     check_connectivity::Bool = true)
 
     # assign timestep_names
     # timestep names are then allocated in a dictionary to map matrix columns
     if length(timestep_names) == 0
-        timestep_names = [string(i) for i in 1:timesteps]
-    elseif length(timestep_names) != timesteps
-        error("timestep_names field must have same length as timesteps")
+        timestep_names = [string(i) for i in 1:time_steps]
+    elseif length(timestep_names) != time_steps
+        error("timestep_names field must have same length as time_steps")
     end
 
     # get the network matrices
@@ -271,80 +272,22 @@ function PowerFlowData(
 
     bus_lookup = aux_network_matrix.lookup[1]
     branch_lookup = aux_network_matrix.lookup[2]
-    bus_type = Vector{PSY.ACBusTypes}(undef, n_buses)
-    bus_angles = zeros(Float64, n_buses)
-    bus_magnitude = zeros(Float64, n_buses)
     temp_bus_map = Dict{Int, String}(
         PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.ACBus, sys)
     )
-
-    _initialize_bus_data!(
-        bus_type,
-        bus_angles,
-        bus_magnitude,
-        temp_bus_map,
-        bus_lookup,
-        sys,
-    )
-
-    # define injection vectors related to the first timestep
-    bus_activepower_injection = zeros(Float64, n_buses)
-    bus_reactivepower_injection = zeros(Float64, n_buses)
-    _get_injections!(
-        bus_activepower_injection,
-        bus_reactivepower_injection,
-        bus_lookup,
-        sys,
-    )
-
-    bus_activepower_withdrawals = zeros(Float64, n_buses)
-    bus_reactivepower_withdrawals = zeros(Float64, n_buses)
-    _get_withdrawals!(
-        bus_activepower_withdrawals,
-        bus_reactivepower_withdrawals,
-        bus_lookup,
+    valid_ix = setdiff(1:n_buses, aux_network_matrix.ref_bus_positions)
+    return make_dc_powerflowdata(
         sys,
-    )
-
-    # initialize data
-    init_1 = zeros(n_buses, timesteps)
-    init_2 = zeros(n_branches, timesteps)
-
-    # define fields as matrices whose number of columns is eqault to the number of timesteps
-    bus_activepower_injection_1 = deepcopy(init_1)
-    bus_reactivepower_injection_1 = deepcopy(init_1)
-    bus_activepower_withdrawals_1 = deepcopy(init_1)
-    bus_reactivepower_withdrawals_1 = deepcopy(init_1)
-    bus_magnitude_1 = zeros(n_buses, 1)
-    bus_angles_1 = deepcopy(init_1)
-    branch_flow_values_1 = deepcopy(init_2)
-
-    # initial values related to first timestep allocated in the first column
-    bus_activepower_injection_1[:, 1] .= bus_activepower_injection
-    bus_reactivepower_injection_1[:, 1] .= bus_reactivepower_injection
-    bus_activepower_withdrawals_1[:, 1] .= bus_activepower_withdrawals
-    bus_reactivepower_withdrawals_1[:, 1] .= bus_reactivepower_withdrawals
-    bus_magnitude_1[:, 1] .= bus_magnitude  # for DC case same value accross all timesteps
-    bus_angles_1[:, 1] .= bus_angles
-    branch_flow_values_1[:, 1] .= zeros(n_branches)
-
-    return PowerFlowData(
-        bus_lookup,
-        branch_lookup,
-        bus_activepower_injection_1,
-        bus_reactivepower_injection_1,
-        bus_activepower_withdrawals_1,
-        bus_reactivepower_withdrawals_1,
-        Vector{Vector{Float64}}(),
-        bus_type,
-        bus_magnitude_1,
-        bus_angles_1,
-        branch_flow_values_1,
-        Dict(zip([i for i in 1:timesteps], timestep_names)),
-        setdiff(1:n_buses, aux_network_matrix.ref_bus_positions),
+        time_steps,
+        timestep_names,
         power_network_matrix,
         aux_network_matrix,
-        Vector{Set{Int}}(),
+        n_buses,
+        n_branches,
+        bus_lookup,
+        branch_lookup,
+        temp_bus_map,
+        valid_ix,
     )
 end
 
@@ -362,28 +305,29 @@ NOTE: use it for DC power flow computations.
 - `sys::PSY.System`:
         container storing the system data to consider in the PowerFlowData
         structure.
-- `timesteps::Int`:
+- `time_steps::Int`:
         number of time periods to consider in the PowerFlowData structure. It
         defines the number of columns of the matrices used to store data.
         Default value = 1.
 - `timestep_names::Vector{String}`:
-        names of the time periods defines by the argmunet "timesteps". Default
+        names of the time periods defines by the argmunet "time_steps". Default
         value = String[].
 """
+
 function PowerFlowData(
     ::PTDFDCPowerFlow,
     sys::PSY.System;
-    timesteps::Int = 1,
+    time_steps::Int = 1,
     timestep_names::Vector{String} = String[],
     check_connectivity::Bool = true)
 
     # assign timestep_names
     # timestep names are then allocated in a dictionary to map matrix columns
-    if timesteps != 0
+    if time_steps != 0
         if length(timestep_names) == 0
-            timestep_names = [string(i) for i in 1:timesteps]
-        elseif length(timestep_names) != timesteps
-            error("timestep_names field must have same length as timesteps")
+            timestep_names = [string(i) for i in 1:time_steps]
+        elseif length(timestep_names) != time_steps
+            error("timestep_names field must have same length as time_steps")
         end
     end
 
@@ -397,80 +341,22 @@ function PowerFlowData(
 
     bus_lookup = power_network_matrix.lookup[1]
     branch_lookup = power_network_matrix.lookup[2]
-    bus_type = Vector{PSY.ACBusTypes}(undef, n_buses)
-    bus_angles = zeros(Float64, n_buses)
-    bus_magnitude = zeros(Float64, n_buses)
     temp_bus_map = Dict{Int, String}(
         PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.Bus, sys)
     )
-
-    _initialize_bus_data!(
-        bus_type,
-        bus_angles,
-        bus_magnitude,
-        temp_bus_map,
-        bus_lookup,
-        sys,
-    )
-
-    # define injection vectors related to the first timestep
-    bus_activepower_injection = zeros(Float64, n_buses)
-    bus_reactivepower_injection = zeros(Float64, n_buses)
-    _get_injections!(
-        bus_activepower_injection,
-        bus_reactivepower_injection,
-        bus_lookup,
-        sys,
-    )
-
-    bus_activepower_withdrawals = zeros(Float64, n_buses)
-    bus_reactivepower_withdrawals = zeros(Float64, n_buses)
-    _get_withdrawals!(
-        bus_activepower_withdrawals,
-        bus_reactivepower_withdrawals,
-        bus_lookup,
+    valid_ix = setdiff(1:n_buses, aux_network_matrix.ref_bus_positions)
+    return make_dc_powerflowdata(
         sys,
-    )
-
-    # initialize data
-    init_1 = zeros(n_buses, timesteps)
-    init_2 = zeros(n_branches, timesteps)
-
-    # define fields as matrices whose number of columns is eqault to the number of timesteps
-    bus_activepower_injection_1 = deepcopy(init_1)
-    bus_reactivepower_injection_1 = deepcopy(init_1)
-    bus_activepower_withdrawals_1 = deepcopy(init_1)
-    bus_reactivepower_withdrawals_1 = deepcopy(init_1)
-    bus_magnitude_1 = zeros(n_buses, 1)
-    bus_angles_1 = deepcopy(init_1)
-    branch_flow_values_1 = deepcopy(init_2)
-
-    # initial values related to first timestep allocated in the first column
-    bus_activepower_injection_1[:, 1] .= bus_activepower_injection
-    bus_reactivepower_injection_1[:, 1] .= bus_reactivepower_injection
-    bus_activepower_withdrawals_1[:, 1] .= bus_activepower_withdrawals
-    bus_reactivepower_withdrawals_1[:, 1] .= bus_reactivepower_withdrawals
-    bus_magnitude_1[:, 1] .= bus_magnitude  # for DC case same value accross all timesteps
-    bus_angles_1[:, 1] .= bus_angles
-    branch_flow_values_1[:, 1] .= zeros(n_branches)
-
-    return PowerFlowData(
-        bus_lookup,
-        branch_lookup,
-        bus_activepower_injection_1,
-        bus_reactivepower_injection_1,
-        bus_activepower_withdrawals_1,
-        bus_reactivepower_withdrawals_1,
-        Vector{Vector{Float64}}(),
-        bus_type,
-        bus_magnitude_1,
-        bus_angles_1,
-        branch_flow_values_1,
-        Dict(zip([i for i in 1:timesteps], timestep_names)),
-        setdiff(1:n_buses, aux_network_matrix.ref_bus_positions),
+        time_steps,
+        timestep_names,
         power_network_matrix,
         aux_network_matrix,
-        Vector{Set{Int}}(),
+        n_buses,
+        n_branches,
+        bus_lookup,
+        branch_lookup,
+        temp_bus_map,
+        valid_ix,
     )
 end
 
@@ -488,28 +374,28 @@ NOTE: use it for DC power flow computations.
 - `sys::PSY.System`:
         container storing the system data to consider in the PowerFlowData
         structure.
-- `timesteps::Int`:
+- `time_steps::Int`:
         number of time periods to consider in the PowerFlowData structure. It
         defines the number of columns of the matrices used to store data.
         Default value = 1.
 - `timestep_names::Vector{String}`:
-        names of the time periods defines by the argmunet "timesteps". Default
+        names of the time periods defines by the argmunet "time_steps". Default
         value = String[].
 """
 function PowerFlowData(
     ::vPTDFDCPowerFlow,
     sys::PSY.System;
-    timesteps::Int = 1,
+    time_steps::Int = 1,
     timestep_names::Vector{String} = String[],
     check_connectivity::Bool = true)
 
     # assign timestep_names
     # timestep names are then allocated in a dictionary to map matrix columns
-    if timesteps != 0
+    if time_steps != 0
         if length(timestep_names) == 0
-            timestep_names = [string(i) for i in 1:timesteps]
-        elseif length(timestep_names) != timesteps
-            error("timestep_names field must have same length as timesteps")
+            timestep_names = [string(i) for i in 1:time_steps]
+        elseif length(timestep_names) != time_steps
+            error("timestep_names field must have same length as time_steps")
         end
     end
 
@@ -523,79 +409,45 @@ function PowerFlowData(
 
     bus_lookup = power_network_matrix.lookup[2]
     branch_lookup = power_network_matrix.lookup[1]
-    bus_type = Vector{PSY.ACBusTypes}(undef, n_buses)
-    bus_angles = zeros(Float64, n_buses)
-    bus_magnitude = zeros(Float64, n_buses)
     temp_bus_map = Dict{Int, String}(
         PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.Bus, sys)
     )
-
-    _initialize_bus_data!(
-        bus_type,
-        bus_angles,
-        bus_magnitude,
-        temp_bus_map,
-        bus_lookup,
+    valid_ix = setdiff(1:n_buses, aux_network_matrix.ref_bus_positions)
+    return make_dc_powerflowdata(
         sys,
-    )
-
-    # define injection vectors related to the first timestep
-    bus_activepower_injection = zeros(Float64, n_buses)
-    bus_reactivepower_injection = zeros(Float64, n_buses)
-    _get_injections!(
-        bus_activepower_injection,
-        bus_reactivepower_injection,
-        bus_lookup,
-        sys,
-    )
-
-    bus_activepower_withdrawals = zeros(Float64, n_buses)
-    bus_reactivepower_withdrawals = zeros(Float64, n_buses)
-    _get_withdrawals!(
-        bus_activepower_withdrawals,
-        bus_reactivepower_withdrawals,
-        bus_lookup,
-        sys,
-    )
-
-    # initialize data
-    init_1 = zeros(n_buses, timesteps)
-    init_2 = zeros(n_branches, timesteps)
-
-    # define fields as matrices whose number of columns is eqault to the number of timesteps
-    bus_activepower_injection_1 = deepcopy(init_1)
-    bus_reactivepower_injection_1 = deepcopy(init_1)
-    bus_activepower_withdrawals_1 = deepcopy(init_1)
-    bus_reactivepower_withdrawals_1 = deepcopy(init_1)
-    bus_magnitude_1 = zeros(n_buses, 1)
-    bus_angles_1 = deepcopy(init_1)
-    branch_flow_values_1 = deepcopy(init_2)
-
-    # initial values related to first timestep allocated in the first column
-    bus_activepower_injection_1[:, 1] .= bus_activepower_injection
-    bus_reactivepower_injection_1[:, 1] .= bus_reactivepower_injection
-    bus_activepower_withdrawals_1[:, 1] .= bus_activepower_withdrawals
-    bus_reactivepower_withdrawals_1[:, 1] .= bus_reactivepower_withdrawals
-    bus_magnitude_1[:, 1] .= bus_magnitude  # for DC case same value accross all timesteps
-    bus_angles_1[:, 1] .= bus_angles
-    branch_flow_values_1[:, 1] .= zeros(n_branches)
-
-    return PowerFlowData(
-        bus_lookup,
-        branch_lookup,
-        bus_activepower_injection_1,
-        bus_reactivepower_injection_1,
-        bus_activepower_withdrawals_1,
-        bus_reactivepower_withdrawals_1,
-        Vector{Vector{Float64}}(),
-        bus_type,
-        bus_magnitude_1,
-        bus_angles_1,
-        branch_flow_values_1,
-        Dict(zip([i for i in 1:timesteps], timestep_names)),
-        setdiff(1:n_buses, aux_network_matrix.ref_bus_positions),
+        time_steps,
+        timestep_names,
         power_network_matrix,
         aux_network_matrix,
-        Vector{Set{Int}}(),
+        n_buses,
+        n_branches,
+        bus_lookup,
+        branch_lookup,
+        temp_bus_map,
+        valid_ix,
     )
 end
+
+"""
+Create an appropriate `PowerFlowContainer` for the given `PowerFlowEvaluationModel` and initialize it from the given `PSY.System`.
+
+# Arguments:
+- `pfem::PowerFlowEvaluationModel`: power flow model to construct a container for (e.g., `DCPowerFlow()`)
+- `sys::PSY.System`: the system from which to initialize the power flow container
+- `time_steps::Int`: number of time periods to consider (default is `1`)
+- `timestep_names::Vector{String}`: names of the time periods defines by the argument "time_steps". Default value is `String[]`.
+- `check_connectivity::Bool`: Perform connectivity check on the network matrix. Default value is `true`.
+"""
+function make_power_flow_container end
+
+make_power_flow_container(pfem::ACPowerFlow, sys::PSY.System; kwargs...) =
+    PowerFlowData(pfem, sys; kwargs...)
+
+make_power_flow_container(pfem::DCPowerFlow, sys::PSY.System; kwargs...) =
+    PowerFlowData(pfem, sys; kwargs...)
+
+make_power_flow_container(pfem::PTDFDCPowerFlow, sys::PSY.System; kwargs...) =
+    PowerFlowData(pfem, sys; kwargs...)
+
+make_power_flow_container(pfem::vPTDFDCPowerFlow, sys::PSY.System; kwargs...) =
+    PowerFlowData(pfem, sys; kwargs...)
diff --git a/src/PowerFlows.jl b/src/PowerFlows.jl
index a7a8df92..1dead14e 100644
--- a/src/PowerFlows.jl
+++ b/src/PowerFlows.jl
@@ -7,6 +7,7 @@ export DCPowerFlow
 export ACPowerFlow
 export PTDFDCPowerFlow
 export vPTDFDCPowerFlow
+export PSSEExportPowerFlow
 export write_results
 export PSSEExporter
 export update_exporter!
@@ -33,12 +34,12 @@ include("common.jl")
 include("definitions.jl")
 include("powerflow_types.jl")
 include("PowerFlowData.jl")
+include("psse_export.jl")
 include("solve_dc_powerflow.jl")
 include("ac_power_flow.jl")
 include("ac_power_flow_jacobian.jl")
 include("nlsolve_ac_powerflow.jl")
 include("post_processing.jl")
-include("psse_export.jl")
 
 # Old PSSE Exporter
 import DataFrames: DataFrame
diff --git a/src/common.jl b/src/common.jl
index 6f33f91c..15bb6603 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -1,10 +1,6 @@
-function get_total_p(l::PSY.PowerLoad)
-    return PSY.get_active_power(l)
-end
-
-function get_total_q(l::PSY.PowerLoad)
-    return PSY.get_reactive_power(l)
-end
+_SingleComponentLoad = Union{PSY.PowerLoad, PSY.ExponentialLoad, PSY.InterruptiblePowerLoad}
+get_total_p(l::_SingleComponentLoad) = PSY.get_active_power(l)
+get_total_q(l::_SingleComponentLoad) = PSY.get_reactive_power(l)
 
 function get_total_p(l::PSY.StandardLoad)
     return PSY.get_constant_active_power(l) +
@@ -18,14 +14,6 @@ function get_total_q(l::PSY.StandardLoad)
            PSY.get_impedance_reactive_power(l)
 end
 
-function get_total_p(l::PSY.ExponentialLoad)
-    return PSY.get_active_power(l)
-end
-
-function get_total_q(l::PSY.ExponentialLoad)
-    return PSY.get_reactive_power(l)
-end
-
 function _get_injections!(
     bus_activepower_injection::Vector{Float64},
     bus_reactivepower_injection::Vector{Float64},
@@ -116,3 +104,133 @@ function my_mul_mt(
     end
     return y
 end
+
+function make_dc_powerflowdata(
+    sys,
+    time_steps,
+    timestep_names,
+    power_network_matrix,
+    aux_network_matrix,
+    n_buses,
+    n_branches,
+    bus_lookup,
+    branch_lookup,
+    temp_bus_map,
+    valid_ix,
+)
+    branch_types = Vector{DataType}(undef, length(branch_lookup))
+    for (ix, b) in enumerate(PNM.get_ac_branches(sys))
+        branch_types[ix] = typeof(b)
+    end
+    bus_reactivepower_bounds = Vector{Vector{Float64}}()
+    timestep_map = Dict(zip([i for i in 1:time_steps], timestep_names))
+    neighbors = Vector{Set{Int}}()
+    return make_powerflowdata(
+        sys,
+        time_steps,
+        power_network_matrix,
+        aux_network_matrix,
+        n_buses,
+        n_branches,
+        bus_lookup,
+        branch_lookup,
+        temp_bus_map,
+        branch_types,
+        bus_reactivepower_bounds,
+        timestep_map,
+        valid_ix,
+        neighbors,
+    )
+end
+
+function make_powerflowdata(
+    sys,
+    time_steps,
+    power_network_matrix,
+    aux_network_matrix,
+    n_buses,
+    n_branches,
+    bus_lookup,
+    branch_lookup,
+    temp_bus_map,
+    branch_types,
+    bus_reactivepower_bounds,
+    timestep_map,
+    valid_ix,
+    neighbors,
+)
+    # TODO: bus_type might need to also be a Matrix since the type can change for a particular scenario
+    bus_type = Vector{PSY.ACBusTypes}(undef, n_buses)
+    bus_angles = zeros(Float64, n_buses)
+    bus_magnitude = zeros(Float64, n_buses)
+
+    _initialize_bus_data!(
+        bus_type,
+        bus_angles,
+        bus_magnitude,
+        temp_bus_map,
+        bus_lookup,
+        sys,
+    )
+
+    # define injection vectors related to the first timestep
+    bus_activepower_injection = zeros(Float64, n_buses)
+    bus_reactivepower_injection = zeros(Float64, n_buses)
+    _get_injections!(
+        bus_activepower_injection,
+        bus_reactivepower_injection,
+        bus_lookup,
+        sys,
+    )
+
+    bus_activepower_withdrawals = zeros(Float64, n_buses)
+    bus_reactivepower_withdrawals = zeros(Float64, n_buses)
+    _get_withdrawals!(
+        bus_activepower_withdrawals,
+        bus_reactivepower_withdrawals,
+        bus_lookup,
+        sys,
+    )
+
+    # initialize data
+    init_1 = zeros(n_buses, time_steps)
+    init_2 = zeros(n_branches, time_steps)
+
+    # define fields as matrices whose number of columns is eqault to the number of time_steps
+    bus_activepower_injection_1 = deepcopy(init_1)
+    bus_reactivepower_injection_1 = deepcopy(init_1)
+    bus_activepower_withdrawals_1 = deepcopy(init_1)
+    bus_reactivepower_withdrawals_1 = deepcopy(init_1)
+    bus_magnitude_1 = deepcopy(init_1)
+    bus_angles_1 = deepcopy(init_1)
+    branch_flow_values_1 = deepcopy(init_2)
+
+    # initial values related to first timestep allocated in the first column
+    bus_activepower_injection_1[:, 1] .= bus_activepower_injection
+    bus_reactivepower_injection_1[:, 1] .= bus_reactivepower_injection
+    bus_activepower_withdrawals_1[:, 1] .= bus_activepower_withdrawals
+    bus_reactivepower_withdrawals_1[:, 1] .= bus_reactivepower_withdrawals
+    bus_magnitude_1[:, 1] .= bus_magnitude
+    bus_angles_1[:, 1] .= bus_angles
+    branch_flow_values_1[:, 1] .= zeros(n_branches)
+
+    return PowerFlowData(
+        bus_lookup,
+        branch_lookup,
+        bus_activepower_injection_1,
+        bus_reactivepower_injection_1,
+        bus_activepower_withdrawals_1,
+        bus_reactivepower_withdrawals_1,
+        bus_reactivepower_bounds,
+        bus_type,
+        branch_types,
+        bus_magnitude_1,
+        bus_angles_1,
+        branch_flow_values_1,
+        timestep_map,
+        valid_ix,
+        power_network_matrix,
+        aux_network_matrix,
+        neighbors,
+    )
+end
diff --git a/src/post_processing.jl b/src/post_processing.jl
index 8671726b..6c500a25 100644
--- a/src/post_processing.jl
+++ b/src/post_processing.jl
@@ -696,7 +696,8 @@ end
 """
 Modify the values in the given `System` to correspond to the given `PowerFlowData` such that
 if a new `PowerFlowData` is constructed from the resulting system it is the same as `data`.
-See also `write_powerflow_solution!`.
+See also `write_powerflow_solution!`. NOTE that this assumes that `data` was initialized
+from `sys` and then solved with no further modifications.
 """
 function update_system!(sys::PSY.System, data::PowerFlowData)
     for bus in PSY.get_components(PSY.Bus, sys)
@@ -719,5 +720,4 @@ function update_system!(sys::PSY.System, data::PowerFlowData)
             PSY.set_angle!(bus, data.bus_angles[data.bus_lookup[PSY.get_number(bus)]])
         end
     end
-    # TODO
 end
diff --git a/src/powerflow_types.jl b/src/powerflow_types.jl
index 03274840..e1b9828a 100644
--- a/src/powerflow_types.jl
+++ b/src/powerflow_types.jl
@@ -1,7 +1,15 @@
-Base.@kwdef struct ACPowerFlow
+abstract type PowerFlowEvaluationModel end
+
+Base.@kwdef struct ACPowerFlow <: PowerFlowEvaluationModel
     check_reactive_power_limits::Bool = false
 end
 
-struct DCPowerFlow end
-struct PTDFDCPowerFlow end
-struct vPTDFDCPowerFlow end
+struct DCPowerFlow <: PowerFlowEvaluationModel end
+struct PTDFDCPowerFlow <: PowerFlowEvaluationModel end
+struct vPTDFDCPowerFlow <: PowerFlowEvaluationModel end
+
+Base.@kwdef struct PSSEExportPowerFlow <: PowerFlowEvaluationModel
+    psse_version::Symbol
+    export_dir::AbstractString
+    write_comments::Bool = false
+end
diff --git a/src/psse_export.jl b/src/psse_export.jl
index 6125415b..bdf5c59b 100644
--- a/src/psse_export.jl
+++ b/src/psse_export.jl
@@ -39,11 +39,13 @@ const PSSE_GROUPS_33 = [
     "Q Record",
 ]
 
+const PSSE_DEFAULT_EXPORT_NAME = "export"
+
 # TODO move this to IS
 """
 A context manager similar to `Logging.with_logger` that sets the system's units to the given
 value, executes the function, then sets them back. Suppresses logging below `Warn` from
-internal calls to `set_units_base_system!`.
+internal calls to `set_units_base_system!`. Not thread safe.
 """
 function with_units(f::Function, sys::System, units::Union{PSY.UnitSystem, String})
     old_units = PSY.get_units_base(sys)
@@ -59,6 +61,49 @@ function with_units(f::Function, sys::System, units::Union{PSY.UnitSystem, Strin
     end
 end
 
+"""
+Make a `deepcopy` of the `System` except replace the time series manager and supplemental
+attribute manager with blank versions so these are not copied.
+"""
+function deepcopy_system_no_time_series_no_supplementals(sys::System)
+    old_time_series_manager = sys.data.time_series_manager
+    old_supplemental_attribute_manager = sys.data.supplemental_attribute_manager
+
+    new_time_series_manager = IS.TimeSeriesManager(
+        IS.InMemoryTimeSeriesStorage(),
+        IS.TimeSeriesMetadataStore(),
+        true,
+    )
+    new_supplemental_attribute_manager = IS.SupplementalAttributeManager()
+
+    sys.data.time_series_manager = new_time_series_manager
+    sys.data.supplemental_attribute_manager = new_supplemental_attribute_manager
+
+    old_refs = Dict{Tuple{DataType, String}, IS.SharedSystemReferences}()
+    for comp in PSY.iterate_components(sys)
+        old_refs[(typeof(comp), PSY.get_name(comp))] =
+            comp.internal.shared_system_references
+        new_refs = IS.SharedSystemReferences(;
+            time_series_manager = new_time_series_manager,
+            supplemental_attribute_manager = new_supplemental_attribute_manager,
+        )
+        IS.set_shared_system_references!(comp, new_refs)
+    end
+
+    new_sys = try
+        deepcopy(sys)
+    finally
+        sys.data.time_series_manager = old_time_series_manager
+        sys.data.supplemental_attribute_manager = old_supplemental_attribute_manager
+
+        for comp in PSY.iterate_components(sys)
+            IS.set_shared_system_references!(comp,
+                old_refs[(typeof(comp), PSY.get_name(comp))])
+        end
+    end
+    return new_sys
+end
+
 """
 Structure to perform an export from a Sienna System, plus optional updates from
 `PowerFlowData`, to the PSS/E format. Construct from a `System` and a PSS/E version, update
@@ -70,20 +115,32 @@ using `update_exporter` with any new data as relevant, and perform the export wi
     flow-related values but may not fundamentally alter the system
   - `psse_version::Symbol`: the version of PSS/E to target, must be one of
     `PSSE_EXPORT_SUPPORTED_VERSIONS`
-  - `write_comments::Bool`: whether to add the customary-but-not-in-spec-annotations after a
-    slash on the first line and at group boundaries
+  - `write_comments::Bool` = false: whether to add the customary-but-not-in-spec-annotations
+    after a slash on the first line and at group boundaries
+  - `name::AbstractString = "export"`: the base name of the export
+  - `step::Union{Nothing, Integer, Tuple{Vararg{Integer}}} = nothing`: optional step number
+    or tuple of step numbers (e.g., step and timestamp within step) to append to the base
+    export name. User is responsible for updating the step.
+  - `overwrite::Bool = false`: `true` to silently overwrite existing exports, `false` to
+    throw an error if existing results are encountered
 """
-mutable struct PSSEExporter
-    # Internal fields are very much subject to change as I iterate on the best way to do
-    # this! For instance, the final version will almost certainly not store an entire System
+mutable struct PSSEExporter <: SystemPowerFlowContainer
     system::PSY.System
     psse_version::Symbol
+    export_dir::AbstractString
     write_comments::Bool
+    name::AbstractString
+    step::Union{Nothing, Integer, Tuple{Vararg{Integer}}}
+    overwrite::Bool
 
     function PSSEExporter(
         base_system::PSY.System,
-        psse_version::Symbol;
+        psse_version::Symbol,
+        export_dir::AbstractString;
         write_comments::Bool = false,
+        name::AbstractString = PSSE_DEFAULT_EXPORT_NAME,
+        step::Union{Nothing, Integer, Tuple{Vararg{Integer}}} = nothing,
+        overwrite::Bool = false,
     )
         (psse_version in PSSE_EXPORT_SUPPORTED_VERSIONS) ||
             throw(
@@ -91,11 +148,14 @@ mutable struct PSSEExporter
                     "PSS/E version $psse_version is not supported, must be one of $PSSE_EXPORT_SUPPORTED_VERSIONS",
                 ),
             )
-        system = deepcopy(base_system)
-        new(system, psse_version, write_comments)
+        system = deepcopy_system_no_time_series_no_supplementals(base_system)
+        mkpath(export_dir)
+        new(system, psse_version, export_dir, write_comments, name, step, overwrite)
     end
 end
 
+supports_multi_period(::PSSEExporter) = false
+
 function _validate_same_system(sys1::PSY.System, sys2::PSY.System)
     return IS.get_uuid(PSY.get_internal(sys1)) == IS.get_uuid(PSY.get_internal(sys2))
 end
@@ -128,7 +188,7 @@ function update_exporter!(exporter::PSSEExporter, data::PSY.System)
             "System passed to update_exporter must be the same system as the one with which the exporter was constructed, just with different values",
         ),
     )
-    exporter.system = deepcopy(data)
+    exporter.system = deepcopy_system_no_time_series_no_supplementals(data)
     return
 end
 
@@ -324,8 +384,16 @@ function _write_raw(
         NAME = _psse_quote_string(bus_name_mapping[PSY.get_name(bus)])
         BASKV = PSY.get_base_voltage(bus)
         IDE = PSSE_BUS_TYPE_MAP[PSY.get_bustype(bus)]
-        AREA = md["area_mapping"][PSY.get_name(PSY.get_area(bus))]
-        ZONE = md["zone_number_mapping"][PSY.get_name(PSY.get_load_zone(bus))]
+        AREA = if isnothing(PSY.get_area(bus))
+            PSSE_DEFAULT
+        else
+            md["area_mapping"][PSY.get_name(PSY.get_area(bus))]
+        end
+        ZONE = if isnothing(PSY.get_load_zone(bus))
+            PSSE_DEFAULT
+        else
+            md["zone_number_mapping"][PSY.get_name(PSY.get_load_zone(bus))]
+        end
         OWNER = PSSE_DEFAULT
         VM = PSY.get_magnitude(bus)
         VA = rad2deg(PSY.get_angle(bus))
@@ -842,24 +910,33 @@ function _write_raw(::Val{T}, io::IO, md::AbstractDict, exporter::PSSEExporter)
     _write_skip_group(io, md, exporter, group_name)
 end
 
+_step_to_string(::Nothing) = ""
+_step_to_string(step::Integer) = "_$step"
+_step_to_string(step::Tuple{Vararg{Integer}}) = "_" * join(step, "_")
+
 "Peform an export from the data contained in a `PSSEExporter` to the PSS/E file format."
 function write_export(
     exporter::PSSEExporter,
-    scenario_name::AbstractString,
-    year::Int,
-    export_location::AbstractString,
+    name::AbstractString;
+    overwrite = false,
 )
+    name *= _step_to_string(exporter.step)
     # Construct paths
-    export_dir = joinpath(export_location, "Raw_Export", scenario_name, string(year))
-    mkpath(export_dir)
-    @info "Exporting to $export_dir"
-    raw_path = joinpath(export_dir, "$scenario_name.raw")
-    md_path = joinpath(export_dir, "$(scenario_name)_metadata.json")
+    export_subdir = joinpath(exporter.export_dir, name)
+    dir_exists = isdir(export_subdir)
+    (dir_exists && !overwrite) && throw(
+        ArgumentError(
+            "Target export directory $(abspath(export_subdir)) already exists; specify `overwrite = true` if it should be overwritten",
+        ),
+    )
+    dir_exists || mkdir(export_subdir)
+    @info "Exporting $name to $export_subdir"
+    raw_path, md_path = get_psse_export_paths(export_subdir)
 
     # Build export files in buffers
     raw = IOBuffer()
     md = OrderedDict()
-    md["case_name"] = "$(scenario_name)_$(year)"
+    md["case_name"] = name
     md["export_settings"] = OrderedDict("sources_as_generators" => true)
     # These mappings are accessed in e.g. _write_bus_data via the metadata
     md["area_mapping"] = _psse_container_numbers(
@@ -886,15 +963,20 @@ function write_export(
     open(file -> JSON3.pretty(file, md), md_path; truncate = true)
 end
 
+write_export(exporter::PSSEExporter; kwargs...) =
+    write_export(
+        exporter,
+        exporter.name;
+        merge(Dict(:overwrite => exporter.overwrite), kwargs)...,
+    )
+
 "Calculate the paths of the (raw, metadata) files that would be written by a certain call to `write_export`"
 function get_psse_export_paths(
-    scenario_name::AbstractString,
-    year::Int,
-    export_location::AbstractString,
+    export_subdir::AbstractString,
 )
-    base_path = joinpath(export_location, "Raw_Export", string(scenario_name), string(year))
-    raw_path = joinpath(base_path, "$scenario_name.raw")
-    metadata_path = joinpath(base_path, "$(scenario_name)_metadata.json")
+    name = last(splitdir(export_subdir))
+    raw_path = joinpath(export_subdir, "$name.raw")
+    metadata_path = joinpath(export_subdir, "$(name)_metadata.json")
     return (raw_path, metadata_path)
 end
 
@@ -956,7 +1038,7 @@ function fix_load_zone_names!(sys::PSY.System, md::Dict)
     # `collect` is necessary due to https://github.com/NREL-Sienna/PowerSystems.jl/issues/1161
     for load_zone in collect(PSY.get_components(PSY.LoadZone, sys))
         old_name = PSY.get_name(load_zone)
-        new_name = lz_map[parse(Int64, old_name)]
+        new_name = get(lz_map, parse(Int64, old_name), old_name)
         (old_name != new_name) && PSY.set_name!(sys, load_zone, new_name)
     end
 end
@@ -1020,3 +1102,15 @@ function PSY.System(raw_path::AbstractString, md::Dict)
     # TODO remap everything else! Should be reading all the keys in `md`
     return sys
 end
+
+# TODO handle kwargs
+make_power_flow_container(pfem::PSSEExportPowerFlow, sys::PSY.System; kwargs...) =
+    PSSEExporter(
+        sys,
+        pfem.psse_version,
+        pfem.export_dir;
+        write_comments = pfem.write_comments,
+        step = (0, 0),
+    )
+
+solve_powerflow!(exporter::PSSEExporter) = write_export(exporter)
diff --git a/src/solve_dc_powerflow.jl b/src/solve_dc_powerflow.jl
index 19e8f9fb..a1bb0f18 100644
--- a/src/solve_dc_powerflow.jl
+++ b/src/solve_dc_powerflow.jl
@@ -11,8 +11,17 @@ const PTDFPowerFlowData = PowerFlowData{
     },
 }
 
-# ? change this to have a more detailed definition ?
-const vPTDFPowerFlowData = PowerFlowData{}
+const vPTDFPowerFlowData = PowerFlowData{
+    PNM.VirtualPTDF{
+        Tuple{Vector{String}, Vector{Int64}},
+        Tuple{Dict{String, Int64}, Dict{Int64, Int64}},
+    },
+    PNM.ABA_Matrix{
+        Tuple{Vector{Int64}, Vector{Int64}},
+        Tuple{Dict{Int64, Int64}, Dict{Int64, Int64}},
+        PNM.KLU.KLUFactorization{Float64, Int64},
+    },
+}
 
 const ABAPowerFlowData = PowerFlowData{
     PNM.ABA_Matrix{
@@ -30,7 +39,7 @@ Evaluates the power flows on each system's branch and updates the PowerFlowData
 
 # Arguments:
 - `data::PTDFPowerFlowData`:
-        PTDFPowerFlowData structure containig all the information related to the system's power flow.
+        PTDFPowerFlowData structure containing all the information related to the system's power flow.
 """
 function solve_powerflow!(
     data::PTDFPowerFlowData,
@@ -50,7 +59,7 @@ Evaluates the power flows on each system's branch and updates the PowerFlowData
 
 # Arguments:
 - `data::vPTDFPowerFlowData`:
-        vPTDFPowerFlowData structure containig all the information related to the system's power flow.
+        vPTDFPowerFlowData structure containing all the information related to the system's power flow.
 """
 function solve_powerflow!(
     data::vPTDFPowerFlowData,
@@ -75,7 +84,7 @@ Evaluates the power flows on each system's branch and updates the PowerFlowData
 
 # Arguments:
 - `data::ABAPowerFlowData`:
-        ABAPowerFlowData structure containig all the information related to the system's power flow.
+        ABAPowerFlowData structure containing all the information related to the system's power flow.
 """
 # DC flow: ABA and BA case
 function solve_powerflow!(
diff --git a/test/runtests.jl b/test/runtests.jl
index 41e6e851..adc16f40 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -125,7 +125,7 @@ function run_tests()
         end
 
         # Testing Topological components of the schema
-        @time @testset "Begin PowerSystems tests" begin
+        @time @testset "Begin PowerFlows tests" begin
             @includetests ARGS
         end
 
diff --git a/test/test_multiperiod_dc_powerflow.jl b/test/test_multiperiod_dc_powerflow.jl
index 7f252e83..f81ef446 100644
--- a/test/test_multiperiod_dc_powerflow.jl
+++ b/test/test_multiperiod_dc_powerflow.jl
@@ -25,10 +25,10 @@
     ##############################################################################
 
     # create structure for multi-period case
-    timesteps = 24
-    data_1 = PowerFlowData(DCPowerFlow(), sys; timesteps = timesteps)
-    data_2 = PowerFlowData(PTDFDCPowerFlow(), sys; timesteps = timesteps)
-    data_3 = PowerFlowData(vPTDFDCPowerFlow(), sys; timesteps = timesteps)
+    time_steps = 24
+    data_1 = PowerFlowData(DCPowerFlow(), sys; time_steps = time_steps)
+    data_2 = PowerFlowData(PTDFDCPowerFlow(), sys; time_steps = time_steps)
+    data_3 = PowerFlowData(vPTDFDCPowerFlow(), sys; time_steps = time_steps)
 
     # allocate data from csv
     injs = Matrix(injections)
diff --git a/test/test_powerflow_data.jl b/test/test_powerflow_data.jl
index 38c14039..75964cb7 100644
--- a/test/test_powerflow_data.jl
+++ b/test/test_powerflow_data.jl
@@ -1,18 +1,20 @@
 @testset "PowerFlowData" begin
     sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
-    PowerFlowData(ACPowerFlow(), sys)
-    PowerFlowData(DCPowerFlow(), sys)
-    PowerFlowData(PTDFDCPowerFlow(), sys)
-    PowerFlowData(vPTDFDCPowerFlow(), sys)
+    @test PowerFlowData(ACPowerFlow(), sys) isa PF.ACPowerFlowData
+    @test PowerFlowData(DCPowerFlow(), sys) isa PF.ABAPowerFlowData
+    @test PowerFlowData(PTDFDCPowerFlow(), sys) isa PF.PTDFPowerFlowData
+    @test PowerFlowData(vPTDFDCPowerFlow(), sys) isa PF.vPTDFPowerFlowData
 end
 
 @testset "PowerFlowData multiperiod" begin
     sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
-    timesteps = 24
+    time_steps = 24
     # TODO: "multiperiod AC still to implement"
-    PowerFlowData(DCPowerFlow(), sys; timesteps = timesteps)
-    PowerFlowData(PTDFDCPowerFlow(), sys; timesteps = timesteps)
-    PowerFlowData(vPTDFDCPowerFlow(), sys; timesteps = timesteps)
+    @test PowerFlowData(DCPowerFlow(), sys; time_steps = time_steps) isa PF.ABAPowerFlowData
+    @test PowerFlowData(PTDFDCPowerFlow(), sys; time_steps = time_steps) isa
+          PF.PTDFPowerFlowData
+    @test PowerFlowData(vPTDFDCPowerFlow(), sys; time_steps = time_steps) isa
+          PF.vPTDFPowerFlowData
 end
 
 @testset "System <-> PowerFlowData round trip" begin
@@ -49,8 +51,11 @@ end
         data_modified,
     )
 
-    # The big one: update_system! with modified PowerFlowData should result in sys_modified
+    # The big one: update_system! with modified PowerFlowData should result in sys_modified,
+    # modulo information that is inherently lost in the PowerFlowData representation
     sys_modify_updated = deepcopy(sys_original)
     PF.update_system!(sys_modify_updated, data_modified)
-    @test IS.compare_values(powerflow_match_fn, sys_modify_updated, sys_modified)
+    sys_mod_redist = deepcopy(sys_modified)
+    PF.update_system!(sys_mod_redist, PowerFlowData(ACPowerFlow(), sys_mod_redist))
+    @test IS.compare_values(powerflow_match_fn, sys_modify_updated, sys_mod_redist)
 end
diff --git a/test/test_psse_export.jl b/test/test_psse_export.jl
index cbbb345c..bc87b3a6 100644
--- a/test/test_psse_export.jl
+++ b/test/test_psse_export.jl
@@ -287,23 +287,23 @@ function read_system_and_metadata(raw_path, metadata_path)
     return sys, md
 end
 
-read_system_and_metadata(scenario_name, year, export_location) = read_system_and_metadata(
-    get_psse_export_paths(scenario_name, year, export_location)...)
+read_system_and_metadata(export_subdir) = read_system_and_metadata(
+    get_psse_export_paths(export_subdir)...)
 
 function test_psse_round_trip(
     sys::System,
     exporter::PSSEExporter,
     scenario_name::AbstractString,
-    year::Int,
     export_location::AbstractString;
     do_power_flow_test = true,
     exclude_reactive_flow = false,
 )
-    raw_path, metadata_path = get_psse_export_paths(scenario_name, year, export_location)
+    raw_path, metadata_path =
+        get_psse_export_paths(joinpath(export_location, scenario_name))
     @test !isfile(raw_path)
     @test !isfile(metadata_path)
 
-    write_export(exporter, scenario_name, year, export_location)
+    write_export(exporter, scenario_name)
     @test isfile(raw_path)
     @test isfile(metadata_path)
 
@@ -356,6 +356,7 @@ end
 # I test so much, my tests have tests
 @testset "Test system comparison utilities" begin
     sys = load_test_system()
+    isnothing(sys) && return
 
     @test compare_systems_wrapper(sys, sys)
     @test compare_systems_wrapper(sys, deepcopy(sys))
@@ -363,21 +364,22 @@ end
 
 @testset "PSSE Exporter with system_240[32].json, v33" begin
     sys = load_test_system()
+    isnothing(sys) && return
 
     # PSS/E version must be one of the supported ones
-    @test_throws ArgumentError PSSEExporter(sys, :vNonexistent)
+    @test_throws ArgumentError PSSEExporter(sys, :vNonexistent, test_psse_export_dir)
 
     # Reimported export should be comparable to original system
-    exporter = PSSEExporter(sys, :v33)
     export_location = joinpath(test_psse_export_dir, "v33", "system_240")
-    test_psse_round_trip(sys, exporter, "basic", 2024, export_location;
+    exporter = PSSEExporter(sys, :v33, export_location)
+    test_psse_round_trip(sys, exporter, "basic", export_location;
         exclude_reactive_flow = true)  # TODO why is reactive flow not matching?
 
     # Exporting the exact same thing again should result in the exact same files
-    write_export(exporter, "basic2", 2024, export_location)
+    write_export(exporter, "basic2")
     test_psse_export_strict_equality(
-        get_psse_export_paths("basic", 2024, export_location)...,
-        get_psse_export_paths("basic2", 2024, export_location)...)
+        get_psse_export_paths(joinpath(export_location, "basic"))...,
+        get_psse_export_paths(joinpath(export_location, "basic2"))...)
 
     # Updating with a completely different system should fail
     different_system = build_system(PSITestSystems, "c_sys5_all_components")
@@ -385,18 +387,19 @@ end
 
     # Updating with the exact same system should result in the exact same files
     update_exporter!(exporter, sys)
-    write_export(exporter, "basic3", 2024, export_location)
+    write_export(exporter, "basic3")
     test_psse_export_strict_equality(
-        get_psse_export_paths("basic", 2024, export_location)...,
-        get_psse_export_paths("basic3", 2024, export_location)...)
+        get_psse_export_paths(joinpath(export_location, "basic"))...,
+        get_psse_export_paths(joinpath(export_location, "basic3"))...)
 
     # Updating with changed value should result in a different reimport (System version)
     sys2 = deepcopy(sys)
     line_to_change = first(get_components(Line, sys2))
     set_rating!(line_to_change, get_rating(line_to_change) * 12345.6)
     update_exporter!(exporter, sys2)
-    write_export(exporter, "basic4", 2024, export_location)
-    reread_sys2, sys2_metadata = read_system_and_metadata("basic4", 2024, export_location)
+    write_export(exporter, "basic4")
+    reread_sys2, sys2_metadata =
+        read_system_and_metadata(joinpath(export_location, "basic4"))
     @test compare_systems_wrapper(sys2, reread_sys2, sys2_metadata)
     @test_logs((:error, r"Mismatch on rate"), (:error, r"values do not match"),
         match_mode = :any, min_level = Logging.Error,
@@ -409,19 +412,19 @@ end
     set_units_base_system!(sys, UnitSystem.SYSTEM_BASE)
 
     # PSS/E version must be one of the supported ones
-    @test_throws ArgumentError PSSEExporter(sys, :vNonexistent)
+    @test_throws ArgumentError PSSEExporter(sys, :vNonexistent, test_psse_export_dir)
 
     # Reimported export should be comparable to original system
-    exporter = PSSEExporter(sys, :v33)
     export_location = joinpath(test_psse_export_dir, "v33", "rts_gmlc")
-    test_psse_round_trip(sys, exporter, "basic", 2024, export_location;
+    exporter = PSSEExporter(sys, :v33, export_location)
+    test_psse_round_trip(sys, exporter, "basic", export_location;
         exclude_reactive_flow = true)  # TODO why is reactive flow not matching?
 
     # Exporting the exact same thing again should result in the exact same files
-    write_export(exporter, "basic2", 2024, export_location)
+    write_export(exporter, "basic2")
     test_psse_export_strict_equality(
-        get_psse_export_paths("basic", 2024, export_location)...,
-        get_psse_export_paths("basic2", 2024, export_location)...)
+        get_psse_export_paths(joinpath(export_location, "basic"))...,
+        get_psse_export_paths(joinpath(export_location, "basic2"))...)
 
     # Updating with a completely different system should fail
     different_system = build_system(PSITestSystems, "c_sys5_all_components")
@@ -429,17 +432,18 @@ end
 
     # Updating with the exact same system should result in the exact same files
     update_exporter!(exporter, sys)
-    write_export(exporter, "basic3", 2024, export_location)
+    write_export(exporter, "basic3")
     test_psse_export_strict_equality(
-        get_psse_export_paths("basic", 2024, export_location)...,
-        get_psse_export_paths("basic3", 2024, export_location)...)
+        get_psse_export_paths(joinpath(export_location, "basic"))...,
+        get_psse_export_paths(joinpath(export_location, "basic3"))...)
 
     # Updating with changed value should result in a different reimport (System version)
     sys2 = deepcopy(sys)
     modify_rts_system!(sys2)
     update_exporter!(exporter, sys2)
-    write_export(exporter, "basic4", 2024, export_location)
-    reread_sys2, sys2_metadata = read_system_and_metadata("basic4", 2024, export_location)
+    write_export(exporter, "basic4")
+    reread_sys2, sys2_metadata =
+        read_system_and_metadata(joinpath(export_location, "basic4"))
     @test compare_systems_wrapper(sys2, reread_sys2, sys2_metadata)
     @test_logs((:error, r"values do not match"),
         (:error, r"Mismatch on active_power"), (:error, r"Mismatch on reactive_power"),
@@ -449,14 +453,15 @@ end
     test_power_flow(sys2, reread_sys2; exclude_reactive_flow = true)  # TODO why is reactive flow not matching?
 
     # Updating with changed value should result in a different reimport (PowerFlowData version)
-    exporter = PSSEExporter(sys, :v33)
+    exporter = PSSEExporter(sys, :v33, export_location)
     pf2 = PowerFlowData(ACPowerFlow(), sys)
     # This modifies the PowerFlowData in the same way that modify_rts_system! modifies the
     # system, so the reimport should be comparable to sys2 from above
     modify_rts_powerflow!(pf2)
     update_exporter!(exporter, pf2)
-    write_export(exporter, "basic5", 2024, export_location)
-    reread_sys3, sys3_metadata = read_system_and_metadata("basic5", 2024, export_location)
+    write_export(exporter, "basic5")
+    reread_sys3, sys3_metadata =
+        read_system_and_metadata(joinpath(export_location, "basic5"))
     @test compare_systems_wrapper(sys2, reread_sys3, sys3_metadata;
         exclude_reactive_power = true)  # TODO why is reactive power not matching?
     @test_logs((:error, r"values do not match"),
@@ -467,8 +472,8 @@ end
     test_power_flow(sys2, reread_sys3; exclude_reactive_flow = true)  # TODO why is reactive flow not matching?
 
     # Exporting with write_comments should be comparable to original system
-    exporter = PSSEExporter(sys, :v33; write_comments = true)
-    test_psse_round_trip(sys, exporter, "basic6", 2024, export_location;
+    exporter = PSSEExporter(sys, :v33, export_location; write_comments = true)
+    test_psse_round_trip(sys, exporter, "basic6", export_location;
         exclude_reactive_flow = true)  # TODO why is reactive flow not matching?
 end