diff --git a/Project.toml b/Project.toml
index 8ae07ff87..3aefc6f55 100644
--- a/Project.toml
+++ b/Project.toml
@@ -31,8 +31,8 @@ ForwardDiff = "~v0.10"
 InfrastructureSystems = "2"
 NLsolve = "4"
 PowerSystems = "4"
-PowerFlows = "^0.7"
-PowerNetworkMatrices = "^0.11"
+PowerFlows = "^0.9"
+PowerNetworkMatrices = "^0.12.1"
 PrettyTables = "1, 2"
 SciMLBase = "2"
 TimerOutputs = "~0.5"
diff --git a/docs/src/tutorials/tutorial_continuation_pf.md b/docs/src/tutorials/tutorial_continuation_pf.md
index 364e30ff9..2904532b7 100644
--- a/docs/src/tutorials/tutorial_continuation_pf.md
+++ b/docs/src/tutorials/tutorial_continuation_pf.md
@@ -77,7 +77,8 @@ for p in P_range
     q_power = power * tan(acos(load_pf))
     set_reactive_power!(load, q_power)
     # Run Power Flow
-    status = solve_ac_powerflow!(sys_static)
+    pf = ACPowerFlow()
+    status = solve_powerflow!(pf, sys_static)
     if !status
         # Finish the loop if the power flow fails
         print("Power Flow failed at p = $(power)")
@@ -143,7 +144,8 @@ for p in P_range
     q_power = power * tan(acos(load_pf))
     set_reactive_power!(load, q_power)
     # Run Power Flow
-    status = solve_ac_powerflow!(sys_static)
+    pf = ACPowerFlow()
+    status = solve_powerflow!(pf, sys_static)
     if !status
         # Finish the loop if the power flow fails
         print("Power Flow failed at p = $(power)")
diff --git a/src/base/nlsolve_wrapper.jl b/src/base/nlsolve_wrapper.jl
index 371e29bcf..c58778bf6 100644
--- a/src/base/nlsolve_wrapper.jl
+++ b/src/base/nlsolve_wrapper.jl
@@ -133,6 +133,7 @@ function _check_residual(
                 Generator = $gen_name, state = $state.
                Residual error is too large to continue")
         else
+            bus_count = get_bus_count(inputs)
             bus_no = ix > bus_count ? ix - bus_count : ix
             component = ix > bus_count ? "imag" : "real"
             error("The initial residual in the state located at $ix has a value of $val.
diff --git a/src/base/simulation_initialization.jl b/src/base/simulation_initialization.jl
index a37786cd8..daefc875a 100644
--- a/src/base/simulation_initialization.jl
+++ b/src/base/simulation_initialization.jl
@@ -11,7 +11,8 @@ function power_flow_solution!(
     sys::PSY.System,
     inputs::SimulationInputs,
 )
-    res = PF.solve_ac_powerflow!(sys)
+    pf = PF.ACPowerFlow()
+    res = PF.solve_powerflow!(pf, sys)
     if !res
         @error("PowerFlow failed to solve")
         return BUILD_FAILED
diff --git a/src/base/simulation_results.jl b/src/base/simulation_results.jl
index 3c9efe157..75ac4a5ef 100644
--- a/src/base/simulation_results.jl
+++ b/src/base/simulation_results.jl
@@ -70,10 +70,17 @@ function post_proc_state_series(
 )
     global_state_index = get_global_index(res)
     if !haskey(global_state_index, ref[1])
-        @error "$(keys(global_state_index))"
-        error("State $(ref[2]) device $(ref[1]) not found in the system. ")
+        available_devices = join(keys(global_state_index), ", ")
+        error("Device $(ref[1]) not found. Available devices: $available_devices")
     end
-    ix = get(global_state_index[ref[1]], ref[2], 0)
+
+    state_dict = global_state_index[ref[1]]
+    if !haskey(state_dict, ref[2])
+        available_states = join(keys(state_dict), ", ")
+        error("State $(ref[2]) not found in device $(ref[1]). Available states: $available_states")
+    end
+
+    ix = state_dict[ref[2]]
     return _post_proc_state_series(get_solution(res), ix, dt, unique_timestamps)
 end
 
diff --git a/src/post_processing/post_proc_generator.jl b/src/post_processing/post_proc_generator.jl
index 8884e084b..08f0aed65 100644
--- a/src/post_processing/post_proc_generator.jl
+++ b/src/post_processing/post_proc_generator.jl
@@ -1018,7 +1018,7 @@ function _mechanical_torque(
     name::String,
     res::SimulationResults,
     dt::Union{Nothing, Float64, Vector{Float64}},
-    unique_timestamps::Bool,
+    unique_timestamps::Bool = true,
 )
     # TODO: This will not plot correctly when changing P_ref in a callback
     ts, _ = _post_proc_state_series(res.solution, 1, dt, unique_timestamps)
@@ -1038,7 +1038,7 @@ function _mechanical_torque(
     name::String,
     res::SimulationResults,
     dt::Union{Nothing, Float64, Vector{Float64}},
-    unique_timestamps::Bool,
+    unique_timestamps::Bool = true,
 )
     # Get params
     Tc = PSY.get_Tc(tg)
@@ -1067,7 +1067,7 @@ function _mechanical_torque(
     name::String,
     res::SimulationResults,
     dt::Union{Nothing, Float64, Vector{Float64}},
-    unique_timestamps::Bool,
+    unique_timestamps::Bool = true,
 )
     # TODO: This will not plot correctly when changing P_ref in a callback
     # Get params
@@ -1097,7 +1097,7 @@ function _mechanical_torque(
     name::String,
     res::SimulationResults,
     dt::Union{Nothing, Float64, Vector{Float64}},
-    unique_timestamps::Bool,
+    unique_timestamps::Bool = true,
 )
     # TODO: This will not plot correctly when changing P_ref in a callback
     # Get params
@@ -1133,7 +1133,7 @@ function _mechanical_torque(
     name::String,
     res::SimulationResults,
     dt::Union{Nothing, Float64, Vector{Float64}},
-    unique_timestamps::Bool,
+    unique_timestamps::Bool = true,
 )
     # Get params
     D_turb = PSY.get_D_turb(tg)
@@ -1154,6 +1154,7 @@ function _mechanical_torque(
     name::String,
     res::SimulationResults,
     dt::Union{Nothing, Float64},
+    unique_timestamps::Bool = true,
 )
     ts, x_a3 = post_proc_state_series(res, (name, :x_a3), dt, unique_timestamps)
     return ts, x_a3
@@ -1167,7 +1168,7 @@ function _mechanical_torque(
     name::String,
     res::SimulationResults,
     dt::Union{Nothing, Float64, Vector{Float64}},
-    unique_timestamps::Bool,
+    unique_timestamps::Bool = true,
 )
     # Get params
     q_nl = PSY.get_q_nl(tg)
@@ -1194,6 +1195,7 @@ function _mechanical_torque(
     name::String,
     res::SimulationResults,
     dt::Union{Nothing, Float64, Vector{Float64}},
+    unique_timestamps::Bool = true,
 )
     # Get params
     D_turb = PSY.get_D_turb(tg)
@@ -1204,9 +1206,9 @@ function _mechanical_torque(
     setpoints = get_setpoints(res)
     ω_ref = setpoints[name]["ω_ref"]
     # Get state results
-    ts, x_g7 = post_proc_state_series(res, (name, :x_g7), dt)
-    _, x_g6 = post_proc_state_series(res, (name, :x_g6), dt)
-    _, ω = post_proc_state_series(res, (name, :ω), dt)
+    ts, x_g7 = post_proc_state_series(res, (name, :x_g7), dt, unique_timestamps)
+    _, x_g6 = post_proc_state_series(res, (name, :x_g6), dt, unique_timestamps)
+    _, ω = post_proc_state_series(res, (name, :ω), dt, unique_timestamps)
     Pe = similar(x_g7)
     for (ix, x7) in enumerate(x_g7)
         x6 = x_g6[ix]
@@ -1225,6 +1227,7 @@ function _mechanical_torque(
     name::String,
     res::SimulationResults,
     dt::Union{Nothing, Float64, Vector{Float64}},
+    unique_timestamps::Bool = true,
 )
     # Get params
     D = PSY.get_D(tg)
@@ -1235,9 +1238,9 @@ function _mechanical_torque(
     ω_ref = setpoints[name]["ω_ref"]
 
     # Get state results
-    ts, x_g7 = post_proc_state_series(res, (name, :x_g7), dt)
-    _, x_g6 = post_proc_state_series(res, (name, :x_g6), dt)
-    _, ω = post_proc_state_series(res, (name, :ω), dt)
+    ts, x_g7 = post_proc_state_series(res, (name, :x_g7), dt, unique_timestamps)
+    _, x_g6 = post_proc_state_series(res, (name, :x_g6), dt, unique_timestamps)
+    _, ω = post_proc_state_series(res, (name, :ω), dt, unique_timestamps)
     Pm = similar(x_g7)
 
     for (ix, x7) in enumerate(x_g7)
diff --git a/test/test_case_5shaft.jl b/test/test_case_5shaft.jl
index 9f5946ed4..f0911632b 100644
--- a/test/test_case_5shaft.jl
+++ b/test/test_case_5shaft.jl
@@ -11,7 +11,8 @@ The fault disconnects a circuit between buses 1 and 2, doubling its impedance.
 ##################################################
 
 threebus_sys = build_system(PSIDTestSystems, "psid_test_threebus_5shaft")
-solve_ac_powerflow!(threebus_sys)
+pf = ACPowerFlow()
+solve_powerflow!(pf, threebus_sys)
 Ybus_fault = get_ybus_fault_threebus_sys(threebus_sys)
 
 ##################################################
diff --git a/test/test_case_anderson.jl b/test/test_case_anderson.jl
index d6361d59a..87c57bfb1 100644
--- a/test/test_case_anderson.jl
+++ b/test/test_case_anderson.jl
@@ -10,7 +10,8 @@ and the generator located in bus 3.
 ##################################################
 
 threebus_sys = build_system(PSIDTestSystems, "psid_test_threebus_anderson")
-solve_ac_powerflow!(threebus_sys)
+pf = ACPowerFlow()
+solve_powerflow!(pf, threebus_sys)
 Ybus_fault = get_ybus_fault_threebus_sys(threebus_sys)
 
 ##################################################
diff --git a/test/test_case_marconato.jl b/test/test_case_marconato.jl
index 26e8cc66d..6d7b7f388 100644
--- a/test/test_case_marconato.jl
+++ b/test/test_case_marconato.jl
@@ -10,7 +10,8 @@ and the generator located in bus 3.
 ##################################################
 
 threebus_sys = build_system(PSIDTestSystems, "psid_test_threebus_marconato")
-solve_ac_powerflow!(threebus_sys)
+pf = ACPowerFlow()
+solve_powerflow!(pf, threebus_sys)
 Ybus_fault = get_ybus_fault_threebus_sys(threebus_sys)
 
 ##################################################
diff --git a/test/test_case_oneDoneQ.jl b/test/test_case_oneDoneQ.jl
index 32898d2ee..3daf6583b 100644
--- a/test/test_case_oneDoneQ.jl
+++ b/test/test_case_oneDoneQ.jl
@@ -10,7 +10,8 @@ and the generator located in bus 3.
 ##################################################
 
 threebus_sys = build_system(PSIDTestSystems, "psid_test_threebus_oneDoneQ")
-solve_ac_powerflow!(threebus_sys)
+pf = ACPowerFlow()
+solve_powerflow!(pf, threebus_sys)
 Ybus_fault = get_ybus_fault_threebus_sys(threebus_sys)
 
 ##################################################
diff --git a/test/test_case_simple_anderson.jl b/test/test_case_simple_anderson.jl
index a75fe5f29..7f0a72855 100644
--- a/test/test_case_simple_anderson.jl
+++ b/test/test_case_simple_anderson.jl
@@ -10,7 +10,8 @@ and the generator located in bus 3.
 ##################################################
 
 threebus_sys = build_system(PSIDTestSystems, "psid_test_threebus_simple_anderson")
-solve_ac_powerflow!(threebus_sys)
+pf = ACPowerFlow()
+solve_powerflow!(pf, threebus_sys)
 Ybus_fault = get_ybus_fault_threebus_sys(threebus_sys)
 
 ##################################################
diff --git a/test/test_case_simple_marconato.jl b/test/test_case_simple_marconato.jl
index 73138b5cd..1ec5aa408 100644
--- a/test/test_case_simple_marconato.jl
+++ b/test/test_case_simple_marconato.jl
@@ -10,7 +10,8 @@ and the generator located in bus 3.
 ##################################################
 
 threebus_sys = build_system(PSIDTestSystems, "psid_test_threebus_simple_marconato")
-solve_ac_powerflow!(threebus_sys)
+pf = ACPowerFlow()
+solve_powerflow!(pf, threebus_sys)
 Ybus_fault = get_ybus_fault_threebus_sys(threebus_sys)
 
 ##################################################