diff --git a/Project.toml b/Project.toml
index 45960a6..4a43da3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,7 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
 JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
+KLU = "ef3ab10e-7fda-4108-b977-705223b18434"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
@@ -22,6 +23,7 @@ DataStructures = "0.18"
 Dates = "1"
 InfrastructureSystems = "2"
 JSON3 = "1"
+KLU = "^0.6"
 LinearAlgebra = "1"
 Logging = "1"
 NLsolve = "4"
diff --git a/src/PowerFlowData.jl b/src/PowerFlowData.jl
index a9d68bd..af1b5bc 100644
--- a/src/PowerFlowData.jl
+++ b/src/PowerFlowData.jl
@@ -180,7 +180,7 @@ NOTE: use it for AC power flow computations.
 WARNING: functions for the evaluation of the multi-period AC PF still to be implemented.
 """
 function PowerFlowData(
-    ::ACPowerFlow,
+    ::ACPowerFlow{<:ACPowerFlowSolverType},
     sys::PSY.System;
     time_steps::Int = 1,
     timestep_names::Vector{String} = String[],
@@ -196,6 +196,8 @@ function PowerFlowData(
         end
     end
 
+    timestep_map = Dict(zip([i for i in 1:time_steps], timestep_names))
+
     # get data for calculations
     power_network_matrix = PNM.Ybus(sys; check_connectivity = check_connectivity)
 
@@ -457,7 +459,11 @@ Create an appropriate `PowerFlowContainer` for the given `PowerFlowEvaluationMod
 """
 function make_power_flow_container end
 
-make_power_flow_container(pfem::ACPowerFlow, sys::PSY.System; kwargs...) =
+make_power_flow_container(
+    pfem::ACPowerFlow{<:ACPowerFlowSolverType},
+    sys::PSY.System;
+    kwargs...,
+) =
     PowerFlowData(pfem, sys; kwargs...)
 
 make_power_flow_container(pfem::DCPowerFlow, sys::PSY.System; kwargs...) =
diff --git a/src/PowerFlows.jl b/src/PowerFlows.jl
index 83562e2..78cbeeb 100644
--- a/src/PowerFlows.jl
+++ b/src/PowerFlows.jl
@@ -1,10 +1,13 @@
 module PowerFlows
 
 export solve_powerflow
-export solve_ac_powerflow!
+export solve_powerflow!
 export PowerFlowData
 export DCPowerFlow
+export NLSolveACPowerFlow
+export KLUACPowerFlow
 export ACPowerFlow
+export ACPowerFlowSolverType
 export PTDFDCPowerFlow
 export vPTDFDCPowerFlow
 export PSSEExportPowerFlow
@@ -20,10 +23,11 @@ import PowerSystems
 import PowerSystems: System, with_units_base
 import LinearAlgebra
 import NLsolve
+import KLU
 import SparseArrays
 import InfrastructureSystems
 import PowerNetworkMatrices
-import SparseArrays: SparseMatrixCSC
+import SparseArrays: SparseMatrixCSC, sparse
 import JSON3
 import DataStructures: OrderedDict
 import Dates
@@ -40,6 +44,7 @@ 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("newton_ac_powerflow.jl")
+include("nlsolve_powerflow.jl")
 include("post_processing.jl")
 end
diff --git a/src/ac_power_flow.jl b/src/ac_power_flow.jl
index 617683c..ee7b2f2 100644
--- a/src/ac_power_flow.jl
+++ b/src/ac_power_flow.jl
@@ -55,16 +55,17 @@ function _calculate_x0(n::Int,
     return x0
 end
 
-function PolarPowerFlow(data::ACPowerFlowData)
-    time_step = 1  # TODO placeholder time_step
+function PolarPowerFlow(data::ACPowerFlowData; time_step::Int64 = 1)
     n_buses = first(size(data.bus_type))
     P_net = zeros(n_buses)
     Q_net = zeros(n_buses)
     for ix in 1:n_buses
         P_net[ix] =
-            data.bus_activepower_injection[ix] - data.bus_activepower_withdrawals[ix]
+            data.bus_activepower_injection[ix, time_step] -
+            data.bus_activepower_withdrawals[ix, time_step]
         Q_net[ix] =
-            data.bus_reactivepower_injection[ix] - data.bus_reactivepower_withdrawals[ix]
+            data.bus_reactivepower_injection[ix, time_step] -
+            data.bus_reactivepower_withdrawals[ix, time_step]
     end
     x0 = _calculate_x0(1,
         data.bus_type[:, time_step],
diff --git a/src/common.jl b/src/common.jl
index 9ea6a41..cd2d1fb 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -159,7 +159,7 @@ function make_powerflowdata(
 )
     bus_type = Vector{PSY.ACBusTypes}(undef, n_buses)
     bus_angles = zeros(Float64, n_buses)
-    bus_magnitude = zeros(Float64, n_buses)
+    bus_magnitude = ones(Float64, n_buses)
 
     _initialize_bus_data!(
         bus_type,
@@ -190,8 +190,9 @@ function make_powerflowdata(
     )
 
     # Shapes to reuse
-    zeros_bus_time = () -> zeros(n_buses, time_steps)
-    zeros_branch_time = () -> zeros(n_branches, time_steps)
+    zeros_bus_time = () -> zeros(Float64, n_buses, time_steps)
+    ones_bus_time = () -> ones(Float64, n_buses, time_steps)
+    zeros_branch_time = () -> zeros(Float64, n_branches, time_steps)
 
     # Define fields as matrices whose number of columns is equal to the number of time_steps
     bus_activepower_injection_1 = zeros_bus_time()
@@ -199,7 +200,7 @@ function make_powerflowdata(
     bus_activepower_withdrawals_1 = zeros_bus_time()
     bus_reactivepower_withdrawals_1 = zeros_bus_time()
     bus_reactivepower_bounds_1 = Matrix{Vector{Float64}}(undef, n_buses, time_steps)
-    bus_magnitude_1 = zeros_bus_time()
+    bus_magnitude_1 = ones_bus_time()
     bus_angles_1 = zeros_bus_time()
 
     # Initial values related to first timestep allocated in the first column
diff --git a/src/definitions.jl b/src/definitions.jl
index d8ecad1..36e53b6 100644
--- a/src/definitions.jl
+++ b/src/definitions.jl
@@ -10,4 +10,7 @@ const DEFAULT_MAX_REDISTRIBUTION_ITERATIONS = 10
 
 const ISAPPROX_ZERO_TOLERANCE = 1e-6
 
+const DEFAULT_NR_MAX_ITER::Int64 = 30   # default maxIter for the NR power flow
+const DEFAULT_NR_TOL::Float64 = 1e-9 # default tolerance for the NR power flow
+
 const AC_PF_KW = []
diff --git a/src/newton_ac_powerflow.jl b/src/newton_ac_powerflow.jl
new file mode 100644
index 0000000..570732f
--- /dev/null
+++ b/src/newton_ac_powerflow.jl
@@ -0,0 +1,722 @@
+"""
+Solves a the power flow into the system and writes the solution into the relevant structs.
+Updates generators active and reactive power setpoints and branches active and reactive
+power flows (calculated in the From - To direction) (see
+[`flow_val`](@ref))
+
+Supports passing NLsolve kwargs in the args. By default shows the solver trace.
+
+Arguments available for `nlsolve`:
+
+  - `get_connectivity::Bool`: Checks if the network is connected. Default true
+  - `method` : See NLSolve.jl documentation for available solvers
+  - `xtol`: norm difference in `x` between two successive iterates under which
+    convergence is declared. Default: `0.0`.
+  - `ftol`: infinite norm of residuals under which convergence is declared.
+    Default: `1e-8`.
+  - `iterations`: maximum number of iterations. Default: `1_000`.
+  - `store_trace`: should a trace of the optimization algorithm's state be
+    stored? Default: `false`.
+  - `show_trace`: should a trace of the optimization algorithm's state be shown
+    on `STDOUT`? Default: `false`.
+  - `extended_trace`: should additifonal algorithm internals be added to the state
+    trace? Default: `false`.
+
+## Examples
+
+```julia
+solve_ac_powerflow!(sys)
+# Passing NLsolve arguments
+solve_ac_powerflow!(sys, method=:newton)
+```
+"""
+function solve_powerflow!(
+    pf::ACPowerFlow{<:ACPowerFlowSolverType},
+    system::PSY.System;
+    time_step::Int64 = 1,
+    kwargs...,
+)
+    #Save per-unit flag
+    settings_unit_cache = deepcopy(system.units_settings.unit_system)
+    #Work in System per unit
+    PSY.set_units_base_system!(system, "SYSTEM_BASE")
+
+    data = PowerFlowData(
+        pf,
+        system;
+        check_connectivity = get(kwargs, :check_connectivity, true),
+    )
+
+    converged, V, Sbus_result = _ac_powereflow(data, pf; time_step = time_step, kwargs...)
+    x = _calc_x(data, V, Sbus_result)
+
+    if converged
+        write_powerflow_solution!(
+            system,
+            x,
+            data,
+            get(kwargs, :maxIter, DEFAULT_NR_MAX_ITER),
+        )
+        @info("PowerFlow solve converged, the results have been stored in the system")
+    else
+        @error("The powerflow solver returned convergence = $(converged)")
+    end
+
+    #Restore original per unit base
+    PSY.set_units_base_system!(system, settings_unit_cache)
+
+    return converged
+end
+
+"""
+Similar to solve_powerflow!(sys) but does not update the system struct with results.
+Returns the results in a dictionary of dataframes.
+
+## Examples
+
+```julia
+res = solve_powerflow(sys)
+# Passing NLsolve arguments
+res = solve_powerflow(sys, method=:newton)
+```
+"""
+function solve_powerflow(
+    pf::ACPowerFlow{<:ACPowerFlowSolverType},
+    system::PSY.System;
+    kwargs...,
+)
+    #Save per-unit flag
+    settings_unit_cache = deepcopy(system.units_settings.unit_system)
+    #Work in System per unit
+    PSY.set_units_base_system!(system, "SYSTEM_BASE")
+
+    data = PowerFlowData(
+        pf,
+        system;
+        check_connectivity = get(kwargs, :check_connectivity, true),
+    )
+
+    converged, V, Sbus_result = _ac_powereflow(data, pf; kwargs...)
+    x = _calc_x(data, V, Sbus_result)
+
+    if converged
+        @info("PowerFlow solve converged, the results are exported in DataFrames")
+        df_results = write_results(pf, system, data, x)
+    else
+        df_results = missing
+        @error("The powerflow solver returned convergence = $(converged)")
+    end
+
+    #Restore original per unit base
+    PSY.set_units_base_system!(system, settings_unit_cache)
+
+    return df_results
+end
+
+# Multiperiod power flow - work in progress
+function solve_powerflow(
+    pf::ACPowerFlow{<:ACPowerFlowSolverType},
+    data::ACPowerFlowData,
+    system::PSY.System;
+    kwargs...,
+)
+    #Save per-unit flag
+    settings_unit_cache = deepcopy(system.units_settings.unit_system)
+    #Work in System per unit
+    PSY.set_units_base_system!(system, "SYSTEM_BASE")
+
+    sorted_time_steps = sort(collect(keys(data.timestep_map)))
+    # preallocate results
+    ts_converged = zeros(Bool, 1, length(sorted_time_steps))
+    ts_V = zeros(Complex{Float64}, length(data.bus_type[:, 1]), length(sorted_time_steps))
+    ts_S = zeros(Complex{Float64}, length(data.bus_type[:, 1]), length(sorted_time_steps))
+
+    results = Dict()
+
+    for t in sorted_time_steps
+        converged, V, Sbus_result =
+            _ac_powereflow(data, pf; time_step = t, kwargs...)
+        ts_converged[1, t] = converged
+        ts_V[:, t] .= V
+        ts_S[:, t] .= Sbus_result
+
+        # temporary implementation that will need to be improved:
+        if converged
+            x = _calc_x(data, V, Sbus_result)
+            results[data.timestep_map[t]] = write_results(pf, system, data, x)
+        else
+            results[data.timestep_map[t]] = missing
+        end
+
+        #todo: implement write_results for multiperiod power flow
+
+        # if converged
+        #     @info("PowerFlow solve converged, the results are exported in DataFrames")
+        #     df_results = write_results(pf, system, data, x)
+        # else
+        #     df_results = missing
+        #     @error("The powerflow solver returned convergence = $(converged)")
+        # end
+    end
+
+    #Restore original per unit base
+    PSY.set_units_base_system!(system, settings_unit_cache)
+
+    # todo:
+    # return df_results
+
+    return results
+end
+
+# Multiperiod power flow - work in progress
+function solve_powerflow!(
+    data::ACPowerFlowData;
+    kwargs...,
+)
+    pf = ACPowerFlow()  # todo: somehow store in data which PF to use (see issue #50)
+
+    sorted_time_steps = sort(collect(keys(data.timestep_map)))
+    # preallocate results
+    ts_converged = zeros(Bool, 1, length(sorted_time_steps))
+    ts_V = zeros(Complex{Float64}, length(data.bus_type[:, 1]), length(sorted_time_steps))
+    ts_S = zeros(Complex{Float64}, length(data.bus_type[:, 1]), length(sorted_time_steps))
+
+    Yft = data.power_network_matrix.data_ft
+    Ytf = data.power_network_matrix.data_tf
+
+    fb = data.power_network_matrix.fb
+    tb = data.power_network_matrix.tb
+
+    for t in sorted_time_steps
+        converged, V, Sbus_result =
+            _ac_powereflow(data, pf; time_step = t, kwargs...)
+        ts_converged[1, t] = converged
+        ts_V[:, t] .= V
+        ts_S[:, t] .= Sbus_result
+
+        ref = findall(
+            x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF,
+            data.bus_type[:, t],
+        )
+        pv = findall(
+            x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV,
+            data.bus_type[:, t],
+        )
+        pq = findall(
+            x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ,
+            data.bus_type[:, t],
+        )
+
+        # temporary implementation that will need to be improved:
+        if converged
+            # write results for REF
+            data.bus_activepower_injection[ref, t] .=
+                real.(Sbus_result[ref]) .+ data.bus_activepower_withdrawals[ref, t]
+            data.bus_reactivepower_injection[ref, t] .=
+                imag.(Sbus_result[ref]) .+ data.bus_reactivepower_withdrawals[ref, t]
+            # write Q results for PV
+            data.bus_reactivepower_injection[pv, t] .=
+                imag.(Sbus_result[pv]) .+ data.bus_reactivepower_withdrawals[pv, t]
+            # results for PQ buses do not need to be updated -> already consistent with inputs
+
+            # write bus bus_types
+            # todo
+
+            # write voltage results
+            data.bus_magnitude[pq, t] .= abs.(V[pq])
+            data.bus_angles[pq, t] .= angle.(V[pq])
+            data.bus_angles[pv, t] .= angle.(V[pv])
+
+        else
+            # todo
+            1 + 2
+        end
+    end
+
+    # write branch flows
+    Sft = ts_V[fb, :] .* conj.(Yft * ts_V)
+    Stf = ts_V[tb, :] .* conj.(Ytf * ts_V)
+
+    data.branch_activepower_flow_from_to .= real.(Sft)
+    data.branch_reactivepower_flow_from_to .= imag.(Sft)
+    data.branch_activepower_flow_to_from .= real.(Stf)
+    data.branch_reactivepower_flow_to_from .= imag.(Stf)
+
+    # todo:
+    # return df_results
+
+    return
+end
+
+function _ac_powereflow(
+    data::ACPowerFlowData,
+    pf::ACPowerFlow{<:ACPowerFlowSolverType};
+    time_step::Int64 = 1,
+    kwargs...,
+)
+    check_reactive_power_limits = get(kwargs, :check_reactive_power_limits, false)
+
+    for _ in 1:MAX_REACTIVE_POWER_ITERATIONS
+        converged, V, Sbus_result =
+            _newton_powerflow(pf, data; time_step = time_step, kwargs...)
+        if !converged || !check_reactive_power_limits ||
+           _check_q_limit_bounds!(data, Sbus_result, time_step)
+            return (converged, V, Sbus_result)
+        end
+    end
+
+    @error("could not enforce reactive power limits after $MAX_REACTIVE_POWER_ITERATIONS")
+    return (converged, V, Sbus_result)
+end
+
+function _check_q_limit_bounds!(
+    data::ACPowerFlowData,
+    Sbus_result::Vector{Complex{Float64}},
+    time_step::Int64,
+)
+    bus_names = data.power_network_matrix.axes[1]
+    within_limits = true
+    for (ix, bt) in enumerate(data.bus_type[:, time_step])
+        if bt == PSY.ACBusTypes.PV
+            Q_gen = imag(Sbus_result[ix])
+        else
+            continue
+        end
+
+        if Q_gen <= data.bus_reactivepower_bounds[ix, time_step][1]
+            @info "Bus $(bus_names[ix]) changed to PSY.ACBusTypes.PQ"
+            within_limits = false
+            data.bus_type[ix, time_step] = PSY.ACBusTypes.PQ
+            data.bus_reactivepower_injection[ix, time_step] =
+                data.bus_reactivepower_bounds[ix, time_step][1]
+        elseif Q_gen >= data.bus_reactivepower_bounds[ix, time_step][2]
+            @info "Bus $(bus_names[ix]) changed to PSY.ACBusTypes.PQ"
+            within_limits = false
+            data.bus_type[ix, time_step] = PSY.ACBusTypes.PQ
+            data.bus_reactivepower_injection[ix, time_step] =
+                data.bus_reactivepower_bounds[ix, time_step][2]
+        else
+            @debug "Within Limits"
+        end
+    end
+    return within_limits
+end
+
+function _solve_powerflow!(
+    pf::ACPowerFlow{<:ACPowerFlowSolverType},
+    data::ACPowerFlowData,
+    check_reactive_power_limits;
+    time_step::Int64 = 1,
+    kwargs...,
+)
+    for _ in 1:MAX_REACTIVE_POWER_ITERATIONS
+        converged, V, Sbus_result =
+            _newton_powerflow(pf, data; time_step = time_step, kwargs...)
+        if !converged || !check_reactive_power_limits ||
+           _check_q_limit_bounds!(data, Sbus_result, time_step)
+            return converged, V, Sbus_result
+        end
+    end
+    # todo: throw error? set converged to false?
+    @error("could not enforce reactive power limits after $MAX_REACTIVE_POWER_ITERATIONS")
+    return converged, V, Sbus_result
+end
+
+function _update_V!(dx::Vector{Float64}, V::Vector{Complex{Float64}}, Vm::Vector{Float64},
+    Va::Vector{Float64},
+    pv::Vector{Int64}, pq::Vector{Int64}, dx_Va_pv::Vector{Int64}, dx_Va_pq::Vector{Int64},
+    dx_Vm_pq::Vector{Int64})
+    for (i, j) in zip(pv, dx_Va_pv)
+        Va[i] -= dx[j]
+    end
+
+    for (i, j) in zip(pq, dx_Va_pq)
+        Va[i] -= dx[j]
+    end
+
+    for (i, j) in zip(pq, dx_Vm_pq)
+        Vm[i] -= dx[j]
+    end
+
+    V .= Vm .* exp.(1im .* Va)
+
+    Vm .= abs.(V)
+    Va .= angle.(V)
+    return
+end
+
+function _update_F!(F::Vector{Float64}, Sbus_result::Vector{Complex{Float64}},
+    mis::Vector{Complex{Float64}},
+    dx_Va_pv::Vector{Int64}, dx_Va_pq::Vector{Int64}, dx_Vm_pq::Vector{Int64},
+    V::Vector{Complex{Float64}}, Ybus::SparseMatrixCSC{Complex{Float64}, Int64},
+    Sbus::Vector{Complex{Float64}},
+    pv::Vector{Int64}, pq::Vector{Int64})
+
+    #mis .= V .* conj.(Ybus * V) .- Sbus
+    LinearAlgebra.mul!(Sbus_result, Ybus, V)
+    Sbus_result .= V .* conj(Sbus_result)
+    mis .= Sbus_result .- Sbus
+
+    for (i, j) in zip(dx_Va_pv, pv)
+        F[i] = real(mis[j])
+    end
+
+    for (i, j) in zip(dx_Va_pq, pq)
+        F[i] = real(mis[j])
+    end
+
+    for (i, j) in zip(dx_Vm_pq, pq)
+        F[i] = imag(mis[j])
+    end
+    return
+end
+
+function _update_dSbus_dV!(rows::Vector{Int64}, cols::Vector{Int64},
+    V::Vector{Complex{Float64}}, Ybus::SparseMatrixCSC{Complex{Float64}, Int64},
+    diagV::LinearAlgebra.Diagonal{Complex{Float64}, Vector{Complex{Float64}}},
+    diagVnorm::LinearAlgebra.Diagonal{Complex{Float64}, Vector{Complex{Float64}}},
+    diagIbus::LinearAlgebra.Diagonal{Complex{Float64}, Vector{Complex{Float64}}},
+    diagIbus_diag::Vector{Complex{Float64}},
+    dSbus_dVa::SparseMatrixCSC{Complex{Float64}, Int64},
+    dSbus_dVm::SparseMatrixCSC{Complex{Float64}, Int64},
+    r_dSbus_dVa::SparseMatrixCSC{Float64, Int64},
+    r_dSbus_dVm::SparseMatrixCSC{Float64, Int64},
+    i_dSbus_dVa::SparseMatrixCSC{Float64, Int64},
+    i_dSbus_dVm::SparseMatrixCSC{Float64, Int64},
+    Ybus_diagVnorm::SparseMatrixCSC{Complex{Float64}, Int64},
+    conj_Ybus_diagVnorm::SparseMatrixCSC{Complex{Float64}, Int64},
+    diagV_conj_Ybus_diagVnorm::SparseMatrixCSC{Complex{Float64}, Int64},
+    conj_diagIbus::LinearAlgebra.Diagonal{Complex{Float64}, Vector{Complex{Float64}}},
+    conj_diagIbus_diagVnorm::LinearAlgebra.Diagonal{
+        Complex{Float64},
+        Vector{Complex{Float64}},
+    },
+    Ybus_diagV::SparseMatrixCSC{Complex{Float64}, Int64},
+    conj_Ybus_diagV::SparseMatrixCSC{Complex{Float64}, Int64})
+    for i in eachindex(V)
+        diagV[i, i] = V[i]
+        diagVnorm[i, i] = V[i] / abs(V[i])
+    end
+
+    # manually calculate the diagIbus matrix
+    LinearAlgebra.mul!(diagIbus_diag, Ybus, V)
+    for i in eachindex(V)
+        diagIbus[i, i] = diagIbus_diag[i]
+    end
+
+    # use the available matrices temporarily to calculate the dSbus_dV matrices
+    # original formula:
+    # dSbus_dVm .= diagV * conj.(Ybus * diagVnorm) + conj.(diagIbus) * diagVnorm
+    # non-allocating version:
+
+    LinearAlgebra.mul!(Ybus_diagVnorm, Ybus, diagVnorm)
+    conj_Ybus_diagVnorm .= conj.(Ybus_diagVnorm)
+    LinearAlgebra.mul!(diagV_conj_Ybus_diagVnorm, diagV, conj_Ybus_diagVnorm)
+    conj_diagIbus .= conj.(diagIbus)
+    LinearAlgebra.mul!(conj_diagIbus_diagVnorm, conj_diagIbus, diagVnorm)
+
+    dSbus_dVm .= diagV_conj_Ybus_diagVnorm
+    LinearAlgebra.axpy!(1, conj_diagIbus_diagVnorm, dSbus_dVm)
+
+    # original formula:
+    # dSbus_dVa .= 1im * diagV * conj.(diagIbus - Ybus * diagV)
+    # non-allocating version:
+    LinearAlgebra.mul!(Ybus_diagV, Ybus, diagV)
+    # Take the conjugate of the result (conj(Ybus * diagV)); conj_diagIbus is already available
+    conj_Ybus_diagV .= conj.(Ybus_diagV)
+
+    # write the result of conj.(diagIbus - Ybus * diagV) in conj_Ybus_diagV
+    LinearAlgebra.axpby!(1, conj_diagIbus, -1, conj_Ybus_diagV)
+
+    # Multiply the result by diagV
+    LinearAlgebra.mul!(dSbus_dVa, diagV, conj_Ybus_diagV)
+
+    # Now multiply by 1im to get the final result
+    #dSbus_dVa .*= 1im
+    LinearAlgebra.mul!(dSbus_dVa, dSbus_dVa, 1im)
+
+    # this loop is slower so we should use vectorize assignments below
+    # for c in cols
+    #     for r in rows
+    #         r_dSbus_dVa[r, c] = real(dSbus_dVa[r, c])
+    #         i_dSbus_dVa[r, c] = imag(dSbus_dVa[r, c])
+    #         r_dSbus_dVm[r, c] = real(dSbus_dVm[r, c])
+    #         i_dSbus_dVm[r, c] = imag(dSbus_dVm[r, c])
+    #     end
+    # end
+
+    # sometimes can allocate so we have to use the for loop above
+    r_dSbus_dVa .= real.(dSbus_dVa)
+    r_dSbus_dVm .= real.(dSbus_dVm)
+    i_dSbus_dVa .= imag.(dSbus_dVa)
+    i_dSbus_dVm .= imag.(dSbus_dVm)
+    return
+end
+
+function _update_submatrix!(
+    A::SparseMatrixCSC,
+    B::SparseMatrixCSC,
+    rows_A::Vector{Int64},
+    cols_A::Vector{Int64},
+    rows_B::Vector{Int64},
+    cols_B::Vector{Int64},
+)
+    for idj in eachindex(cols_A)
+        for idi in eachindex(rows_A)
+            A[rows_A[idi], cols_A[idj]] = B[rows_B[idi], cols_B[idj]]
+        end
+    end
+    return
+end
+
+function _update_J!(J::SparseMatrixCSC,
+    r_dSbus_dVa::SparseMatrixCSC,
+    r_dSbus_dVm::SparseMatrixCSC,
+    i_dSbus_dVa::SparseMatrixCSC,
+    i_dSbus_dVm::SparseMatrixCSC,
+    pvpq::Vector{Int64},
+    pq::Vector{Int64},
+    j_pvpq::Vector{Int64},
+    j_pq::Vector{Int64},
+)
+    _update_submatrix!(J, r_dSbus_dVa, j_pvpq, j_pvpq, pvpq, pvpq)
+    _update_submatrix!(J, r_dSbus_dVm, j_pvpq, j_pq, pvpq, pq)
+    _update_submatrix!(J, i_dSbus_dVa, j_pq, j_pvpq, pq, pvpq)
+    _update_submatrix!(J, i_dSbus_dVm, j_pq, j_pq, pq, pq)
+    return
+end
+
+function _calc_x(
+    data::ACPowerFlowData,
+    V::Vector{Complex{Float64}},
+    Sbus_result::Vector{Complex{Float64}};
+    time_step::Int64 = 1,
+)
+    Vm = abs.(V)
+    Va = angle.(V)
+    n_buses = length(V)
+    # mock the expected x format, where the values depend on the type of the bus:
+    x = zeros(Float64, 2 * n_buses)
+    for (ix, bt) in enumerate(data.bus_type[:, time_step])
+        if bt == PSY.ACBusTypes.REF
+            # When bustype == REFERENCE PSY.Bus, state variables are Active and Reactive Power Generated
+            x[2 * ix - 1] = real(Sbus_result[ix]) + data.bus_activepower_withdrawals[ix]
+            x[2 * ix] = imag(Sbus_result[ix]) + data.bus_reactivepower_withdrawals[ix]
+        elseif bt == PSY.ACBusTypes.PV
+            # When bustype == PV PSY.Bus, state variables are Reactive Power Generated and Voltage Angle
+            x[2 * ix - 1] = imag(Sbus_result[ix]) + data.bus_reactivepower_withdrawals[ix]
+            x[2 * ix] = Va[ix]
+        elseif bt == PSY.ACBusTypes.PQ
+            # When bustype == PQ PSY.Bus, state variables are Voltage Magnitude and Voltage Angle
+            x[2 * ix - 1] = Vm[ix]
+            x[2 * ix] = Va[ix]
+        end
+    end
+    return x
+end
+
+function _preallocate_J(
+    rows::Vector{Int64},
+    cols::Vector{Int64},
+    pvpq::Vector{Int64},
+    pq::Vector{Int64},
+)
+    J_block = sparse(rows, cols, Float64(0), maximum(rows), maximum(cols))
+    J = [J_block[pvpq, pvpq] J_block[pvpq, pq]; J_block[pq, pvpq] J_block[pq, pq]]
+    return J
+end
+
+function _newton_powerflow(
+    pf::ACPowerFlow{KLUACPowerFlow},
+    data::ACPowerFlowData;
+    time_step::Int64 = 1,
+    kwargs...,
+)
+    # Fetch maxIter and tol from kwargs, or use defaults if not provided
+    maxIter = get(kwargs, :maxIter, DEFAULT_NR_MAX_ITER)
+    tol = get(kwargs, :tol, DEFAULT_NR_TOL)
+    i = 0
+
+    Ybus = data.power_network_matrix.data
+
+    # Find indices for each bus type
+    ref = findall(
+        x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF,
+        data.bus_type[:, time_step],
+    )
+    pv = findall(
+        x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV,
+        data.bus_type[:, time_step],
+    )
+    pq = findall(
+        x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ,
+        data.bus_type[:, time_step],
+    )
+    pvpq = [pv; pq]
+
+    # nref = length(ref)
+    npv = length(pv)
+    npq = length(pq)
+    npvpq = npv + npq
+
+    Vm = data.bus_magnitude[:, time_step]
+    # prevent unfeasible starting values for Vm; for pv and ref buses we cannot do this:
+    Vm[pq] .= clamp.(Vm[pq], 0.9, 1.1)
+    Va = data.bus_angles[:, time_step]
+    V = zeros(Complex{Float64}, length(Vm))
+    V .= Vm .* exp.(1im .* Va)
+
+    # early return if only ref buses present - no need to solve the power flow
+    converged = npvpq == 0
+    if converged
+        # if only ref buses present, we do not need to enter the power flow loop
+        Sbus_result = V .* conj(Ybus * V)
+        return (converged, V, Sbus_result)
+    end
+
+    # we need to define lookups for mappings of pv, pq buses onto the internal J indexing
+    pvpq_lookup = zeros(Int64, maximum([ref; pvpq]) + 1)
+    pvpq_lookup[pvpq] .= 1:npvpq
+    pq_lookup = zeros(Int64, maximum([ref; pvpq]) + 1)
+    pq_lookup[pq] .= 1:npq
+
+    # define the internal J indexing using the lookup arrays
+    j_pvpq = pvpq_lookup[pvpq]
+    j_pq = npvpq .+ pq_lookup[pq]
+
+    # indices for updating of V
+    dx_Va_pv = Vector{Int64}([1:npv...])
+    dx_Va_pq = Vector{Int64}([(npv + 1):(npv + npq)...])
+    dx_Vm_pq = Vector{Int64}([(npv + npq + 1):(npv + 2 * npq)...])
+
+    Sbus =
+        data.bus_activepower_injection[:, time_step] -
+        data.bus_activepower_withdrawals[:, time_step] +
+        1im * (
+            data.bus_reactivepower_injection[:, time_step] -
+            data.bus_reactivepower_withdrawals[:, time_step]
+        )
+
+    # Pre-allocate mis and F and create views for the respective real and imaginary sections of the arrays:
+    mis = zeros(Complex{Float64}, length(V))
+    Sbus_result = zeros(Complex{Float64}, length(V))
+    F = zeros(Float64, npvpq + npq)
+    mis .= V .* conj.(Ybus * V) .- Sbus
+    F .= [real(mis[pvpq]); imag(mis[pq])]
+
+    # preallocate Jacobian matrix and arrays for calculating dSbus_dVa, dSbus_dVm
+    rows, cols = SparseArrays.findnz(Ybus)
+
+    #diagV = sparse(1:n_buses, 1:n_buses, V)
+    diagV = LinearAlgebra.Diagonal(V)
+    diagIbus_diag = zeros(Complex{Float64}, size(V, 1))
+    diagIbus = LinearAlgebra.Diagonal(diagIbus_diag)
+
+    diagVnorm = LinearAlgebra.Diagonal(V ./ abs.(V))
+
+    Ybus_diagVnorm = sparse(rows, cols, Complex{Float64}(0))
+    conj_Ybus_diagVnorm = sparse(rows, cols, Complex{Float64}(0))
+    diagV_conj_Ybus_diagVnorm = sparse(rows, cols, Complex{Float64}(0))
+    conj_diagIbus = conj.(diagIbus)
+    conj_diagIbus_diagVnorm = conj.(diagIbus)
+    Ybus_diagV = sparse(rows, cols, Complex{Float64}(0))
+    conj_Ybus_diagV = sparse(rows, cols, Complex{Float64}(0))
+
+    dSbus_dVm = sparse(rows, cols, Complex{Float64}(0))
+    dSbus_dVa = sparse(rows, cols, Complex{Float64}(0))
+    r_dSbus_dVa = sparse(rows, cols, Float64(0))
+    r_dSbus_dVm = sparse(rows, cols, Float64(0))
+    i_dSbus_dVa = sparse(rows, cols, Float64(0))
+    i_dSbus_dVm = sparse(rows, cols, Float64(0))
+
+    # maybe use this in the future?
+    # pvpq_rows = pvpq_lookup[rows][pvpq_lookup[rows] .!= 0]
+    # pvpq_cols = pvpq_lookup[cols][pvpq_lookup[cols] .!= 0]
+    # pq_rows = pq_lookup[rows][pq_lookup[rows] .!= 0] 
+    # pq_cols = pq_lookup[cols][pq_lookup[cols] .!= 0]
+
+    J = _preallocate_J(rows, cols, pvpq, pq)
+
+    # preallocate the KLU factorization object - symbolic object only
+    colptr = KLU.decrement(J.colptr)
+    rowval = KLU.decrement(J.rowval)
+    n = size(J, 1)
+    factor_J = KLU.KLUFactorization(n, colptr, rowval, J.nzval)
+    KLU.klu_analyze!(factor_J)
+    rf = Ref(factor_J.common)
+    # factorization for the numeric object does not work here:
+    # factor_J._numeric = KLU.klu_l_factor(colptr, rowval, J.nzval, factor_J._symbolic, rf)
+
+    while i < maxIter && !converged
+        i += 1
+
+        _update_dSbus_dV!(rows, cols, V, Ybus, diagV, diagVnorm, diagIbus, diagIbus_diag,
+            dSbus_dVa, dSbus_dVm, r_dSbus_dVa, r_dSbus_dVm, i_dSbus_dVa, i_dSbus_dVm,
+            Ybus_diagVnorm, conj_Ybus_diagVnorm, diagV_conj_Ybus_diagVnorm,
+            conj_diagIbus, conj_diagIbus_diagVnorm, Ybus_diagV, conj_Ybus_diagV)
+
+        # todo: improve pvpq, pq, j_pvpq, j_pq (use more specific indices)
+        _update_J!(
+            J,
+            r_dSbus_dVa,
+            r_dSbus_dVm,
+            i_dSbus_dVa,
+            i_dSbus_dVm,
+            pvpq,
+            pq,
+            j_pvpq,
+            j_pq,
+        )
+
+        # Workaround for the issue with KLU.klu_l_factor
+        # background: KLU.klu_l_factor does not work properly with the preallocated J matrix with dummy values
+        # the workaround is to initialize the numeric object here in the loop once and then refactorize the matrix in the loop inplace
+        if i == 1
+            # works when J values are ok:
+            factor_J._numeric =
+                KLU.klu_l_factor(colptr, rowval, J.nzval, factor_J._symbolic, rf)
+        end
+
+        try
+            # factorize the numeric object of KLU inplace, while reusing the symbolic object
+            KLU.klu_l_refactor(
+                colptr,
+                rowval,
+                J.nzval,
+                factor_J._symbolic,
+                factor_J._numeric,
+                rf,
+            )
+
+            # solve inplace - the results are written to F, so that we must use F instead of dx for updating V
+            KLU.klu_l_solve(
+                factor_J._symbolic,
+                factor_J._numeric,
+                size(F, 1),
+                size(F, 2),
+                F,
+                rf,
+            )
+
+        catch e
+            @error("KLU factorization failed: $e")
+            return (converged, V, Sbus_result)
+        end
+
+        # KLU.solve! overwrites F with the solution instead of returning it as dx, so -F is used here to update V
+        _update_V!(F, V, Vm, Va, pv, pq, dx_Va_pv, dx_Va_pq, dx_Vm_pq)
+
+        # here F is mismatch again
+        _update_F!(F, Sbus_result, mis, dx_Va_pv, dx_Va_pq, dx_Vm_pq, V, Ybus, Sbus, pv, pq)
+
+        converged = LinearAlgebra.norm(F, Inf) < tol
+    end
+
+    if !converged
+        @error("The powerflow solver with KLU did not converge after $i iterations")
+    else
+        @info("The powerflow solver with KLU converged after $i iterations")
+    end
+
+    return (converged, V, Sbus_result)
+end
diff --git a/src/nlsolve_ac_powerflow.jl b/src/nlsolve_ac_powerflow.jl
deleted file mode 100644
index ffce769..0000000
--- a/src/nlsolve_ac_powerflow.jl
+++ /dev/null
@@ -1,166 +0,0 @@
-const _SOLVE_AC_POWERFLOW_KWARGS = Set([:check_reactive_power_limits, :check_connectivity])
-"""
-Solves a the power flow into the system and writes the solution into the relevant structs.
-Updates generators active and reactive power setpoints and branches active and reactive
-power flows (calculated in the From - To direction) (see
-[`flow_val`](@ref))
-
-Supports passing NLsolve kwargs in the args. By default shows the solver trace.
-
-Arguments available for `nlsolve`:
-
-  - `get_connectivity::Bool`: Checks if the network is connected. Default true
-  - `method` : See NLSolve.jl documentation for available solvers
-  - `xtol`: norm difference in `x` between two successive iterates under which
-    convergence is declared. Default: `0.0`.
-  - `ftol`: infinite norm of residuals under which convergence is declared.
-    Default: `1e-8`.
-  - `iterations`: maximum number of iterations. Default: `1_000`.
-  - `store_trace`: should a trace of the optimization algorithm's state be
-    stored? Default: `false`.
-  - `show_trace`: should a trace of the optimization algorithm's state be shown
-    on `STDOUT`? Default: `false`.
-  - `extended_trace`: should additifonal algorithm internals be added to the state
-    trace? Default: `false`.
-
-## Examples
-
-```julia
-solve_ac_powerflow!(sys)
-# Passing NLsolve arguments
-solve_ac_powerflow!(sys, method=:newton)
-```
-"""
-function solve_ac_powerflow!(system::PSY.System; kwargs...)
-    #Save per-unit flag
-    settings_unit_cache = deepcopy(system.units_settings.unit_system)
-    #Work in System per unit
-    PSY.set_units_base_system!(system, "SYSTEM_BASE")
-    check_reactive_power_limits = get(kwargs, :check_reactive_power_limits, false)
-    data = PowerFlowData(
-        ACPowerFlow(; check_reactive_power_limits = check_reactive_power_limits),
-        system;
-        check_connectivity = get(kwargs, :check_connectivity, true),
-    )
-    max_iterations = DEFAULT_MAX_REDISTRIBUTION_ITERATIONS
-    solver_kwargs = filter(p -> !(p.first in _SOLVE_AC_POWERFLOW_KWARGS), kwargs)
-    res = _solve_powerflow!(data, check_reactive_power_limits; solver_kwargs...)
-    if res.f_converged
-        write_powerflow_solution!(system, res.zero, data, max_iterations)
-        @info("PowerFlow solve converged, the results have been stored in the system")
-        #Restore original per unit base
-        PSY.set_units_base_system!(system, settings_unit_cache)
-        return res.f_converged
-    end
-    @error("The powerflow solver returned convergence = $(res.f_converged)")
-    PSY.set_units_base_system!(system, settings_unit_cache)
-    return res.f_converged
-end
-
-"""
-Similar to solve_powerflow!(sys) but does not update the system struct with results.
-Returns the results in a dictionary of dataframes.
-
-## Examples
-
-```julia
-res = solve_powerflow(sys)
-# Passing NLsolve arguments
-res = solve_powerflow(sys, method=:newton)
-```
-"""
-function solve_powerflow(
-    pf::ACPowerFlow,
-    system::PSY.System;
-    kwargs...,
-)
-    #Save per-unit flag
-    settings_unit_cache = deepcopy(system.units_settings.unit_system)
-    #Work in System per unit
-    PSY.set_units_base_system!(system, "SYSTEM_BASE")
-    data = PowerFlowData(
-        pf,
-        system;
-        check_connectivity = get(kwargs, :check_connectivity, true),
-    )
-
-    res = _solve_powerflow!(data, pf.check_reactive_power_limits; kwargs...)
-
-    if res.f_converged
-        @info("PowerFlow solve converged, the results are exported in DataFrames")
-        df_results = write_results(pf, system, data, res.zero)
-        #Restore original per unit base
-        PSY.set_units_base_system!(system, settings_unit_cache)
-        return df_results
-    end
-    @error("The powerflow solver returned convergence = $(res.f_converged)")
-    PSY.set_units_base_system!(system, settings_unit_cache)
-    return res.f_converged
-end
-
-function _check_q_limit_bounds!(
-    data::ACPowerFlowData,
-    zero::Vector{Float64},
-    time_step::Int,
-)
-    bus_names = data.power_network_matrix.axes[1]
-    within_limits = true
-    for (ix, b) in enumerate(data.bus_type)
-        if b == PSY.ACBusTypes.PV
-            Q_gen = zero[2 * ix - 1]
-        else
-            continue
-        end
-
-        if Q_gen <= data.bus_reactivepower_bounds[ix, time_step][1]
-            @info "Bus $(bus_names[ix]) changed to PSY.ACBusTypes.PQ"
-            within_limits = false
-            data.bus_type[ix, time_step] = PSY.ACBusTypes.PQ
-            data.bus_reactivepower_injection[ix, time_step] =
-                data.bus_reactivepower_bounds[ix, time_step][1]
-        elseif Q_gen >= data.bus_reactivepower_bounds[ix, time_step][2]
-            @info "Bus $(bus_names[ix]) changed to PSY.ACBusTypes.PQ"
-            within_limits = false
-            data.bus_type[ix, time_step] = PSY.ACBusTypes.PQ
-            data.bus_reactivepower_injection[ix, time_step] =
-                data.bus_reactivepower_bounds[ix, time_step][2]
-        else
-            @debug "Within Limits"
-        end
-    end
-    return within_limits
-end
-
-function _solve_powerflow!(
-    data::ACPowerFlowData,
-    check_reactive_power_limits;
-    time_step = 1,  # TODO placeholder time_step
-    nlsolve_kwargs...,
-)
-    if check_reactive_power_limits
-        for _ in 1:MAX_REACTIVE_POWER_ITERATIONS
-            res = _nlsolve_powerflow(data; nlsolve_kwargs...)
-            if res.f_converged
-                if _check_q_limit_bounds!(data, res.zero, time_step)
-                    return res
-                end
-            else
-                return res
-            end
-        end
-    else
-        return _nlsolve_powerflow(data; nlsolve_kwargs...)
-    end
-end
-
-function _nlsolve_powerflow(data::ACPowerFlowData; nlsolve_kwargs...)
-    pf = PolarPowerFlow(data)
-    J = PowerFlows.PolarPowerFlowJacobian(data, pf.x0)
-
-    df = NLsolve.OnceDifferentiable(pf, J, pf.x0, pf.residual, J.Jv)
-    res = NLsolve.nlsolve(df, pf.x0; nlsolve_kwargs...)
-    if !res.f_converged
-        @error("The powerflow solver returned convergence = $(res.f_converged)")
-    end
-    return res
-end
diff --git a/src/nlsolve_powerflow.jl b/src/nlsolve_powerflow.jl
new file mode 100644
index 0000000..c04e5fc
--- /dev/null
+++ b/src/nlsolve_powerflow.jl
@@ -0,0 +1,66 @@
+const _NLSOLVE_AC_POWERFLOW_KWARGS =
+    Set([:check_reactive_power_limits, :check_connectivity])
+function _newton_powerflow(
+    pf::ACPowerFlow{NLSolveACPowerFlow},
+    data::ACPowerFlowData;
+    time_step::Int64 = 1,  # not implemented for NLSolve and not used
+    nlsolve_kwargs...,
+)
+    if time_step != 1
+        throw(
+            ArgumentError(
+                "Multiperiod power flow not implemented for NLSolve AC power flow",
+            ),
+        )
+    end
+
+    nlsolve_solver_kwargs =
+        filter(p -> !(p.first in _NLSOLVE_AC_POWERFLOW_KWARGS), nlsolve_kwargs)
+
+    pf = PolarPowerFlow(data)
+    J = PowerFlows.PolarPowerFlowJacobian(data, pf.x0)
+
+    df = NLsolve.OnceDifferentiable(pf, J, pf.x0, pf.residual, J.Jv)
+    res = NLsolve.nlsolve(df, pf.x0; nlsolve_solver_kwargs...)
+    if !res.f_converged
+        @error(
+            "The powerflow solver NLSolve did not converge (returned convergence = $(res.f_converged))"
+        )
+    end
+    V = _calc_V(data, res.zero; time_step = time_step)
+    Sbus_result = V .* conj(data.power_network_matrix.data * V)
+    return (res.f_converged, V, Sbus_result)
+end
+
+function _calc_V(
+    data::ACPowerFlowData,
+    x::Vector{Float64};
+    time_step::Int64 = 1,
+)
+    n_buses = length(x) ÷ 2  # Since x has 2 elements per bus (real and imaginary)
+    V = zeros(Complex{Float64}, n_buses)
+    Vm_data = data.bus_magnitude[:, time_step]
+    Va_data = data.bus_angles[:, time_step]
+
+    # Extract values for Vm and Va from x
+    for (ix, b) in enumerate(data.bus_type)
+        if b == PSY.ACBusTypes.REF
+            # For REF bus, we have active and reactive power
+            Vm = Vm_data[ix]
+            Va = Va_data[ix]
+            V[ix] = Vm * exp(im * Va)
+        elseif b == PSY.ACBusTypes.PV
+            # For PV bus, we have reactive power and voltage angle
+            Vm = Vm_data[ix]
+            Va = x[2 * ix]
+            V[ix] = Vm * exp(im * Va)  # Rebuild voltage from magnitude and angle
+        elseif b == PSY.ACBusTypes.PQ
+            # For PQ bus, we have voltage magnitude and voltage angle
+            Vm = x[2 * ix - 1]
+            Va = x[2 * ix]
+            V[ix] = Vm * exp(im * Va)  # Rebuild voltage from magnitude and angle
+        end
+    end
+
+    return V
+end
diff --git a/src/post_processing.jl b/src/post_processing.jl
index ebea6d8..9faa256 100644
--- a/src/post_processing.jl
+++ b/src/post_processing.jl
@@ -617,10 +617,11 @@ dictionary will therefore feature just one key linked to one DataFrame.
         vector containing the reults for one single time-period.
 """
 function write_results(
-    ::ACPowerFlow,
+    ::ACPowerFlow{<:ACPowerFlowSolverType},
     sys::PSY.System,
     data::ACPowerFlowData,
-    result::Vector{Float64},
+    result::Union{Vector{Float64}, Matrix{Float64}};
+    time_step::Int64 = 1,
 )
     @info("Voltages are exported in pu. Powers are exported in MW/MVAr.")
     buses = sort!(collect(PSY.get_components(PSY.Bus, sys)); by = x -> PSY.get_number(x))
@@ -634,27 +635,27 @@ function write_results(
     P_load_vect = fill(0.0, N_BUS)
     Q_load_vect = fill(0.0, N_BUS)
 
-    for (ix, bustype) in enumerate(data.bus_type)
-        P_load_vect[ix] = data.bus_activepower_withdrawals[ix] * sys_basepower
-        Q_load_vect[ix] = data.bus_reactivepower_withdrawals[ix] * sys_basepower
+    for (ix, bustype) in enumerate(data.bus_type[:, time_step])
+        P_load_vect[ix] = data.bus_activepower_withdrawals[ix, time_step] * sys_basepower
+        Q_load_vect[ix] = data.bus_reactivepower_withdrawals[ix, time_step] * sys_basepower
         P_admittance, Q_admittance = _get_fixed_admittance_power(sys, buses[ix], result, ix)
         P_load_vect[ix] += P_admittance * sys_basepower
         Q_load_vect[ix] += Q_admittance * sys_basepower
         if bustype == PSY.ACBusTypes.REF
-            Vm_vect[ix] = data.bus_magnitude[ix]
-            θ_vect[ix] = data.bus_angles[ix]
+            Vm_vect[ix] = data.bus_magnitude[ix, time_step]
+            θ_vect[ix] = data.bus_angles[ix, time_step]
             P_gen_vect[ix] = result[2 * ix - 1] * sys_basepower
             Q_gen_vect[ix] = result[2 * ix] * sys_basepower
         elseif bustype == PSY.ACBusTypes.PV
-            Vm_vect[ix] = data.bus_magnitude[ix]
+            Vm_vect[ix] = data.bus_magnitude[ix, time_step]
             θ_vect[ix] = result[2 * ix]
-            P_gen_vect[ix] = data.bus_activepower_injection[ix] * sys_basepower
+            P_gen_vect[ix] = data.bus_activepower_injection[ix, time_step] * sys_basepower
             Q_gen_vect[ix] = result[2 * ix - 1] * sys_basepower
         elseif bustype == PSY.ACBusTypes.PQ
             Vm_vect[ix] = result[2 * ix - 1]
             θ_vect[ix] = result[2 * ix]
-            P_gen_vect[ix] = data.bus_activepower_injection[ix] * sys_basepower
-            Q_gen_vect[ix] = data.bus_reactivepower_injection[ix] * sys_basepower
+            P_gen_vect[ix] = data.bus_activepower_injection[ix, time_step] * sys_basepower
+            Q_gen_vect[ix] = data.bus_reactivepower_injection[ix, time_step] * sys_basepower
         end
     end
 
@@ -708,25 +709,63 @@ if a new `PowerFlowData` is constructed from the resulting system it is the same
 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)
+function update_system!(sys::PSY.System, data::PowerFlowData; time_step = 1)
     for bus in PSY.get_components(PSY.Bus, sys)
-        if bus.bustype == PSY.ACBusTypes.REF
+        bus_number = data.bus_lookup[PSY.get_number(bus)]
+        bus_type = data.bus_type[bus_number, time_step]  # use this instead of bus.bustype to account for PV -> PQ
+        if bus_type == PSY.ACBusTypes.REF
             # For REF bus, voltage and angle are fixed; update active and reactive
-            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)]]
+            P_gen = data.bus_activepower_injection[bus_number, time_step]
+            Q_gen = data.bus_reactivepower_injection[bus_number, time_step]
             _power_redistribution_ref(sys, P_gen, Q_gen, bus,
                 DEFAULT_MAX_REDISTRIBUTION_ITERATIONS)
-        elseif bus.bustype == PSY.ACBusTypes.PV
+        elseif bus_type == PSY.ACBusTypes.PV
             # For PV bus, active and voltage are fixed; update reactive and angle
-            Q_gen = data.bus_reactivepower_injection[data.bus_lookup[PSY.get_number(bus)]]
+            Q_gen = data.bus_reactivepower_injection[bus_number, time_step]
             _reactive_power_redistribution_pv(sys, Q_gen, bus,
                 DEFAULT_MAX_REDISTRIBUTION_ITERATIONS)
-            PSY.set_angle!(bus, data.bus_angles[data.bus_lookup[PSY.get_number(bus)]])
-        elseif bus.bustype == PSY.ACBusTypes.PQ
+            PSY.set_angle!(bus, data.bus_angles[bus_number, time_step])
+        elseif bus_type == PSY.ACBusTypes.PQ
             # For PQ bus, active and reactive are fixed; update voltage and angle
-            Vm = data.bus_magnitude[data.bus_lookup[PSY.get_number(bus)]]
+            Vm = data.bus_magnitude[bus_number, time_step]
             PSY.set_magnitude!(bus, Vm)
-            PSY.set_angle!(bus, data.bus_angles[data.bus_lookup[PSY.get_number(bus)]])
+            PSY.set_angle!(bus, data.bus_angles[bus_number, time_step])
+            # if it used to be a PV bus, also set the Q value:
+            if bus.bustype == PSY.ACBusTypes.PV
+                Q_gen = data.bus_reactivepower_injection[bus_number, time_step]
+                _reactive_power_redistribution_pv(sys, Q_gen, bus,
+                    DEFAULT_MAX_REDISTRIBUTION_ITERATIONS)
+                # now both the Q and the Vm, Va are correct for this kind of buses
+            end
+        end
+    end
+end
+
+# This cannot work because we do not know the redistribution keys from the power values of devices 
+# to the power values of buses. Leaving this here anyways to document this issue.
+function set_system_time_step(sys::PSY.System, data::PowerFlowData; time_step = 1)
+    @error("This function cannot work")
+    for bus in PSY.get_components(PSY.Bus, sys)
+        bus_number = data.bus_lookup[PSY.get_number(bus)]
+        bus_type = data.bus_type[bus_number, time_step]
+        sys_bus_type = sys.bustype # here we should be using the system bus type to know which inputs to modify
+
+        # for REF bus, we don't need to set anything
+        if sys_bus_type == PSY.ACBusTypes.PV
+            # we need to modify the P value
+            # TODO if we also modify time-series data for voltage set-points, also set that value. 
+            #  However, we cannot do this easily because it can happen that the PV -> PQ transformation 
+            #  caused a different Vm value as the result, while the set-point Vm is still the same. 
+            #  We need a solution for this.
+            P_gen = data.bus_activepower_injection[bus_number, time_step]
+            Q_gen = 0.0
+            # TODO write time-series value of the bus back to devices (not possible)
+        elseif sys_bus_type == PSY.ACBusTypes.PQ
+            P_gen = data.bus_activepower_injection[bus_number, time_step]
+            Q_gen = data.bus_reactivepower_injection[bus_number, time_step]
+            P_load = data.bus_activepower_withdrawals[bus_number, time_step]
+            Q_load = data.bus_reactivepower_withdrawals[bus_number, time_step]
+            # TODO write time-series value of the bus back to devices (not possible)
         end
     end
 end
diff --git a/src/powerflow_types.jl b/src/powerflow_types.jl
index e1b9828..5f69023 100644
--- a/src/powerflow_types.jl
+++ b/src/powerflow_types.jl
@@ -1,9 +1,21 @@
 abstract type PowerFlowEvaluationModel end
+abstract type ACPowerFlowSolverType end
 
-Base.@kwdef struct ACPowerFlow <: PowerFlowEvaluationModel
+struct KLUACPowerFlow <: ACPowerFlowSolverType end
+struct NLSolveACPowerFlow <: ACPowerFlowSolverType end
+
+Base.@kwdef struct ACPowerFlow{ACSolver <: ACPowerFlowSolverType} <:
+                   PowerFlowEvaluationModel
     check_reactive_power_limits::Bool = false
 end
 
+# Create a constructor for ACPowerFlow that defaults to KLUACPowerFlow
+function ACPowerFlow(ACSolver::Type{<:ACPowerFlowSolverType} = KLUACPowerFlow;
+    check_reactive_power_limits::Bool = false,
+)
+    return ACPowerFlow{ACSolver}(check_reactive_power_limits)
+end
+
 struct DCPowerFlow <: PowerFlowEvaluationModel end
 struct PTDFDCPowerFlow <: PowerFlowEvaluationModel end
 struct vPTDFDCPowerFlow <: PowerFlowEvaluationModel end
diff --git a/test/Project.toml b/test/Project.toml
index 7ab27d5..f931fbf 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -13,6 +13,7 @@ PowerFlows = "94fada2c-fd9a-4e89-8d82-81405f5cb4f6"
 PowerNetworkMatrices = "bed98974-b02a-5e2f-9fe0-a103f5c450dd"
 PowerSystemCaseBuilder = "f00506e0-b84f-492a-93c2-c0a9afc4364e"
 PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 3f8827f..590459c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -10,6 +10,8 @@ using CSV
 using DataFrames
 using JSON3
 using DataStructures
+import SparseArrays
+import SparseArrays: SparseMatrixCSC, sparse
 
 const IS = InfrastructureSystems
 const PSB = PowerSystemCaseBuilder
@@ -33,6 +35,9 @@ MAIN_DIR = dirname(@__DIR__)
 
 include("test_utils/common.jl")
 include("test_utils/psse_results_compare.jl")
+#include("test_utils/legacy_pf.jl")
+#Base.eval(PF, Meta.parse("include(\"test/test_utils/legacy_pf.jl\")"))
+Base.eval(PowerFlows, :(include("./test_utils/legacy_pf.jl")))
 
 LOG_FILE = "power-systems.log"
 
diff --git a/test/test_multiperiod_ac_powerflow.jl b/test/test_multiperiod_ac_powerflow.jl
new file mode 100644
index 0000000..dfef418
--- /dev/null
+++ b/test/test_multiperiod_ac_powerflow.jl
@@ -0,0 +1,53 @@
+# Work in progress
+@testset "MULTI-PERIOD power flows evaluation: NR" begin
+    # get system
+    sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
+    injections = CSV.read(
+        joinpath(MAIN_DIR, "test", "test_data", "c_sys14_injections.csv"),
+        DataFrame;
+        header = 0,
+    )
+    withdrawals = CSV.read(
+        joinpath(MAIN_DIR, "test", "test_data", "c_sys14_withdrawals.csv"),
+        DataFrame;
+        header = 0,
+    )
+    flows = CSV.read(
+        joinpath(MAIN_DIR, "test", "test_data", "c_sys14_flows.csv"),
+        DataFrame;
+        header = 0,
+    )
+    angles = CSV.read(
+        joinpath(MAIN_DIR, "test", "test_data", "c_sys14_angles.csv"),
+        DataFrame;
+        header = 0,
+    )
+
+    ##############################################################################
+
+    # create structure for multi-period case
+    pf = ACPowerFlow()
+    time_steps = 24
+    data = PowerFlowData(pf, sys; time_steps = time_steps)
+
+    # allocate data from csv
+    injs = Matrix(injections)
+    withs = Matrix(withdrawals)
+
+    data.bus_activepower_injection .= deepcopy(injs)
+    data.bus_activepower_withdrawals .= deepcopy(withs)
+
+    # get power flows with NR KLU method and write results
+    results = solve_powerflow(pf, data, sys)
+
+    # check results
+    # for t in 1:length(data.timestep_map)
+    #     res_t = solve_powerflow(pf, sys; time_step=t)  # does not work - ts data not set in sys
+    #     flow_ft = res_t["flow_results"].P_from_to
+    #     flow_tf = res_t["flow_results"].P_to_from
+    #     ts_flow_ft = results[data.timestep_map[t]]["flow_results"].P_from_to
+    #     ts_flow_tf = results[data.timestep_map[t]]["flow_results"].P_to_from
+    #     @test isapprox(ts_flow_ft, flow_ft, atol = 1e-9)
+    #     @test isapprox(ts_flow_tf, flow_tf, atol = 1e-9)
+    # end
+end
diff --git a/test/test_nlsolve_powerflow.jl b/test/test_newton_ac_powerflow.jl
similarity index 53%
rename from test/test_nlsolve_powerflow.jl
rename to test/test_newton_ac_powerflow.jl
index b07ae71..beed697 100644
--- a/test/test_nlsolve_powerflow.jl
+++ b/test/test_newton_ac_powerflow.jl
@@ -1,6 +1,10 @@
-@testset "NLsolve Power Flow 14-Bus testing" begin
-    sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
-    set_units_base_system!(sys, UnitSystem.SYSTEM_BASE)
+
+@testset "AC Power Flow 14-Bus testing" for ACSolver in
+                                            (
+    NLSolveACPowerFlow,
+    KLUACPowerFlow,
+    PowerFlows.LUACPowerFlow,
+)
     result_14 = [
         2.3255081760423684
         -0.15529254415401786
@@ -31,79 +35,100 @@
         1.0213119628726421
         -0.2803812119374241
     ]
-    data = PowerFlows.PowerFlowData(ACPowerFlow(), sys; check_connectivity = true)
+
+    sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
+    set_units_base_system!(sys, UnitSystem.SYSTEM_BASE)
+    pf = ACPowerFlow{ACSolver}()
+    data = PowerFlows.PowerFlowData(pf, sys; check_connectivity = true)
     #Compare results between finite diff methods and Jacobian method
-    res1 = PowerFlows._solve_powerflow!(data, false)
-    @test LinearAlgebra.norm(result_14 - res1.zero) <= 1e-6
+    converged1, V1, S1 = PowerFlows._solve_powerflow!(pf, data, false)
+    x1 = PowerFlows._calc_x(data, V1, S1)
+    @test LinearAlgebra.norm(result_14 - x1, Inf) <= 1e-6
 
-    # Test that solve_ac_powerflow! succeeds
+    # Test that solve_powerflow! succeeds
     solved1 = deepcopy(sys)
-    @test solve_ac_powerflow!(solved1)
+    @test solve_powerflow!(pf, solved1)
 
     # Test that passing check_reactive_power_limits=false is the default and violates limits
     solved2 = deepcopy(sys)
-    @test solve_ac_powerflow!(solved2; check_reactive_power_limits = false)
+    @test solve_powerflow!(pf, solved2; check_reactive_power_limits = false)
     @test IS.compare_values(solved1, solved2)
     @test get_reactive_power(get_component(ThermalStandard, solved2, "Bus8")) >
           get_reactive_power_limits(get_component(ThermalStandard, solved2, "Bus8")).max
 
     # Test that passing check_reactive_power_limits=true fixes that
     solved3 = deepcopy(sys)
-    @test solve_ac_powerflow!(solved3; check_reactive_power_limits = true)
+    @test solve_powerflow!(pf, solved3; check_reactive_power_limits = true)
     @test get_reactive_power(get_component(ThermalStandard, solved3, "Bus8")) <=
           get_reactive_power_limits(get_component(ThermalStandard, solved3, "Bus8")).max
 
     # Test Newton method
-    @test solve_ac_powerflow!(deepcopy(sys); method = :newton)
+    @test solve_powerflow!(pf, deepcopy(sys); method = :newton)
 
     # Test enforcing the reactive power limits in closer detail
     set_reactive_power!(get_component(PowerLoad, sys, "Bus4"), 0.0)
-    data = PowerFlows.PowerFlowData(ACPowerFlow(), sys; check_connectivity = true)
-    res2 = PowerFlows._solve_powerflow!(data, true)
-    @test LinearAlgebra.norm(result_14 - res2.zero) >= 1e-6
-    @test 1.08 <= res2.zero[15] <= 1.09
+    data = PowerFlows.PowerFlowData(pf, sys; check_connectivity = true)
+    converged2, V2, S2 = PowerFlows._solve_powerflow!(pf, data, true)
+    x2 = PowerFlows._calc_x(data, V2, S2)
+    @test LinearAlgebra.norm(result_14 - x2, Inf) >= 1e-6
+    @test 1.08 <= x2[15] <= 1.09
 end
 
-@testset "NLsolve Power Flow 14-Bus Line Configurations" begin
+@testset "AC Power Flow 14-Bus Line Configurations" for ACSolver in
+                                                        (
+    NLSolveACPowerFlow,
+    KLUACPowerFlow,
+    PowerFlows.LUACPowerFlow,
+)
     sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
-    base_res = solve_powerflow(ACPowerFlow(), sys)
+    pf = ACPowerFlow{ACSolver}()
+    base_res = solve_powerflow(pf, sys)
     branch = first(PSY.get_components(Line, sys))
     dyn_branch = DynamicBranch(branch)
     add_component!(sys, dyn_branch)
-    @test dyn_pf = solve_ac_powerflow!(sys)
-    dyn_pf = solve_powerflow(ACPowerFlow(), sys)
-    @test LinearAlgebra.norm(dyn_pf["bus_results"].Vm - base_res["bus_results"].Vm) <=
+    @test dyn_pf = solve_powerflow!(pf, sys)
+    dyn_pf = solve_powerflow(pf, sys)
+    @test LinearAlgebra.norm(dyn_pf["bus_results"].Vm - base_res["bus_results"].Vm, Inf) <=
           1e-6
 
     sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
     line = get_component(Line, sys, "Line4")
     PSY.set_available!(line, false)
-    solve_ac_powerflow!(sys)
+    solve_powerflow!(pf, sys)
     @test PSY.get_active_power_flow(line) == 0.0
     test_bus = get_component(PSY.Bus, sys, "Bus 4")
-    @test isapprox(PSY.get_magnitude(test_bus), 1.002; atol = 1e-3)
+    @test isapprox(PSY.get_magnitude(test_bus), 1.002; atol = 1e-3, rtol = 0)
 
     sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
     line = get_component(Line, sys, "Line4")
     PSY.set_available!(line, false)
-    res = solve_powerflow(ACPowerFlow(), sys)
+    res = solve_powerflow(pf, sys)
     @test res["flow_results"].P_from_to[4] == 0.0
     @test res["flow_results"].P_to_from[4] == 0.0
 end
 
-@testset "NLsolve Power Flow 3-Bus Fixed FixedAdmittance testing" begin
+@testset "AC Power Flow 3-Bus Fixed FixedAdmittance testing" for ACSolver in (
+    NLSolveACPowerFlow,
+    KLUACPowerFlow, PowerFlows.LUACPowerFlow,
+)
     p_gen_matpower_3bus = [20.3512373930753, 100.0, 100.0]
     q_gen_matpower_3bus = [45.516916781567232, 10.453799727283879, -31.992561631394636]
     sys_3bus = PSB.build_system(PSB.PSYTestSystems, "psse_3bus_gen_cls_sys")
     bus_103 = get_component(PSY.Bus, sys_3bus, "BUS 3")
     fix_shunt = PSY.FixedAdmittance("FixAdmBus3", true, bus_103, 0.0 + 0.2im)
     add_component!(sys_3bus, fix_shunt)
-    df = solve_powerflow(ACPowerFlow(), sys_3bus)
+    pf = ACPowerFlow{ACSolver}()
+    df = solve_powerflow(pf, sys_3bus)
     @test isapprox(df["bus_results"].P_gen, p_gen_matpower_3bus, atol = 1e-4)
     @test isapprox(df["bus_results"].Q_gen, q_gen_matpower_3bus, atol = 1e-4)
 end
 
-@testset "NLsolve Power Flow convergence fail testing" begin
+@testset "AC Power Flow convergence fail testing" for ACSolver in
+                                                      (
+    NLSolveACPowerFlow,
+    KLUACPowerFlow,
+    PowerFlows.LUACPowerFlow,
+)
     pf_sys5_re = PSB.build_system(PSB.PSITestSystems, "c_sys5_re"; add_forecasts = false)
     remove_component!(Line, pf_sys5_re, "1")
     remove_component!(Line, pf_sys5_re, "2")
@@ -111,15 +136,22 @@ end
     PSY.set_x!(br, 20.0)
     PSY.set_r!(br, 2.0)
 
+    pf = ACPowerFlow{ACSolver}()
+
     # This is a negative test. The data passed for sys5_re is known to be infeasible.
     @test_logs(
         (:error, "The powerflow solver returned convergence = false"),
         match_mode = :any,
-        @test !solve_ac_powerflow!(pf_sys5_re)
+        @test !solve_powerflow!(pf, pf_sys5_re)
     )
 end
 
-@testset "Test 240 Case PSS/e results" begin
+@testset "AC Test 240 Case PSS/e results" for ACSolver in
+                                              (
+    NLSolveACPowerFlow,
+    KLUACPowerFlow,
+    PowerFlows.LUACPowerFlow,
+)
     file = joinpath(
         TEST_FILES_DIR,
         "test_data",
@@ -134,9 +166,11 @@ end
     pf_bus_result_file = joinpath(TEST_FILES_DIR, "test_data", "pf_bus_results.csv")
     pf_gen_result_file = joinpath(TEST_FILES_DIR, "test_data", "pf_gen_results.csv")
 
-    pf = solve_ac_powerflow!(system)
-    @test pf
-    pf_result_df = solve_powerflow(ACPowerFlow(), system)
+    pf = ACPowerFlow{ACSolver}()
+
+    pf1 = solve_powerflow!(pf, system)
+    @test pf1
+    pf_result_df = solve_powerflow(pf, system)
 
     v_diff, angle_diff, number = psse_bus_results_compare(pf_bus_result_file, pf_result_df)
     p_diff, q_diff, names = psse_gen_results_compare(pf_gen_result_file, system)
@@ -152,7 +186,12 @@ end
     @test norm(q_diff, 2) / length(q_diff) < DIFF_L2_TOLERANCE
 end
 
-@testset "Multiple sources at ref" begin
+@testset "AC Multiple sources at ref" for ACSolver in
+                                          (
+    NLSolveACPowerFlow,
+    KLUACPowerFlow,
+    PowerFlows.LUACPowerFlow,
+)
     sys = System(100.0)
     b = ACBus(;
         number = 1,
@@ -186,16 +225,22 @@ end
         X_th = 1e-5,
     )
     add_component!(sys, s2)
-    @test solve_ac_powerflow!(sys)
+    pf = ACPowerFlow{ACSolver}()
+    @test solve_powerflow!(pf, sys)
 
     #Create power mismatch, test for error
     set_active_power!(get_component(Source, sys, "source_1"), -0.4)
     @test_throws ErrorException(
         "Sources do not match P and/or Q requirements for reference bus.",
-    ) solve_ac_powerflow!(sys)
+    ) solve_powerflow!(pf, sys)
 end
 
-@testset "AC PowerFlow with Multiple sources at PV" begin
+@testset "AC PowerFlow with Multiple sources at PV" for ACSolver in
+                                                        (
+    NLSolveACPowerFlow,
+    KLUACPowerFlow,
+    PowerFlows.LUACPowerFlow,
+)
     sys = System(100.0)
     b1 = ACBus(;
         number = 1,
@@ -265,16 +310,24 @@ end
     )
     add_component!(sys, s3)
 
-    @test solve_ac_powerflow!(sys)
+    pf = ACPowerFlow{ACSolver}()
+
+    @test solve_powerflow!(pf, sys)
 
     #Create power mismatch, test for error
     set_reactive_power!(get_component(Source, sys, "source_3"), -0.5)
-    @test_throws ErrorException("Sources do not match Q requirements for PV bus.") solve_ac_powerflow!(
+    @test_throws ErrorException("Sources do not match Q requirements for PV bus.") solve_powerflow!(
+        pf,
         sys,
     )
 end
 
-@testset "AC PowerFlow Source + non-source at Ref" begin
+@testset "AC PowerFlow Source + non-source at Ref" for ACSolver in
+                                                       (
+    NLSolveACPowerFlow,
+    KLUACPowerFlow,
+    PowerFlows.LUACPowerFlow,
+)
     sys = System(100.0)
     b = ACBus(;
         number = 1,
@@ -320,7 +373,9 @@ end
     )
     add_component!(sys, g1)
 
-    @test solve_ac_powerflow!(sys)
+    pf = ACPowerFlow{ACSolver}()
+
+    @test solve_powerflow!(pf, sys)
     @test isapprox(
         get_active_power(get_component(Source, sys, "source_1")),
         0.5;
@@ -333,7 +388,12 @@ end
     )
 end
 
-@testset "AC PowerFlow Source + non-source at PV" begin
+@testset "AC PowerFlow Source + non-source at PV" for ACSolver in
+                                                      (
+    NLSolveACPowerFlow,
+    KLUACPowerFlow,
+    PowerFlows.LUACPowerFlow,
+)
     sys = System(100.0)
     b1 = ACBus(;
         number = 1,
@@ -414,7 +474,9 @@ end
     )
     add_component!(sys, g1)
 
-    @test solve_ac_powerflow!(sys)
+    pf = ACPowerFlow{ACSolver}()
+
+    @test solve_powerflow!(pf, sys)
     @test isapprox(
         get_active_power(get_component(Source, sys, "source_2")),
         0.5;
@@ -426,3 +488,196 @@ end
         atol = 0.001,
     )
 end
+
+# in this test, the following aspects are checked:
+# 1. The results of the power flow are consistent for the KLU and NLSolve solvers
+# 2. The results of the power flow are consistent for the KLU solver and the legacy implementation
+# 3. The Jacobian matrix is the same for the KLU solver and the legacy implementation
+@testset "Compare larger grid results KLU vs NLSolve" begin
+    sys = build_system(MatpowerTestSystems, "matpower_ACTIVSg2000_sys")
+
+    pf_default = ACPowerFlow()
+    pf_klu = ACPowerFlow(KLUACPowerFlow)
+    pf_nlsolve = ACPowerFlow(NLSolveACPowerFlow)
+
+    PSY.set_units_base_system!(sys, "SYSTEM_BASE")
+    data = PowerFlowData(
+        pf_default,
+        sys;
+        check_connectivity = true)
+
+    time_step = 1
+
+    ref = findall(
+        x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF,
+        data.bus_type[:, time_step],
+    )
+    pv = findall(
+        x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV,
+        data.bus_type[:, time_step],
+    )
+    pq = findall(
+        x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ,
+        data.bus_type[:, time_step],
+    )
+    pvpq = [pv; pq]
+
+    npvpq = length(pvpq)
+    npq = length(pq)
+
+    # we need to define lookups for mappings of pv, pq buses onto the internal J indexing
+    pvpq_lookup = zeros(Int64, maximum([ref; pvpq]) + 1)
+    pvpq_lookup[pvpq] .= 1:npvpq
+    pq_lookup = zeros(Int64, maximum([ref; pvpq]) + 1)
+    pq_lookup[pq] .= 1:npq
+
+    # define the internal J indexing using the lookup arrays
+    j_pvpq = pvpq_lookup[pvpq]
+    j_pq = npvpq .+ pq_lookup[pq]
+
+    Vm0 = data.bus_magnitude[:]
+    Va0 = data.bus_angles[:]
+    V0 = Vm0 .* exp.(1im * Va0)
+
+    Ybus = data.power_network_matrix.data
+    rows, cols = SparseArrays.findnz(Ybus)
+
+    diagV = LinearAlgebra.Diagonal(V0)
+    diagIbus_diag = zeros(Complex{Float64}, size(V0, 1))
+    diagIbus = LinearAlgebra.Diagonal(diagIbus_diag)
+
+    diagVnorm = LinearAlgebra.Diagonal(V0 ./ abs.(V0))
+
+    Ybus_diagVnorm = sparse(rows, cols, Complex{Float64}(0))
+    conj_Ybus_diagVnorm = sparse(rows, cols, Complex{Float64}(0))
+    diagV_conj_Ybus_diagVnorm = sparse(rows, cols, Complex{Float64}(0))
+    conj_diagIbus = conj.(diagIbus)
+    conj_diagIbus_diagVnorm = conj.(diagIbus)
+    Ybus_diagV = sparse(rows, cols, Complex{Float64}(0))
+    conj_Ybus_diagV = sparse(rows, cols, Complex{Float64}(0))
+
+    dSbus_dVm = sparse(rows, cols, Complex{Float64}(0))
+    dSbus_dVa = sparse(rows, cols, Complex{Float64}(0))
+    r_dSbus_dVa = sparse(rows, cols, Float64(0))
+    r_dSbus_dVm = sparse(rows, cols, Float64(0))
+    i_dSbus_dVa = sparse(rows, cols, Float64(0))
+    i_dSbus_dVm = sparse(rows, cols, Float64(0))
+
+    J_block = sparse(rows, cols, Float64(0), maximum(rows), maximum(cols), unique)
+    J0_KLU = [J_block[pvpq, pvpq] J_block[pvpq, pq]; J_block[pq, pvpq] J_block[pq, pq]]
+    PF._update_dSbus_dV!(rows, cols, V0, Ybus, diagV, diagVnorm, diagIbus, diagIbus_diag,
+        dSbus_dVa, dSbus_dVm, r_dSbus_dVa, r_dSbus_dVm, i_dSbus_dVa, i_dSbus_dVm,
+        Ybus_diagVnorm, conj_Ybus_diagVnorm, diagV_conj_Ybus_diagVnorm, conj_diagIbus,
+        conj_diagIbus_diagVnorm, Ybus_diagV, conj_Ybus_diagV)
+    PF._update_J!(
+        J0_KLU,
+        r_dSbus_dVa,
+        r_dSbus_dVm,
+        i_dSbus_dVa,
+        i_dSbus_dVm,
+        pvpq,
+        pq,
+        j_pvpq,
+        j_pq,
+    )
+
+    dSbus_dVa0_LU, dSbus_dVm0_LU = PF._legacy_dSbus_dV(V0, Ybus)
+    J0_LU = PF._legacy_J(dSbus_dVa, dSbus_dVm, pvpq, pq)
+
+    @test all(isapprox.(J0_LU, J0_KLU, rtol = 0, atol = 1e-12))
+    @test all(isapprox.(J0_LU.nzval, J0_KLU.nzval, rtol = 0, atol = 1e-12))
+    @test J0_KLU.rowval == J0_LU.rowval
+    @test J0_KLU.colptr == J0_LU.colptr
+
+    res_default = solve_powerflow(pf_default, sys)  # must be the same as KLU
+    res_klu = solve_powerflow(pf_klu, sys)
+    res_nlsolve = solve_powerflow(pf_nlsolve, sys)
+
+    @test all(
+        isapprox.(
+            res_klu["bus_results"][!, :Vm],
+            res_default["bus_results"][!, :Vm],
+            rtol = 0,
+            atol = 1e-12,
+        ),
+    )
+    @test all(
+        isapprox.(
+            res_klu["bus_results"][!, :θ],
+            res_default["bus_results"][!, :θ],
+            rtol = 0,
+            atol = 1e-12,
+        ),
+    )
+
+    @test all(
+        isapprox.(
+            res_klu["bus_results"][!, :Vm],
+            res_nlsolve["bus_results"][!, :Vm],
+            rtol = 0,
+            atol = 1e-8,
+        ),
+    )
+    @test all(
+        isapprox.(
+            res_klu["bus_results"][!, :θ],
+            res_nlsolve["bus_results"][!, :θ],
+            rtol = 0,
+            atol = 1e-8,
+        ),
+    )
+
+    # test against legacy implementation
+    pf_legacy = ACPowerFlow(PowerFlows.LUACPowerFlow)
+    res_legacy = solve_powerflow(pf_legacy, sys)
+
+    @test all(
+        isapprox.(
+            res_klu["bus_results"][!, :Vm],
+            res_legacy["bus_results"][!, :Vm],
+            rtol = 0,
+            atol = 1e-12,
+        ),
+    )
+    @test all(
+        isapprox.(
+            res_klu["bus_results"][!, :θ],
+            res_legacy["bus_results"][!, :θ],
+            rtol = 0,
+            atol = 1e-12,
+        ),
+    )
+
+    Vm1_KLU = res_klu["bus_results"][!, :Vm]
+    Va1_KLU = res_klu["bus_results"][!, :θ]
+    V1_KLU = Vm1_KLU .* exp.(1im * Va1_KLU)
+
+    Vm1_LU = res_legacy["bus_results"][!, :Vm]
+    Va1_LU = res_legacy["bus_results"][!, :θ]
+    V1_LU = Vm1_LU .* exp.(1im * Va1_LU)
+
+    J1_KLU = [J_block[pvpq, pvpq] J_block[pvpq, pq]; J_block[pq, pvpq] J_block[pq, pq]]
+    PF._update_dSbus_dV!(rows, cols, V1_KLU, Ybus, diagV, diagVnorm, diagIbus,
+        diagIbus_diag, dSbus_dVa, dSbus_dVm, r_dSbus_dVa, r_dSbus_dVm, i_dSbus_dVa,
+        i_dSbus_dVm, Ybus_diagVnorm, conj_Ybus_diagVnorm, diagV_conj_Ybus_diagVnorm,
+        conj_diagIbus, conj_diagIbus_diagVnorm, Ybus_diagV, conj_Ybus_diagV)
+    PF._update_J!(
+        J1_KLU,
+        r_dSbus_dVa,
+        r_dSbus_dVm,
+        i_dSbus_dVa,
+        i_dSbus_dVm,
+        pvpq,
+        pq,
+        j_pvpq,
+        j_pq,
+    )
+
+    dSbus_dVa1_LU, dSbus_dVm1_LU = PF._legacy_dSbus_dV(V1_LU, Ybus)
+    J1_LU = PF._legacy_J(dSbus_dVa1_LU, dSbus_dVm1_LU, pvpq, pq)
+
+    @test all(isapprox.(J1_LU, J1_KLU, rtol = 0, atol = 1e-10))
+    @test all(isapprox.(J1_LU.nzval, J1_KLU.nzval, rtol = 0, atol = 1e-10))
+    @test J1_KLU.rowval == J1_LU.rowval
+    @test J1_KLU.colptr == J1_LU.colptr
+end
diff --git a/test/test_powerflow_data.jl b/test/test_powerflow_data.jl
index 75964cb..2d45a14 100644
--- a/test/test_powerflow_data.jl
+++ b/test/test_powerflow_data.jl
@@ -1,6 +1,7 @@
 @testset "PowerFlowData" begin
     sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
-    @test PowerFlowData(ACPowerFlow(), sys) isa PF.ACPowerFlowData
+    @test PowerFlowData(ACPowerFlow{NLSolveACPowerFlow}(), sys) isa PF.ACPowerFlowData
+    @test PowerFlowData(ACPowerFlow{KLUACPowerFlow}(), 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
@@ -17,23 +18,24 @@ end
           PF.vPTDFPowerFlowData
 end
 
-@testset "System <-> PowerFlowData round trip" begin
+@testset "System <-> PowerFlowData round trip" for ACSolver in
+                                                   (NLSolveACPowerFlow, KLUACPowerFlow)
     # 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)
+    data_original = PowerFlowData(ACPowerFlow{ACSolver}(), sys_original)
 
     sys_modified = deepcopy(sys_original)
     modify_rts_system!(sys_modified)
-    data_modified = PowerFlowData(ACPowerFlow(), sys_original)
+    data_modified = PowerFlowData(ACPowerFlow{ACSolver}(), 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)
+    data_null_updated = PowerFlowData(ACPowerFlow{ACSolver}(), sys_null_updated)
     @test IS.compare_values(powerflow_match_fn, data_null_updated, data_original)
 
     # Modified versions should not be the same as unmodified versions
@@ -47,7 +49,7 @@ end
     # Constructing PowerFlowData from modified system should result in data_modified
     @test IS.compare_values(
         powerflow_match_fn,
-        PowerFlowData(ACPowerFlow(), sys_modified),
+        PowerFlowData(ACPowerFlow{ACSolver}(), sys_modified),
         data_modified,
     )
 
@@ -56,6 +58,9 @@ end
     sys_modify_updated = deepcopy(sys_original)
     PF.update_system!(sys_modify_updated, data_modified)
     sys_mod_redist = deepcopy(sys_modified)
-    PF.update_system!(sys_mod_redist, PowerFlowData(ACPowerFlow(), sys_mod_redist))
+    PF.update_system!(
+        sys_mod_redist,
+        PowerFlowData(ACPowerFlow{ACSolver}(), 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 d6abd32..802eee3 100644
--- a/test/test_psse_export.jl
+++ b/test/test_psse_export.jl
@@ -222,9 +222,14 @@ function compare_systems_loosely(sys1::PSY.System, sys2::PSY.System;
     return result
 end
 
-function test_power_flow(sys1::System, sys2::System; exclude_reactive_flow = false)
-    result1 = solve_powerflow(ACPowerFlow(), sys1)
-    result2 = solve_powerflow(ACPowerFlow(), sys2)
+function test_power_flow(
+    pf::ACPowerFlow{<:ACPowerFlowSolverType},
+    sys1::System,
+    sys2::System;
+    exclude_reactive_flow = false,
+)
+    result1 = solve_powerflow(pf, sys1)
+    result2 = solve_powerflow(pf, sys2)
     reactive_power_tol =
         exclude_reactive_flow ? nothing : POWERFLOW_COMPARISON_TOLERANCE
     @test compare_df_within_tolerance("bus_results", result1["bus_results"],
@@ -246,6 +251,7 @@ read_system_and_metadata(export_subdir) = read_system_and_metadata(
     get_psse_export_paths(export_subdir)...)
 
 function test_psse_round_trip(
+    pf::ACPowerFlow{<:ACPowerFlowSolverType},
     sys::System,
     exporter::PSSEExporter,
     scenario_name::AbstractString,
@@ -265,7 +271,7 @@ function test_psse_round_trip(
     sys2, sys2_metadata = read_system_and_metadata(raw_path, metadata_path)
     @test compare_systems_loosely(sys, sys2)
     do_power_flow_test &&
-        test_power_flow(sys, sys2; exclude_reactive_flow = exclude_reactive_flow)
+        test_power_flow(pf, sys, sys2; exclude_reactive_flow = exclude_reactive_flow)
 end
 
 "Test that the two raw files are exactly identical and the two metadata files parse to identical JSON"
@@ -318,17 +324,21 @@ end
     @test compare_systems_loosely(sys, deepcopy(sys))
 end
 
-@testset "PSSE Exporter with system_240[32].json, v33" begin
+@testset "PSSE Exporter with system_240[32].json, v33" for (ACSolver, folder_name) in (
+    (NLSolveACPowerFlow, "system_240_NLSolve"),
+    (KLUACPowerFlow, "system_240_KLU"),
+)
     sys = load_test_system()
+    pf = ACPowerFlow{ACSolver}()
     isnothing(sys) && return
 
     # PSS/E version must be one of the supported ones
     @test_throws ArgumentError PSSEExporter(sys, :vNonexistent, test_psse_export_dir)
 
     # Reimported export should be comparable to original system
-    export_location = joinpath(test_psse_export_dir, "v33", "system_240")
+    export_location = joinpath(test_psse_export_dir, "v33", folder_name)
     exporter = PSSEExporter(sys, :v33, export_location)
-    test_psse_round_trip(sys, exporter, "basic", export_location;
+    test_psse_round_trip(pf, 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
@@ -360,20 +370,24 @@ end
     @test_logs((:error, r"values do not match"),
         match_mode = :any, min_level = Logging.Error,
         compare_systems_loosely(sys, reread_sys2))
-    test_power_flow(sys2, reread_sys2; exclude_reactive_flow = true)  # TODO why is reactive flow not matching?
+    test_power_flow(pf, sys2, reread_sys2; exclude_reactive_flow = true)  # TODO why is reactive flow not matching?
 end
 
-@testset "PSSE Exporter with RTS_GMLC_DA_sys, v33" begin
+@testset "PSSE Exporter with RTS_GMLC_DA_sys, v33" for (ACSolver, folder_name) in (
+    (NLSolveACPowerFlow, "rts_gmlc_NLSolve"),
+    (KLUACPowerFlow, "rts_gmlc_KLU"),
+)
     sys = create_pf_friendly_rts_gmlc()
+    pf = ACPowerFlow{ACSolver}()
     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_psse_export_dir)
 
     # Reimported export should be comparable to original system
-    export_location = joinpath(test_psse_export_dir, "v33", "rts_gmlc")
+    export_location = joinpath(test_psse_export_dir, "v33", folder_name)
     exporter = PSSEExporter(sys, :v33, export_location)
-    test_psse_round_trip(sys, exporter, "basic", export_location;
+    test_psse_round_trip(pf, 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
@@ -404,11 +418,11 @@ end
     @test_logs((:error, r"values do not match"),
         match_mode = :any, min_level = Logging.Error,
         compare_systems_loosely(sys, reread_sys2))
-    test_power_flow(sys2, reread_sys2; exclude_reactive_flow = true)  # TODO why is reactive flow not matching?
+    test_power_flow(pf, 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, export_location)
-    pf2 = PowerFlowData(ACPowerFlow(), sys)
+    pf2 = PowerFlowData(pf, 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)
@@ -421,11 +435,11 @@ end
     @test_logs((:error, r"values do not match"),
         match_mode = :any, min_level = Logging.Error,
         compare_systems_loosely(sys, reread_sys3))
-    test_power_flow(sys2, reread_sys3; exclude_reactive_flow = true)  # TODO why is reactive flow not matching?
+    test_power_flow(pf, 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, export_location; write_comments = true)
-    test_psse_round_trip(sys, exporter, "basic6", export_location;
+    test_psse_round_trip(pf, sys, exporter, "basic6", export_location;
         exclude_reactive_flow = true)  # TODO why is reactive flow not matching?
 end
 
diff --git a/test/test_utils/legacy_pf.jl b/test/test_utils/legacy_pf.jl
new file mode 100644
index 0000000..5436b33
--- /dev/null
+++ b/test/test_utils/legacy_pf.jl
@@ -0,0 +1,120 @@
+
+import PowerFlows
+
+struct LUACPowerFlow <: ACPowerFlowSolverType end  # Only for testing, a basic implementation using LinearAlgebra.lu, allocates a lot of memory
+
+# this function is for testing purposes only
+function _legacy_dSbus_dV(
+    V::Vector{Complex{Float64}},
+    Ybus::SparseMatrixCSC{Complex{Float64}, Int64},
+)
+    diagV = LinearAlgebra.Diagonal(V)
+    diagVnorm = LinearAlgebra.Diagonal(V ./ abs.(V))
+    diagIbus = LinearAlgebra.Diagonal(Ybus * V)
+    dSbus_dVm = diagV * conj.(Ybus * diagVnorm) + conj.(diagIbus) * diagVnorm
+    dSbus_dVa = 1im * diagV * conj.(diagIbus - Ybus * diagV)
+    return dSbus_dVa, dSbus_dVm
+end
+
+# this function is for testing purposes only
+function _legacy_J(
+    dSbus_dVa::SparseMatrixCSC{Complex{Float64}, Int64},
+    dSbus_dVm::SparseMatrixCSC{Complex{Float64}, Int64},
+    pvpq::Vector{Int64},
+    pq::Vector{Int64},
+)
+    j11 = real(dSbus_dVa[pvpq, pvpq])
+    j12 = real(dSbus_dVm[pvpq, pq])
+    j21 = imag(dSbus_dVa[pq, pvpq])
+    j22 = imag(dSbus_dVm[pq, pq])
+    J = sparse([j11 j12; j21 j22])
+    return J
+end
+
+# legacy NR implementation - here we do not care about allocations, we use this function only for testing purposes
+function _newton_powerflow(
+    pf::ACPowerFlow{PowerFlows.LUACPowerFlow},
+    data::PowerFlows.ACPowerFlowData;
+    time_step::Int64 = 1,
+    kwargs...,
+)
+    # Fetch maxIter and tol from kwargs, or use defaults if not provided
+    maxIter = get(kwargs, :maxIter, DEFAULT_NR_MAX_ITER)
+    tol = get(kwargs, :tol, DEFAULT_NR_TOL)
+    i = 0
+
+    Ybus = data.power_network_matrix.data
+
+    # Find indices for each bus type
+    #ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type)
+    pv = findall(
+        x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV,
+        data.bus_type[:, time_step],
+    )
+    pq = findall(
+        x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ,
+        data.bus_type[:, time_step],
+    )
+    pvpq = [pv; pq]
+
+    #nref = length(ref)
+    npv = length(pv)
+    npq = length(pq)
+    npvpq = npv + npq
+    n_buses = length(data.bus_type)
+
+    Vm = data.bus_magnitude[:, time_step]
+    # prevent unfeasible starting values for Vm; for pv and ref buses we cannot do this:
+    Vm[pq] .= clamp.(Vm[pq], 0.9, 1.1)
+    Va = data.bus_angles[:, time_step]
+    V = zeros(Complex{Float64}, length(Vm))
+    V .= Vm .* exp.(1im .* Va)
+
+    # pre-allocate dx
+    dx = zeros(Float64, npv + 2 * npq)
+
+    Ybus = data.power_network_matrix.data
+
+    Sbus =
+        data.bus_activepower_injection[:, time_step] -
+        data.bus_activepower_withdrawals[:, time_step] +
+        1im * (
+            data.bus_reactivepower_injection[:, time_step] -
+            data.bus_reactivepower_withdrawals[:, time_step]
+        )
+
+    mis = V .* conj.(Ybus * V) .- Sbus
+    F = [real(mis[pvpq]); imag(mis[pq])]
+
+    converged = npvpq == 0
+
+    while i < maxIter && !converged
+        i += 1
+        dSbus_dVa, dSbus_dVm = _legacy_dSbus_dV(V, Ybus)
+        J = _legacy_J(dSbus_dVa, dSbus_dVm, pvpq, pq)
+
+        # using a different factorization that KLU for testing
+        factor_J = LinearAlgebra.lu(J)
+        dx .= factor_J \ F
+
+        Va[pv] .-= dx[1:npv]
+        Va[pq] .-= dx[(npv + 1):(npv + npq)]
+        Vm[pq] .-= dx[(npv + npq + 1):(npv + 2 * npq)]
+        V .= Vm .* exp.(1im .* Va)
+        Vm .= abs.(V)
+        Va .= angle.(V)
+
+        mis = V .* conj.(Ybus * V) .- Sbus
+        F .= [real(mis[pvpq]); imag(mis[pq])]
+
+        converged = LinearAlgebra.norm(F, Inf) < tol
+    end
+
+    if !converged
+        @error("The powerflow solver with KLU did not converge after $i iterations")
+    else
+        @info("The powerflow solver with KLU converged after $i iterations")
+    end
+    Sbus_result = V .* conj(Ybus * V)
+    return (converged, V, Sbus_result)
+end
diff --git a/test/test_utils/psse_results_compare.jl b/test/test_utils/psse_results_compare.jl
index 361cda8..dfe162b 100644
--- a/test/test_utils/psse_results_compare.jl
+++ b/test/test_utils/psse_results_compare.jl
@@ -1,4 +1,4 @@
-function psse_bus_results_compare(file_name, pf_results)
+function psse_bus_results_compare(file_name::String, pf_results::Dict)
     pf_result_bus = CSV.read(file_name, DataFrame)
 
     v_diff = Float64[]
@@ -15,6 +15,10 @@ function psse_bus_results_compare(file_name, pf_results)
     return v_diff, angle_diff, number
 end
 
+function psse_bus_results_compare(file_name::String, pf_results::Missing)
+    throw(ArgumentError("pf_results are missing - calculation failed"))
+end
+
 function psse_gen_results_compare(file_name, system::PSY.System)
     base_power = get_base_power(system)
     pf_result_gen = CSV.read(file_name, DataFrame)