diff --git a/src/post_processing.jl b/src/post_processing.jl
index 9e759f19..ee5b73e6 100644
--- a/src/post_processing.jl
+++ b/src/post_processing.jl
@@ -678,3 +678,30 @@ function write_results(
 
     return Dict("bus_results" => bus_df, "flow_results" => branch_df)
 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!`.
+"""
+function update_system!(sys::PSY.System, data::PowerFlowData)
+    for bus in PSY.get_components(PSY.Bus, sys)
+        if bus.bustype == PSY.ACBusTypes.REF
+            P_gen = data.bus_activepower_injection[data.bus_lookup[PSY.get_number(bus)]]
+            Q_gen = data.bus_reactivepower_injection[data.bus_lookup[PSY.get_number(bus)]]
+            _power_redistribution_ref(sys, P_gen, Q_gen, bus)
+            # elseif bus.bustype == PSY.ACBusTypes.PV
+            #     # TODO. This is the write_powerflow_solution! logic:
+            #     Q_gen = result[2 * ix - 1]
+            #     bus.angle = result[2 * ix]
+            #     _reactive_power_redistribution_pv(sys, Q_gen, bus)
+            # elseif bus.bustype == PSY.ACBusTypes.PQ
+            #     # TODO. This is the write_powerflow_solution! logic:
+            #     Vm = result[2 * ix - 1]
+            #     θ = result[2 * ix]
+            #     PSY.set_magnitude!(bus, Vm)
+            #     PSY.set_angle!(bus, θ)
+        end
+    end
+    # TODO
+end
diff --git a/test/test_powerflow_data.jl b/test/test_powerflow_data.jl
index 1f352156..03c124ed 100644
--- a/test/test_powerflow_data.jl
+++ b/test/test_powerflow_data.jl
@@ -14,3 +14,40 @@ end
     PowerFlowData(PTDFDCPowerFlow(), sys; timesteps = timesteps)
     PowerFlowData(vPTDFDCPowerFlow(), sys; timesteps = timesteps)
 end
+
+@testset "System <-> PowerFlowData round trip" begin
+    # TODO currently only tested with ACPowerFlow
+    # TODO test that update_system! errors if the PowerFlowData doesn't correspond to the system
+
+    sys_original = build_system(PSISystems, "RTS_GMLC_DA_sys")
+    data_original = PowerFlowData(ACPowerFlow(), sys_original)
+
+    sys_modified = deepcopy(sys_original)
+    modify_rts_system!(sys_modified)
+    data_modified = PowerFlowData(ACPowerFlow(), sys_original)
+    modify_rts_powerflow!(data_modified)
+
+    # update_system! with unmodified PowerFlowData should result in system that yields unmodified PowerFlowData
+    # (NOTE does NOT necessarily yield original system due to power redistribution)
+    sys_null_updated = deepcopy(sys_original)
+    PF.update_system!(sys_null_updated, data_original)
+    data_null_updated = PowerFlowData(ACPowerFlow(), sys_null_updated)
+    @test IS.compare_values(data_null_updated, data_original; match_fn = powerflow_match_fn)
+
+    # Modified versions should not be the same as unmodified versions
+    @test !@test_logs((:error, r"values do not match"),
+        match_mode = :any, min_level = Logging.Error,
+        IS.compare_values(data_original, data_modified; match_fn = powerflow_match_fn))
+    @test !@test_logs((:error, r"values do not match"),
+        match_mode = :any, min_level = Logging.Error,
+        IS.compare_values(sys_original, sys_modified; match_fn = powerflow_match_fn))
+
+    # Constructing PowerFlowData from modified system should result in data_modified
+    @test IS.compare_values(PowerFlowData(ACPowerFlow(), sys_modified), data_modified;
+        match_fn = powerflow_match_fn)
+
+    # The big one: update_system! with modified PowerFlowData should result in sys_modified
+    sys_modify_updated = deepcopy(sys_original)
+    PF.update_system!(sys_modify_updated, data_modified)
+    @test IS.compare_values(sys_modify_updated, sys_modified)
+end
diff --git a/test/test_psse_export.jl b/test/test_psse_export.jl
index abcc0604..9518e2b3 100644
--- a/test/test_psse_export.jl
+++ b/test/test_psse_export.jl
@@ -1,6 +1,3 @@
-const SYSTEM_REIMPORT_COMPARISON_TOLERANCE = 1e-10
-const POWERFLOW_COMPARISON_TOLERANCE = 1e-9
-
 test_psse_export_dir = joinpath(TEST_FILES_DIR, "test_psse_exports")  # at some point could move this to temp files
 isdir(test_psse_export_dir) && rm(test_psse_export_dir; recursive = true)
 
@@ -231,13 +228,6 @@ function compare_systems_wrapper(sys1::System, sys2::System, sys2_metadata = not
     return first_result && second_result
 end
 
-# TODO temporary hack, see https://github.com/NREL-Sienna/PowerFlows.jl/issues/39
-function PowerSystems.get_reactive_power_limits(gen::RenewableNonDispatch)
-    gen_pf = get_power_factor(gen)
-    gen_q = get_max_active_power(gen) * sqrt((1 / gen_pf^2) - 1)
-    return (min = 0.0, max = gen_q)
-end
-
 function test_power_flow(sys1::System, sys2::System)
     result1 = solve_powerflow(ACPowerFlow(), sys1)
     result2 = solve_powerflow(ACPowerFlow(), sys2)
diff --git a/test/test_utils/common.jl b/test/test_utils/common.jl
new file mode 100644
index 00000000..c8932249
--- /dev/null
+++ b/test/test_utils/common.jl
@@ -0,0 +1,36 @@
+const SYSTEM_REIMPORT_COMPARISON_TOLERANCE = 1e-10
+const POWERFLOW_COMPARISON_TOLERANCE = 1e-9
+
+powerflow_match_fn(
+    a::T,
+    b::T,
+) where {T <: Union{AbstractFloat, AbstractArray{<:AbstractFloat}}} =
+    isapprox(a, b; atol = POWERFLOW_COMPARISON_TOLERANCE) || IS.isequivalent(a, b)
+powerflow_match_fn(a, b) = IS.isequivalent(a, b)
+
+# TODO temporary hack, see https://github.com/NREL-Sienna/PowerFlows.jl/issues/39
+function PowerSystems.get_reactive_power_limits(gen::RenewableNonDispatch)
+    gen_pf = get_power_factor(gen)
+    gen_q = get_max_active_power(gen) * sqrt((1 / gen_pf^2) - 1)
+    return (min = 0.0, max = gen_q)
+end
+
+"Take RTS_GMLC_DA_sys and make some changes to it that are fully captured in the PowerFlowData(ACPowerFlow(), ...)"
+function modify_rts_system!(sys::System)
+    # For REF bus, voltage and angle are fixed; update active and reactive
+    ref_bus = get_component(ACBus, sys, "Arne")  # bus number 113
+    @assert get_bustype(ref_bus) == ACBusTypes.REF
+    # NOTE: we are not testing the correctness of _power_redistribution_ref here, it is used on both sides of the test
+    PF._power_redistribution_ref(sys, 2.4375, 0.1875, ref_bus)
+
+    # For PV bus, active and voltage are fixed; update reactive and angle
+
+    # For PQ bus, active and reactive are fixed; update voltage and angle
+
+end
+
+"Make the same changes to the PowerFlowData that modify_rts_system! makes to the System"
+function modify_rts_powerflow!(data::PowerFlowData)
+    data.bus_activepower_injection[data.bus_lookup[113]] -= 1.0
+    data.bus_reactivepower_injection[data.bus_lookup[113]] -= 1.0
+end