diff --git a/Project.toml b/Project.toml
index 63c022a2..c8cd422b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -10,7 +10,6 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 PowerNetworkMatrices = "bed98974-b02a-5e2f-9fe0-a103f5c450dd"
 PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"
-SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
diff --git a/docs/src/modeler_guide/power_flow.md b/docs/src/modeler_guide/power_flow.md
index 77679e6b..eea10016 100644
--- a/docs/src/modeler_guide/power_flow.md
+++ b/docs/src/modeler_guide/power_flow.md
@@ -67,7 +67,7 @@ For instance, initializing dynamic simulations. Also, because [`run_powerflow!`]
 are also available for [`run_powerflow!`](@ref)
 
 ````@example generated_power_flow
-run_powerflow!(system_data; finite_diff = true, method = :newton)
+run_powerflow!(system_data; method = :newton)
 ````
 
 After running the power flow command this are the values of the
diff --git a/src/PowerFlowData.jl b/src/PowerFlowData.jl
index 325ad17b..669123ac 100644
--- a/src/PowerFlowData.jl
+++ b/src/PowerFlowData.jl
@@ -1,72 +1,92 @@
 """
-Structure containing all the data required for the evaluation of the power 
+Structure containing all the data required for the evaluation of the power
 flows and angles, as well as these ones.
 
 # Arguments:
 - `bus_lookup::Dict{Int, Int}`:
-        dictionary linking the system's bus number with the rows of either 
+        dictionary linking the system's bus number with the rows of either
         "power_network_matrix" or "aux_network_matrix".
 - `branch_lookup::Dict{String, Int}`:
-        dictionary linking the branch name with the column name of either the 
+        dictionary linking the branch name with the column name of either the
         "power_network_matrix" or "aux_network_matrix".
 - `bus_activepower_injection::Matrix{Float64}`:
-        "(b, t)" matrix containing the bus active power injection. b: number of 
+        "(b, t)" matrix containing the bus active power injection. b: number of
         buses, t: number of time period.
 - `bus_reactivepower_injection::Matrix{Float64}`:
-        "(b, t)" matrix containing the bus reactive power injection. b: number 
+        "(b, t)" matrix containing the bus reactive power injection. b: number
         of buses, t: number of time period.
 - `bus_activepower_withdrawals::Matrix{Float64}`:
-        "(b, t)" matrix containing the bus reactive power withdrawals. b: 
+        "(b, t)" matrix containing the bus reactive power withdrawals. b:
         number of buses, t: number of time period.
 - `bus_reactivepower_withdrawals::Matrix{Float64}`:
-        "(b, t)" matrix containing the bus reactive power withdrawals. b: 
+        "(b, t)" matrix containing the bus reactive power withdrawals. b:
         number of buses, t: number of time period.
+- `bus_reactivepower_bounds::Vector{Float64}`: Upper and Lower bounds for the reactive supply
+        at each bus.
 - `bus_type::Vector{PSY.ACBusTypes}`:
-        vector containing type of buses present in the system, ordered 
+        vector containing type of buses present in the system, ordered
         according to "bus_lookup".
 - `bus_magnitude::Matrix{Float64}`:
-        "(b, t)" matrix containing the bus magnitudes, ordered according to 
+        "(b, t)" matrix containing the bus magnitudes, ordered according to
         "bus_lookup". b: number of buses, t: number of time period.
 - `bus_angles::Matrix{Float64}`:
-        "(b, t)" matrix containing the bus angles, ordered according to 
+        "(b, t)" matrix containing the bus angles, ordered according to
         "bus_lookup". b: number of buses, t: number of time period.
 - `branch_flow_values::Matrix{Float64}`:
-        "(br, t)" matrix containing the power flows, ordered according to 
+        "(br, t)" matrix containing the power flows, ordered according to
         "branch_lookup". br: number of branches, t: number of time period.
 - `timestep_map::Dict{Int, S}`:
-        dictonary mapping the number of the time periods (corresponding to the 
-        column number of the previosly mentioned matrices) and their actual 
-        names.
+        dictonary mapping the number of the time periods (corresponding to the
+        column number of the previosly mentioned matrices) and their names.
 - `valid_ix::Vector{Int}`:
-        vector containing the indeces related to those buses that are not slack
-        ones.
+        vector containing the indeces of not slack buses
 - `power_network_matrix::M`:
-        matrix used for the evaluation of either the power flows or bus angles, 
+        matrix used for the evaluation of either the power flows or bus angles,
         depending on the method considered.
 - `aux_network_matrix::N`:
-        matrix used for the evaluation of either the power flows or bus angles, 
+        matrix used for the evaluation of either the power flows or bus angles,
         depending on the method considered.
-
+- `neighbors::Vector{Set{Int}}`: Vector with the sets of adjacent buses.
 """
-struct PowerFlowData{M <: PNM.PowerNetworkMatrix, N, S <: Union{String, Char}}
+struct PowerFlowData{
+    M <: PNM.PowerNetworkMatrix,
+    N <: Union{PNM.PowerNetworkMatrix, Nothing},
+}
     bus_lookup::Dict{Int, Int}
     branch_lookup::Dict{String, Int}
     bus_activepower_injection::Matrix{Float64}
     bus_reactivepower_injection::Matrix{Float64}
     bus_activepower_withdrawals::Matrix{Float64}
     bus_reactivepower_withdrawals::Matrix{Float64}
+    bus_reactivepower_bounds::Vector{Vector{Float64}}
     bus_type::Vector{PSY.ACBusTypes}
     bus_magnitude::Matrix{Float64}
     bus_angles::Matrix{Float64}
     branch_flow_values::Matrix{Float64}
-    timestep_map::Dict{Int, S}
+    timestep_map::Dict{Int, String}
     valid_ix::Vector{Int}
     power_network_matrix::M
     aux_network_matrix::N
+    neighbors::Vector{Set{Int}}
 end
 
 # AC Power Flow Data
 # TODO -> MULTI PERIOD: AC Power Flow Data
+function _calculate_neighbors(
+    Yb::PNM.Ybus{
+        Tuple{Vector{Int64}, Vector{Int64}},
+        Tuple{Dict{Int64, Int64}, Dict{Int64, Int64}},
+    },
+)
+    I, J, V = SparseArrays.findnz(Yb.data)
+    neighbors = [Set{Int}([i]) for i in 1:length(Yb.axes[1])]
+    for nz in eachindex(V)
+        push!(neighbors[I[nz]], J[nz])
+        push!(neighbors[J[nz]], I[nz])
+    end
+    return neighbors
+end
+
 """
 Function for the definition of the PowerFlowData strucure given the System
 data, number of time periods to consider and their names.
@@ -80,20 +100,23 @@ NOTE: use it for AC power flow computations.
         container storing the system data to consider in the PowerFlowData
         structure.
 - `timesteps::Int`:
-        number of time periods to consider in the PowerFlowData structure. It 
+        number of time periods to consider in the PowerFlowData structure. It
         defines the number of columns of the matrices used to store data.
         Default value = 1.
 - `timestep_names::Vector{String}`:
         names of the time periods defines by the argmunet "timesteps". Default
         value = String[].
+- `check_connectivity::Bool`:
+        Perform connectivity check on the network matrix. Default value = true.
 
 WARNING: functions for the evaluation of the multi-period AC OPF still to be implemented.
 """
 function PowerFlowData(
     ::ACPowerFlow,
-    sys::PSY.System,
+    sys::PSY.System;
     timesteps::Int = 1,
-    timestep_names::Vector{String} = String[])
+    timestep_names::Vector{String} = String[],
+    check_connectivity::Bool = true)
 
     # assign timestep_names
     # timestep names are then allocated in a dictionary to map matrix columns
@@ -106,7 +129,7 @@ function PowerFlowData(
     end
 
     # get data for calculations
-    power_network_matrix = PNM.Ybus(sys)
+    power_network_matrix = PNM.Ybus(sys; check_connectivity = check_connectivity)
 
     # get number of buses and branches
     n_buses = length(axes(power_network_matrix, 1))
@@ -119,6 +142,8 @@ function PowerFlowData(
     bus_lookup = power_network_matrix.lookup[2]
     branch_lookup =
         Dict{String, Int}(PSY.get_name(b) => ix for (ix, b) in enumerate(branches))
+
+    # TODO: bus_type might need to also be a Matrix since the type can change for a particular scenario
     bus_type = Vector{PSY.ACBusTypes}(undef, n_buses)
     bus_angles = zeros(Float64, n_buses)
     bus_magnitude = zeros(Float64, n_buses)
@@ -126,7 +151,7 @@ function PowerFlowData(
         PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.Bus, sys)
     )
 
-    for (ix, bus_no) in bus_lookup
+    for (bus_no, ix) in bus_lookup
         bus_name = temp_bus_map[bus_no]
         bus = PSY.get_component(PSY.Bus, sys, bus_name)
         bus_type[ix] = PSY.get_bustype(bus)
@@ -156,18 +181,20 @@ function PowerFlowData(
         sys,
     )
 
-    # initialize data
-    init_1 = zeros(n_buses, timesteps)
-    init_2 = zeros(n_branches, timesteps)
-
     # define fields as matrices whose number of columns is eqault to the number of timesteps
-    bus_activepower_injection_1 = deepcopy(init_1)
-    bus_reactivepower_injection_1 = deepcopy(init_1)
-    bus_activepower_withdrawals_1 = deepcopy(init_1)
-    bus_reactivepower_withdrawals_1 = deepcopy(init_1)
-    bus_magnitude_1 = deepcopy(init_1)
-    bus_angles_1 = deepcopy(init_1)
-    branch_flow_values_1 = deepcopy(init_2)
+    bus_activepower_injection_1 = zeros(n_buses, timesteps)
+    bus_reactivepower_injection_1 = zeros(n_buses, timesteps)
+    bus_activepower_withdrawals_1 = zeros(n_buses, timesteps)
+    bus_reactivepower_withdrawals_1 = zeros(n_buses, timesteps)
+    bus_magnitude_1 = zeros(n_buses, timesteps)
+    bus_angles_1 = zeros(n_buses, timesteps)
+    branch_flow_values_1 = zeros(n_branches, timesteps)
+
+    bus_reactivepower_bounds = Vector{Vector{Float64}}(undef, n_buses)
+    for i in 1:n_buses
+        bus_reactivepower_bounds[i] = [0.0, 0.0]
+    end
+    _get_reactive_power_bound!(bus_reactivepower_bounds, bus_lookup, sys)
 
     # initial values related to first timestep allocated in the first column
     bus_activepower_injection_1[:, 1] .= bus_activepower_injection
@@ -185,14 +212,16 @@ function PowerFlowData(
         bus_reactivepower_injection_1,
         bus_activepower_withdrawals_1,
         bus_reactivepower_withdrawals_1,
+        bus_reactivepower_bounds,
         bus_type,
         bus_magnitude_1,
         bus_angles_1,
         branch_flow_values_1,
-        Dict(zip([1], "1")),
+        Dict(1 => "1"),
         setdiff(1:n_buses, ref_bus_positions),
         power_network_matrix,
         nothing,
+        _calculate_neighbors(power_network_matrix),
     )
 end
 
@@ -205,24 +234,27 @@ NOTE: use it for DC power flow computations.
 
 # Arguments:
 - `::DCPowerFlow`:
-        use DCPowerFlow() to store the ABA matrix as power_network_matrix and 
+        use DCPowerFlow() to store the ABA matrix as power_network_matrix and
         the BA matrix as aux_network_matrix.
 - `sys::PSY.System`:
         container storing the system data to consider in the PowerFlowData
         structure.
 - `timesteps::Int`:
-        number of time periods to consider in the PowerFlowData structure. It 
+        number of time periods to consider in the PowerFlowData structure. It
         defines the number of columns of the matrices used to store data.
         Default value = 1.
 - `timestep_names::Vector{String}`:
         names of the time periods defines by the argmunet "timesteps". Default
         value = String[].
+- `check_connectivity::Bool`:
+        Perform connectivity check on the network matrix. Default value = true.
 """
 function PowerFlowData(
     ::DCPowerFlow,
-    sys::PSY.System,
+    sys::PSY.System;
     timesteps::Int = 1,
-    timestep_names::Vector{String} = String[])
+    timestep_names::Vector{String} = String[],
+    check_connectivity::Bool = true)
 
     # assign timestep_names
     # timestep names are then allocated in a dictionary to map matrix columns
@@ -293,8 +325,6 @@ function PowerFlowData(
     bus_angles_1 = deepcopy(init_1)
     branch_flow_values_1 = deepcopy(init_2)
 
-    # @show size(bus_activepower_injection_1)
-
     # initial values related to first timestep allocated in the first column
     bus_activepower_injection_1[:, 1] .= bus_activepower_injection
     bus_reactivepower_injection_1[:, 1] .= bus_reactivepower_injection
@@ -311,6 +341,7 @@ function PowerFlowData(
         bus_reactivepower_injection_1,
         bus_activepower_withdrawals_1,
         bus_reactivepower_withdrawals_1,
+        Vector{Vector{Float64}}(),
         bus_type,
         bus_magnitude_1,
         bus_angles_1,
@@ -319,6 +350,7 @@ function PowerFlowData(
         setdiff(1:n_buses, aux_network_matrix.ref_bus_positions),
         power_network_matrix,
         aux_network_matrix,
+        Vector{Set{Int}}(),
     )
 end
 
@@ -331,13 +363,13 @@ NOTE: use it for DC power flow computations.
 
 # Arguments:
 - `::PTDFDCPowerFlow`:
-        use PTDFDCPowerFlow() to store the PTDF matrix as power_network_matrix 
+        use PTDFDCPowerFlow() to store the PTDF matrix as power_network_matrix
         and the ABA matrix as aux_network_matrix.
 - `sys::PSY.System`:
         container storing the system data to consider in the PowerFlowData
         structure.
 - `timesteps::Int`:
-        number of time periods to consider in the PowerFlowData structure. It 
+        number of time periods to consider in the PowerFlowData structure. It
         defines the number of columns of the matrices used to store data.
         Default value = 1.
 - `timestep_names::Vector{String}`:
@@ -346,9 +378,10 @@ NOTE: use it for DC power flow computations.
 """
 function PowerFlowData(
     ::PTDFDCPowerFlow,
-    sys::PSY.System,
+    sys::PSY.System;
     timesteps::Int = 1,
-    timestep_names::Vector{String} = String[])
+    timestep_names::Vector{String} = String[],
+    check_connectivity::Bool = true)
 
     # assign timestep_names
     # timestep names are then allocated in a dictionary to map matrix columns
@@ -437,6 +470,7 @@ function PowerFlowData(
         bus_reactivepower_injection_1,
         bus_activepower_withdrawals_1,
         bus_reactivepower_withdrawals_1,
+        Vector{Vector{Float64}}(),
         bus_type,
         bus_magnitude_1,
         bus_angles_1,
@@ -445,6 +479,7 @@ function PowerFlowData(
         setdiff(1:n_buses, aux_network_matrix.ref_bus_positions),
         power_network_matrix,
         aux_network_matrix,
+        Vector{Set{Int}}(),
     )
 end
 
@@ -457,13 +492,13 @@ NOTE: use it for DC power flow computations.
 
 # Arguments:
 - `::PTDFDCPowerFlow`:
-        use vPTDFDCPowerFlow() to store the Virtual PTDF matrix as 
+        use vPTDFDCPowerFlow() to store the Virtual PTDF matrix as
         power_network_matrix and the ABA matrix as aux_network_matrix.
 - `sys::PSY.System`:
         container storing the system data to consider in the PowerFlowData
         structure.
 - `timesteps::Int`:
-        number of time periods to consider in the PowerFlowData structure. It 
+        number of time periods to consider in the PowerFlowData structure. It
         defines the number of columns of the matrices used to store data.
         Default value = 1.
 - `timestep_names::Vector{String}`:
@@ -472,9 +507,10 @@ NOTE: use it for DC power flow computations.
 """
 function PowerFlowData(
     ::vPTDFDCPowerFlow,
-    sys::PSY.System,
+    sys::PSY.System;
     timesteps::Int = 1,
-    timestep_names::Vector{String} = String[])
+    timestep_names::Vector{String} = String[],
+    check_connectivity::Bool = true)
 
     # assign timestep_names
     # timestep names are then allocated in a dictionary to map matrix columns
@@ -563,6 +599,7 @@ function PowerFlowData(
         bus_reactivepower_injection_1,
         bus_activepower_withdrawals_1,
         bus_reactivepower_withdrawals_1,
+        Vector{Vector{Float64}}(),
         bus_type,
         bus_magnitude_1,
         bus_angles_1,
@@ -571,5 +608,6 @@ function PowerFlowData(
         setdiff(1:n_buses, aux_network_matrix.ref_bus_positions),
         power_network_matrix,
         aux_network_matrix,
+        Vector{Set{Int}}(),
     )
 end
diff --git a/src/PowerFlows.jl b/src/PowerFlows.jl
index 3f00c08e..85000007 100644
--- a/src/PowerFlows.jl
+++ b/src/PowerFlows.jl
@@ -1,7 +1,7 @@
 module PowerFlows
 
 export solve_powerflow
-export solve_powerflow!
+export solve_ac_powerflow!
 export PowerFlowData
 export DCPowerFlow
 export ACPowerFlow
@@ -14,7 +14,6 @@ import PowerSystems
 import LinearAlgebra
 import NLsolve
 import SparseArrays
-import SharedArrays: SharedArray
 import InfrastructureSystems
 import PowerNetworkMatrices
 import SparseArrays: SparseMatrixCSC
@@ -28,6 +27,8 @@ include("definitions.jl")
 include("powerflow_types.jl")
 include("PowerFlowData.jl")
 include("solve_dc_powerflow.jl")
+include("ac_power_flow.jl")
+include("ac_power_flow_jacobian.jl")
 include("nlsolve_ac_powerflow.jl")
 include("post_processing.jl")
 
diff --git a/src/ac_power_flow.jl b/src/ac_power_flow.jl
new file mode 100644
index 00000000..bb5c3b3a
--- /dev/null
+++ b/src/ac_power_flow.jl
@@ -0,0 +1,160 @@
+const ACPowerFlowData = PowerFlowData{
+    PNM.Ybus{
+        Tuple{Vector{Int64}, Vector{Int64}},
+        Tuple{Dict{Int64, Int64}, Dict{Int64, Int64}},
+    },
+    Nothing,
+}
+
+struct PolarPowerFlow{F, D}
+    Pf::F
+    data::D
+    residual::Vector{Float64}
+    P_net::Vector{Float64}
+    Q_net::Vector{Float64}
+    x0::Vector{Float64}
+end
+
+function _calculate_x0(n::Int,
+    bus_types::Vector{PSY.ACBusTypes},
+    bus_angles::Matrix{Float64},
+    bus_magnitude::Matrix{Float64},
+    bus_activepower_injection::Matrix{Float64},
+    bus_reactivepower_injection::Matrix{Float64},
+    bus_activepower_withdrawals::Matrix{Float64},
+    bus_reactivepower_withdrawals::Matrix{Float64})
+    n_buses = length(bus_types)
+    x0 = Vector{Float64}(undef, 2 * n_buses)
+    for i_n in 1:n
+        state_variable_count = 1
+        for (ix, b) in enumerate(bus_types)
+            if b == PSY.ACBusTypes.REF
+                x0[state_variable_count, i_n] =
+                    bus_activepower_injection[ix, i_n] -
+                    bus_activepower_withdrawals[ix, i_n]
+                x0[state_variable_count + 1, i_n] =
+                    bus_reactivepower_injection[ix, i_n] -
+                    bus_reactivepower_withdrawals[ix, i_n]
+                state_variable_count += 2
+            elseif b == PSY.ACBusTypes.PV
+                x0[state_variable_count, i_n] =
+                    bus_reactivepower_injection[ix, i_n] -
+                    bus_reactivepower_withdrawals[ix, i_n]
+                x0[state_variable_count + 1, i_n] = bus_angles[ix, i_n]
+                state_variable_count += 2
+            elseif b == PSY.ACBusTypes.PQ
+                x0[state_variable_count, i_n] = bus_magnitude[ix, i_n]
+                x0[state_variable_count + 1, i_n] = bus_angles[ix, i_n]
+                state_variable_count += 2
+            else
+                throw(ArgumentError("$b not recognized as a bustype"))
+            end
+        end
+        @assert state_variable_count - 1 == n_buses * 2
+    end
+    return x0
+end
+
+function PolarPowerFlow(data::ACPowerFlowData)
+    n_buses = length(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]
+        Q_net[ix] =
+            data.bus_reactivepower_injection[ix] - data.bus_reactivepower_withdrawals[ix]
+    end
+    x0 = _calculate_x0(1,
+        data.bus_type,
+        data.bus_angles,
+        data.bus_magnitude,
+        data.bus_activepower_injection,
+        data.bus_reactivepower_injection,
+        data.bus_activepower_withdrawals,
+        data.bus_reactivepower_withdrawals)
+    pf_function =
+        (res::Vector{Float64}, X::Vector{Float64}) -> polar_pf!(res, X, P_net, Q_net, data)
+    res = similar(x0)
+    pf_function(res, x0)
+
+    if sum(res) > 10 * (n_buses * 2)
+        _, ix = findmax(res)
+        bus_no = data.bus_lookup[ix]
+        @warn "Initial guess provided results in a large initial residual. Largest residual at bus $bus_no"
+    end
+
+    return PolarPowerFlow(
+        pf_function,
+        data,
+        res,
+        P_net,
+        Q_net,
+        x0,
+    )
+end
+
+function (pf::PolarPowerFlow)(res::Vector{Float64}, x::Vector{Float64})
+    pf.Pf(pf.residual, x)
+    copyto!(res, pf.residual)
+    return
+end
+
+function (pf::PolarPowerFlow)(x::Vector{Float64})
+    pf.Pf(pf.residual, x)
+    return
+end
+
+function polar_pf!(
+    F::Vector{Float64},
+    X::Vector{Float64},
+    P_net::Vector{Float64},
+    Q_net::Vector{Float64},
+    data::ACPowerFlowData,
+)
+    Yb = data.power_network_matrix.data
+    n_buses = length(data.bus_type)
+    for (ix, b) in enumerate(data.bus_type)
+        if b == PSY.ACBusTypes.REF
+            # When bustype == REFERENCE PSY.Bus, state variables are Active and Reactive Power Generated
+            P_net[ix] = X[2 * ix - 1] - data.bus_activepower_withdrawals[ix]
+            Q_net[ix] = X[2 * ix] - data.bus_reactivepower_withdrawals[ix]
+        elseif b == PSY.ACBusTypes.PV
+            # When bustype == PV PSY.Bus, state variables are Reactive Power Generated and Voltage Angle
+            Q_net[ix] = X[2 * ix - 1] - data.bus_reactivepower_withdrawals[ix]
+            data.bus_angles[ix] = X[2 * ix]
+        elseif b == PSY.ACBusTypes.PQ
+            # When bustype == PQ PSY.Bus, state variables are Voltage Magnitude and Voltage Angle
+            data.bus_magnitude[ix] = X[2 * ix - 1]
+            data.bus_angles[ix] = X[2 * ix]
+        end
+    end
+
+    Vm = data.bus_magnitude
+    θ = data.bus_angles
+    # F is active and reactive power balance equations at all buses
+    for ix_f in 1:n_buses
+        S_re = -P_net[ix_f]
+        S_im = -Q_net[ix_f]
+        for ix_t in data.neighbors[ix_f]
+            gb = real(Yb[ix_f, ix_t])
+            bb = imag(Yb[ix_f, ix_t])
+            if ix_f == ix_t
+                S_re += Vm[ix_f] * Vm[ix_t] * gb
+                S_im += -Vm[ix_f] * Vm[ix_t] * bb
+            else
+                S_re +=
+                    Vm[ix_f] *
+                    Vm[ix_t] *
+                    (gb * cos(θ[ix_f] - θ[ix_t]) + bb * sin(θ[ix_f] - θ[ix_t]))
+                S_im +=
+                    Vm[ix_f] *
+                    Vm[ix_t] *
+                    (gb * sin(θ[ix_f] - θ[ix_t]) - bb * cos(θ[ix_f] - θ[ix_t]))
+            end
+        end
+        F[2 * ix_f - 1] = S_re
+        F[2 * ix_f] = S_im
+    end
+    return
+end
diff --git a/src/ac_power_flow_jacobian.jl b/src/ac_power_flow_jacobian.jl
new file mode 100644
index 00000000..542346c1
--- /dev/null
+++ b/src/ac_power_flow_jacobian.jl
@@ -0,0 +1,220 @@
+struct PolarPowerFlowJacobian{F,
+    T <: SparseArrays.SparseMatrixCSC{Float64, Int32},
+}
+    Jf::F
+    Jv::T
+end
+
+function (Jacobianf::PolarPowerFlowJacobian)(
+    J::SparseArrays.SparseMatrixCSC{Float64, Int32},
+    x::Vector{Float64},
+)
+    Jacobianf.Jf(Jacobianf.Jv, x)
+    copyto!(J, Jacobianf.Jv)
+    return
+end
+
+function PolarPowerFlowJacobian(data::ACPowerFlowData, x0::Vector{Float64})
+    j0 = _create_jacobian_matrix_structure(data)
+    F0 =
+        (J::SparseArrays.SparseMatrixCSC{Float64, Int32}, x::Vector{Float64}) ->
+            jsp!(J, x, data)
+    F0(j0, x0)
+    return PolarPowerFlowJacobian(F0, j0)
+end
+
+function _create_jacobian_matrix_structure(data::ACPowerFlowData)
+    #Create Jacobian structure
+    J0_I = Int32[]
+    J0_J = Int32[]
+    J0_V = Float64[]
+
+    for ix_f in eachindex(data.bus_type)
+        F_ix_f_r = 2 * ix_f - 1
+        F_ix_f_i = 2 * ix_f
+        for ix_t in data.neighbors[ix_f]
+            X_ix_t_fst = 2 * ix_t - 1
+            X_ix_t_snd = 2 * ix_t
+            nb = data.bus_type[ix_t]
+            #Set to 0.0 only on connected buses
+            if nb == PSY.ACBusTypes.REF
+                if ix_f == ix_t
+                    #Active PF w/r Local Active Power
+                    push!(J0_I, F_ix_f_r)
+                    push!(J0_J, X_ix_t_fst)
+                    push!(J0_V, 0.0)
+                    #Rective PF w/r Local Reactive Power
+                    push!(J0_I, F_ix_f_i)
+                    push!(J0_J, X_ix_t_snd)
+                    push!(J0_V, 0.0)
+                end
+            elseif nb == PSY.ACBusTypes.PV
+                #Active PF w/r Angle
+                push!(J0_I, F_ix_f_r)
+                push!(J0_J, X_ix_t_snd)
+                push!(J0_V, 0.0)
+                #Reactive PF w/r Angle
+                push!(J0_I, F_ix_f_i)
+                push!(J0_J, X_ix_t_snd)
+                push!(J0_V, 0.0)
+                if ix_f == ix_t
+                    #Reactive PF w/r Local Reactive Power
+                    push!(J0_I, F_ix_f_i)
+                    push!(J0_J, X_ix_t_fst)
+                    push!(J0_V, 0.0)
+                end
+            elseif nb == PSY.ACBusTypes.PQ
+                #Active PF w/r VoltageMag
+                push!(J0_I, F_ix_f_r)
+                push!(J0_J, X_ix_t_fst)
+                push!(J0_V, 0.0)
+                #Reactive PF w/r VoltageMag
+                push!(J0_I, F_ix_f_i)
+                push!(J0_J, X_ix_t_fst)
+                push!(J0_V, 0.0)
+                #Active PF w/r Angle
+                push!(J0_I, F_ix_f_r)
+                push!(J0_J, X_ix_t_snd)
+                push!(J0_V, 0.0)
+                #Reactive PF w/r Angle
+                push!(J0_I, F_ix_f_i)
+                push!(J0_J, X_ix_t_snd)
+                push!(J0_V, 0.0)
+            end
+        end
+    end
+    J0 = SparseArrays.sparse(J0_I, J0_J, J0_V)
+    return J0
+end
+
+function jsp!(
+    J::SparseArrays.SparseMatrixCSC{Float64, Int32},
+    ::Vector{Float64},
+    data::ACPowerFlowData,
+)
+    Yb = data.power_network_matrix.data
+    Vm = data.bus_magnitude
+    θ = data.bus_angles
+    for (ix_f, b) in enumerate(data.bus_type)
+        F_ix_f_r = 2 * ix_f - 1
+        F_ix_f_i = 2 * ix_f
+
+        for ix_t in data.neighbors[ix_f]
+            X_ix_t_fst = 2 * ix_t - 1
+            X_ix_t_snd = 2 * ix_t
+            nb = data.bus_type[ix_t]
+            if nb == PSY.ACBusTypes.REF
+                # State variables are Active and Reactive Power Generated
+                # F[2*i-1] := p[i] = p_flow[i] + p_load[i] - x[2*i-1]
+                # F[2*i] := q[i] = q_flow[i] + q_load[i] - x[2*i]
+                # x does not appears in p_flow and q_flow
+                if ix_f == ix_t
+                    J[F_ix_f_r, X_ix_t_fst] = -1.0
+                    J[F_ix_f_i, X_ix_t_snd] = -1.0
+                end
+            elseif nb == PSY.ACBusTypes.PV
+                # State variables are Reactive Power Generated and Voltage Angle
+                # F[2*i-1] := p[i] = p_flow[i] + p_load[i] - p_gen[i]
+                # F[2*i] := q[i] = q_flow[i] + q_load[i] - x[2*i]
+                # x[2*i] (associated with q_gen) does not appear in q_flow
+                # x[2*i] (associated with q_gen) does not appear in the active power balance
+                if ix_f == ix_t
+                    #Jac: Reactive PF against local active power
+                    J[F_ix_f_i, X_ix_t_fst] = -1.0
+                    #Jac: Active PF against same Angle: θ[ix_f] =  θ[ix_t]
+                    J[F_ix_f_r, X_ix_t_snd] =
+                        Vm[ix_f] * sum(
+                            Vm[k] * (
+                                real(Yb[ix_f, k]) * -sin(θ[ix_f] - θ[k]) +
+                                imag(Yb[ix_f, k]) * cos(θ[ix_f] - θ[k])
+                            ) for k in data.neighbors[ix_f] if k != ix_f
+                        )
+                    #Jac: Reactive PF against same Angle: θ[ix_f] = θ[ix_t]
+                    J[F_ix_f_i, X_ix_t_snd] =
+                        Vm[ix_f] * sum(
+                            Vm[k] * (
+                                real(Yb[ix_f, k]) * cos(θ[ix_f] - θ[k]) -
+                                imag(Yb[ix_f, k]) * -sin(θ[ix_f] - θ[k])
+                            ) for k in data.neighbors[ix_f] if k != ix_f
+                        )
+                else
+                    g_ij = real(Yb[ix_f, ix_t])
+                    b_ij = imag(Yb[ix_f, ix_t])
+                    #Jac: Active PF against other angles θ[ix_t]
+                    J[F_ix_f_r, X_ix_t_snd] =
+                        Vm[ix_f] *
+                        Vm[ix_t] *
+                        (g_ij * sin(θ[ix_f] - θ[ix_t]) + b_ij * -cos(θ[ix_f] - θ[ix_t]))
+                    #Jac: Reactive PF against other angles θ[ix_t]
+                    J[F_ix_f_i, X_ix_t_snd] =
+                        Vm[ix_f] *
+                        Vm[ix_t] *
+                        (g_ij * -cos(θ[ix_f] - θ[ix_t]) - b_ij * sin(θ[ix_f] - θ[ix_t]))
+                end
+            elseif nb == PSY.ACBusTypes.PQ
+                # State variables are Voltage Magnitude and Voltage Angle
+                # Everything appears in everything
+                if ix_f == ix_t
+                    #Jac: Active PF against same voltage magnitude Vm[ix_f]
+                    J[F_ix_f_r, X_ix_t_fst] =
+                        2 * real(Yb[ix_f, ix_t]) * Vm[ix_f] + sum(
+                            Vm[k] * (
+                                real(Yb[ix_f, k]) * cos(θ[ix_f] - θ[k]) +
+                                imag(Yb[ix_f, k]) * sin(θ[ix_f] - θ[k])
+                            ) for k in data.neighbors[ix_f] if k != ix_f
+                        )
+                    #Jac: Active PF against same angle θ[ix_f]
+                    J[F_ix_f_r, X_ix_t_snd] =
+                        Vm[ix_f] * sum(
+                            Vm[k] * (
+                                real(Yb[ix_f, k]) * -sin(θ[ix_f] - θ[k]) +
+                                imag(Yb[ix_f, k]) * cos(θ[ix_f] - θ[k])
+                            ) for k in data.neighbors[ix_f] if k != ix_f
+                        )
+
+                    #Jac: Reactive PF against same voltage magnitude Vm[ix_f]
+                    J[F_ix_f_i, X_ix_t_fst] =
+                        -2 * imag(Yb[ix_f, ix_t]) * Vm[ix_f] + sum(
+                            Vm[k] * (
+                                real(Yb[ix_f, k]) * sin(θ[ix_f] - θ[k]) -
+                                imag(Yb[ix_f, k]) * cos(θ[ix_f] - θ[k])
+                            ) for k in data.neighbors[ix_f] if k != ix_f
+                        )
+                    #Jac: Reactive PF against same angle θ[ix_f]
+                    J[F_ix_f_i, X_ix_t_snd] =
+                        Vm[ix_f] * sum(
+                            Vm[k] * (
+                                real(Yb[ix_f, k]) * cos(θ[ix_f] - θ[k]) -
+                                imag(Yb[ix_f, k]) * -sin(θ[ix_f] - θ[k])
+                            ) for k in data.neighbors[ix_f] if k != ix_f
+                        )
+                else
+                    g_ij = real(Yb[ix_f, ix_t])
+                    b_ij = imag(Yb[ix_f, ix_t])
+                    #Jac: Active PF w/r to different voltage magnitude Vm[ix_t]
+                    J[F_ix_f_r, X_ix_t_fst] =
+                        Vm[ix_f] *
+                        (g_ij * cos(θ[ix_f] - θ[ix_t]) + b_ij * sin(θ[ix_f] - θ[ix_t]))
+                    #Jac: Active PF w/r to different angle θ[ix_t]
+                    J[F_ix_f_r, X_ix_t_snd] =
+                        Vm[ix_f] *
+                        Vm[ix_t] *
+                        (g_ij * sin(θ[ix_f] - θ[ix_t]) + b_ij * -cos(θ[ix_f] - θ[ix_t]))
+
+                    #Jac: Reactive PF w/r to different voltage magnitude Vm[ix_t]
+                    J[F_ix_f_i, X_ix_t_fst] =
+                        Vm[ix_f] *
+                        (g_ij * sin(θ[ix_f] - θ[ix_t]) - b_ij * cos(θ[ix_f] - θ[ix_t]))
+                    #Jac: Reactive PF w/r to different angle θ[ix_t]
+                    J[F_ix_f_i, X_ix_t_snd] =
+                        Vm[ix_f] *
+                        Vm[ix_t] *
+                        (g_ij * -cos(θ[ix_f] - θ[ix_t]) - b_ij * sin(θ[ix_f] - θ[ix_t]))
+                end
+            else
+                error("Undefined Conditional")
+            end
+        end
+    end
+    return
+end
diff --git a/src/common.jl b/src/common.jl
index 1c326d74..6ef897b8 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -54,12 +54,33 @@ function _get_withdrawals!(
         !PSY.get_available(l) && continue
         bus = PSY.get_bus(l)
         bus_ix = bus_lookup[PSY.get_number(bus)]
-        bus_activepower_withdrawals[bus_ix] += PSY.get_active_power(l)
-        bus_reactivepower_withdrawals[bus_ix] += PSY.get_reactive_power(l)
+        bus_activepower_withdrawals[bus_ix] += get_total_p(l)
+        bus_reactivepower_withdrawals[bus_ix] += get_total_q(l)
     end
     return
 end
 
+# TODO: Might need changes if we have SwitchedAdmittances
+function _get_reactive_power_bound!(
+    bus_reactivepower_bounds::Vector{Vector{Float64}},
+    bus_lookup::Dict{Int, Int},
+    sys::PSY.System)
+    sources = PSY.get_components(d -> !isa(d, PSY.ElectricLoad), PSY.StaticInjection, sys)
+    for source in sources
+        !PSY.get_available(source) && continue
+        bus = PSY.get_bus(source)
+        bus_ix = bus_lookup[PSY.get_number(bus)]
+        reactive_power_limits = PSY.get_reactive_power_limits(source)
+        if reactive_power_limits !== nothing
+            bus_reactivepower_bounds[bus_ix][1] += min(0, reactive_power_limits.min)
+            bus_reactivepower_bounds[bus_ix][2] += max(0, reactive_power_limits.max)
+        else
+            @warn("Reactive Power Bounds at Bus $(PSY.get_name(bus)) set to (-Inf, Inf)")
+            bus_reactivepower_bounds[bus_ix][1] = -Inf
+            bus_reactivepower_bounds[bus_ix][2] = Inf
+        end
+    end
+end
 ##############################################################################
 # Matrix Methods #############################################################
 
diff --git a/src/definitions.jl b/src/definitions.jl
index a59f1b27..03af2999 100644
--- a/src/definitions.jl
+++ b/src/definitions.jl
@@ -5,5 +5,8 @@ const MINIMAL_ACCEPTABLE_NLSOLVE_F_TOLERANCE = :1e-3
 const MAX_INIT_RESIDUAL = 10.0
 const MAX_NLSOLVE_INTERATIONS = 10
 const BOUNDS_TOLERANCE = 1e-6
+const MAX_REACTIVE_POWER_ITERATIONS = 10
 
 const ISAPPROX_ZERO_TOLERANCE = 1e-6
+
+const AC_PF_KW = []
diff --git a/src/nlsolve_ac_powerflow.jl b/src/nlsolve_ac_powerflow.jl
index 747864ce..c2d1ed20 100644
--- a/src/nlsolve_ac_powerflow.jl
+++ b/src/nlsolve_ac_powerflow.jl
@@ -1,16 +1,14 @@
-
 """
 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 solving using Finite Differences Method (instead of using analytic Jacobian)
-by setting finite_diff = true.
 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`.
@@ -30,16 +28,20 @@ Arguments available for `nlsolve`:
 solve_powerflow!(sys)
 # Passing NLsolve arguments
 solve_powerflow!(sys, method=:newton)
-# Using Finite Differences
-solve_powerflow!(sys, finite_diff=true)
 ```
 """
-function solve_powerflow!(system::PSY.System; finite_diff = false, kwargs...)
+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")
-    res = _solve_powerflow(system, finite_diff; kwargs...)
+    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),
+    )
+    res = _solve_powerflow!(data, check_reactive_power_limits; kwargs...)
     if res.f_converged
         write_powerflow_solution!(system, res.zero)
         @info("PowerFlow solve converged, the results have been stored in the system")
@@ -62,24 +64,28 @@ Returns the results in a dictionary of dataframes.
 res = solve_powerflow(sys)
 # Passing NLsolve arguments
 res = solve_powerflow(sys, method=:newton)
-# Using Finite Differences
-res = solve_powerflow(sys, finite_diff=true)
 ```
 """
 function solve_powerflow(
-    ::ACPowerFlow,
+    pf::ACPowerFlow,
     system::PSY.System;
-    finite_diff = false,
     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")
-    res = _solve_powerflow(system, finite_diff; kwargs...)
+    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(ACPowerFlow(), system, res.zero)
+        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
@@ -89,333 +95,62 @@ function solve_powerflow(
     return res.f_converged
 end
 
-function _solve_powerflow(system::PSY.System, finite_diff::Bool; kwargs...)
-    ybus_kw = [:check_connectivity, :connectivity_method]
-    ybus_kwargs = (k for k in kwargs if first(k) ∈ ybus_kw)
-    kwargs = (k for k in kwargs if first(k) ∉ ybus_kw)
-
-    buses = sort!(collect(PSY.get_components(PSY.Bus, system)); by = x -> PSY.get_number(x))
-    N_BUS = length(buses)
-
-    # assumes the ordering in YPSY.Bus is the same as in the buses.
-    Yb = PNM.Ybus(system; ybus_kwargs...).data
-    a = collect(1:N_BUS)
-    I, J, V = SparseArrays.findnz(Yb)
-    neighbors = [Set{Int}([i]) for i in 1:N_BUS]
-    for nz in eachindex(V)
-        push!(neighbors[I[nz]], J[nz])
-        push!(neighbors[J[nz]], I[nz])
-    end
-    x0 = zeros(N_BUS * 2)
-
-    #Create Jacobian structure
-    J0_I = Int[]
-    J0_J = Int[]
-    J0_V = Float64[]
-
-    for ix_f in a
-        F_ix_f_r = 2 * ix_f - 1
-        F_ix_f_i = 2 * ix_f
-
-        for ix_t in neighbors[ix_f]
-            X_ix_t_fst = 2 * ix_t - 1
-            X_ix_t_snd = 2 * ix_t
-            b = PSY.get_bustype(buses[ix_t])
-            #Set to 0.0 only on connected buses
-            if b == PSY.ACBusTypes.REF
-                if ix_f == ix_t
-                    #Active PF w/r Local Active Power
-                    push!(J0_I, F_ix_f_r)
-                    push!(J0_J, X_ix_t_fst)
-                    push!(J0_V, 0.0)
-                    #Rective PF w/r Local Reactive Power
-                    push!(J0_I, F_ix_f_i)
-                    push!(J0_J, X_ix_t_snd)
-                    push!(J0_V, 0.0)
-                end
-            elseif b == PSY.ACBusTypes.PV
-                #Active PF w/r Angle
-                push!(J0_I, F_ix_f_r)
-                push!(J0_J, X_ix_t_snd)
-                push!(J0_V, 0.0)
-                #Reactive PF w/r Angle
-                push!(J0_I, F_ix_f_i)
-                push!(J0_J, X_ix_t_snd)
-                push!(J0_V, 0.0)
-                if ix_f == ix_t
-                    #Reactive PF w/r Local Reactive Power
-                    push!(J0_I, F_ix_f_i)
-                    push!(J0_J, X_ix_t_fst)
-                    push!(J0_V, 0.0)
-                end
-            elseif b == PSY.ACBusTypes.PQ
-                #Active PF w/r VoltageMag
-                push!(J0_I, F_ix_f_r)
-                push!(J0_J, X_ix_t_fst)
-                push!(J0_V, 0.0)
-                #Reactive PF w/r VoltageMag
-                push!(J0_I, F_ix_f_i)
-                push!(J0_J, X_ix_t_fst)
-                push!(J0_V, 0.0)
-                #Active PF w/r Angle
-                push!(J0_I, F_ix_f_r)
-                push!(J0_J, X_ix_t_snd)
-                push!(J0_V, 0.0)
-                #Reactive PF w/r Angle
-                push!(J0_I, F_ix_f_i)
-                push!(J0_J, X_ix_t_snd)
-                push!(J0_V, 0.0)
-            end
-        end
-    end
-    J0 = SparseArrays.sparse(J0_I, J0_J, J0_V)
-
-    # Use vectors to cache data for closure
-    # These should be read only
-    P_GEN_BUS = fill(0.0, N_BUS)
-    Q_GEN_BUS = fill(0.0, N_BUS)
-    P_LOAD_BUS = fill(0.0, N_BUS)
-    Q_LOAD_BUS = fill(0.0, N_BUS)
-    P_net = fill(0.0, N_BUS)
-    Q_net = fill(0.0, N_BUS)
-    Vm = fill(0.0, N_BUS)
-    θ = fill(0.0, N_BUS)
-
-    state_variable_count = 1
-    sources =
-        PSY.get_components(d -> !isa(d, PSY.ElectricLoad), PSY.StaticInjection, system)
-
-    for (ix, b) in enumerate(buses)
-        bus_angle = PSY.get_angle(b)::Float64
-        Vm[ix] = bus_voltage_magnitude = PSY.get_magnitude(b)::Float64
-        P_GEN_BUS[ix] = 0.0
-        Q_GEN_BUS[ix] = 0.0
-        PSY.get_ext(b)["neighbors"] = neighbors[ix]
-        for gen in sources
-            !PSY.get_available(gen) && continue
-            if gen.bus == b
-                P_GEN_BUS[ix] += PSY.get_active_power(gen)
-                Q_GEN_BUS[ix] += PSY.get_reactive_power(gen)
-            end
-        end
-
-        P_LOAD_BUS[ix], Q_LOAD_BUS[ix] = _get_load_data(system, b)
-        if b.bustype == PSY.ACBusTypes.REF
-            injection_components = PSY.get_components(
-                d -> PSY.get_available(d) && PSY.get_bus(d) == b,
-                PSY.StaticInjection,
-                system,
-            )
-            isempty(injection_components) && throw(
-                IS.ConflictingInputsError(
-                    "The slack bus does not have any injection component. Power Flow can not proceed",
-                ),
-            )
-            θ[ix] = PSY.get_angle(b)::Float64
-            x0[state_variable_count] = P_GEN_BUS[ix]
-            x0[state_variable_count + 1] = Q_GEN_BUS[ix]
-            state_variable_count += 2
-        elseif b.bustype == PSY.ACBusTypes.PV
-            x0[state_variable_count] = Q_GEN_BUS[ix]
-            x0[state_variable_count + 1] = bus_angle
-            state_variable_count += 2
-        elseif b.bustype == PSY.ACBusTypes.PQ
-            x0[state_variable_count] = bus_voltage_magnitude
-            x0[state_variable_count + 1] = bus_angle
-            state_variable_count += 2
+function _check_q_limit_bounds!(data::ACPowerFlowData, zero::Vector{Float64})
+    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
-            throw(ArgumentError("PSY.Bustype not recognized"))
-        end
-    end
-
-    @assert state_variable_count - 1 == N_BUS * 2
-    bus_types = PSY.get_bustype.(buses)
-    function pf!(F::Vector{Float64}, X::Vector{Float64})
-        for (ix, b) in enumerate(bus_types)
-            if b == PSY.ACBusTypes.REF
-                # When bustype == REFERENCE PSY.Bus, state variables are Active and Reactive Power Generated
-                P_net[ix] = X[2 * ix - 1] - P_LOAD_BUS[ix]
-                Q_net[ix] = X[2 * ix] - Q_LOAD_BUS[ix]
-            elseif b == PSY.ACBusTypes.PV
-                # When bustype == PV PSY.Bus, state variables are Reactive Power Generated and Voltage Angle
-                P_net[ix] = P_GEN_BUS[ix] - P_LOAD_BUS[ix]
-                Q_net[ix] = X[2 * ix - 1] - Q_LOAD_BUS[ix]
-                θ[ix] = X[2 * ix]
-            elseif b == PSY.ACBusTypes.PQ
-                # When bustype == PQ PSY.Bus, state variables are Voltage Magnitude and Voltage Angle
-                P_net[ix] = P_GEN_BUS[ix] - P_LOAD_BUS[ix]
-                Q_net[ix] = Q_GEN_BUS[ix] - Q_LOAD_BUS[ix]
-                Vm[ix] = X[2 * ix - 1]
-                θ[ix] = X[2 * ix]
-            end
+            continue
         end
 
-        # F is active and reactive power balance equations at all buses
-        for ix_f in a
-            S_re = -P_net[ix_f]
-            S_im = -Q_net[ix_f]
-            for ix_t in neighbors[ix_f]
-                gb = real(Yb[ix_f, ix_t])
-                bb = imag(Yb[ix_f, ix_t])
-                if ix_f == ix_t
-                    S_re += Vm[ix_f] * Vm[ix_t] * gb
-                    S_im += -Vm[ix_f] * Vm[ix_t] * bb
-                else
-                    S_re +=
-                        Vm[ix_f] *
-                        Vm[ix_t] *
-                        (gb * cos(θ[ix_f] - θ[ix_t]) + bb * sin(θ[ix_f] - θ[ix_t]))
-                    S_im +=
-                        Vm[ix_f] *
-                        Vm[ix_t] *
-                        (gb * sin(θ[ix_f] - θ[ix_t]) - bb * cos(θ[ix_f] - θ[ix_t]))
-                end
-            end
-            F[2 * ix_f - 1] = S_re
-            F[2 * ix_f] = S_im
+        if Q_gen <= data.bus_reactivepower_bounds[ix][1]
+            @info "Bus $(bus_names[ix]) changed to PSY.ACBusTypes.PQ"
+            within_limits = false
+            data.bus_type[ix] = PSY.ACBusTypes.PQ
+            data.bus_reactivepower_injection[ix] = data.bus_reactivepower_bounds[ix][1]
+        elseif Q_gen >= data.bus_reactivepower_bounds[ix][2]
+            @info "Bus $(bus_names[ix]) changed to PSY.ACBusTypes.PQ"
+            within_limits = false
+            data.bus_type[ix] = PSY.ACBusTypes.PQ
+            data.bus_reactivepower_injection[ix] = data.bus_reactivepower_bounds[ix][2]
+        else
+            @debug "Within Limits"
         end
     end
+    return within_limits
+end
 
-    function jsp!(J::SparseArrays.SparseMatrixCSC{Float64, Int}, X::Vector{Float64})
-        for ix_f in a
-            F_ix_f_r = 2 * ix_f - 1
-            F_ix_f_i = 2 * ix_f
-
-            for ix_t in neighbors[ix_f]
-                X_ix_t_fst = 2 * ix_t - 1
-                X_ix_t_snd = 2 * ix_t
-                b = bus_types[ix_t]
-
-                if b == PSY.ACBusTypes.REF
-                    # State variables are Active and Reactive Power Generated
-                    # F[2*i-1] := p[i] = p_flow[i] + p_load[i] - x[2*i-1]
-                    # F[2*i] := q[i] = q_flow[i] + q_load[i] - x[2*i]
-                    # x does not appears in p_flow and q_flow
-                    if ix_f == ix_t
-                        J[F_ix_f_r, X_ix_t_fst] = -1.0
-                        J[F_ix_f_i, X_ix_t_snd] = -1.0
-                    end
-                elseif b == PSY.ACBusTypes.PV
-                    # State variables are Reactive Power Generated and Voltage Angle
-                    # F[2*i-1] := p[i] = p_flow[i] + p_load[i] - p_gen[i]
-                    # F[2*i] := q[i] = q_flow[i] + q_load[i] - x[2*i]
-                    # x[2*i] (associated with q_gen) does not appear in q_flow
-                    # x[2*i] (associated with q_gen) does not appear in the active power balance
-                    if ix_f == ix_t
-                        #Jac: Reactive PF against local active power
-                        J[F_ix_f_i, X_ix_t_fst] = -1.0
-                        #Jac: Active PF against same Angle: θ[ix_f] =  θ[ix_t]
-                        J[F_ix_f_r, X_ix_t_snd] =
-                            Vm[ix_f] * sum(
-                                Vm[k] * (
-                                    real(Yb[ix_f, k]) * -sin(θ[ix_f] - θ[k]) +
-                                    imag(Yb[ix_f, k]) * cos(θ[ix_f] - θ[k])
-                                ) for k in neighbors[ix_f] if k != ix_f
-                            )
-                        #Jac: Reactive PF against same Angle: θ[ix_f] = θ[ix_t]
-                        J[F_ix_f_i, X_ix_t_snd] =
-                            Vm[ix_f] * sum(
-                                Vm[k] * (
-                                    real(Yb[ix_f, k]) * cos(θ[ix_f] - θ[k]) -
-                                    imag(Yb[ix_f, k]) * -sin(θ[ix_f] - θ[k])
-                                ) for k in neighbors[ix_f] if k != ix_f
-                            )
-                    else
-                        g_ij = real(Yb[ix_f, ix_t])
-                        b_ij = imag(Yb[ix_f, ix_t])
-                        #Jac: Active PF against other angles θ[ix_t]
-                        J[F_ix_f_r, X_ix_t_snd] =
-                            Vm[ix_f] *
-                            Vm[ix_t] *
-                            (g_ij * sin(θ[ix_f] - θ[ix_t]) + b_ij * -cos(θ[ix_f] - θ[ix_t]))
-                        #Jac: Reactive PF against other angles θ[ix_t]
-                        J[F_ix_f_i, X_ix_t_snd] =
-                            Vm[ix_f] *
-                            Vm[ix_t] *
-                            (g_ij * -cos(θ[ix_f] - θ[ix_t]) - b_ij * sin(θ[ix_f] - θ[ix_t]))
-                    end
-                elseif b == PSY.ACBusTypes.PQ
-                    # State variables are Voltage Magnitude and Voltage Angle
-                    # Everything appears in everything
-                    if ix_f == ix_t
-                        #Jac: Active PF against same voltage magnitude Vm[ix_f]
-                        J[F_ix_f_r, X_ix_t_fst] =
-                            2 * real(Yb[ix_f, ix_t]) * Vm[ix_f] + sum(
-                                Vm[k] * (
-                                    real(Yb[ix_f, k]) * cos(θ[ix_f] - θ[k]) +
-                                    imag(Yb[ix_f, k]) * sin(θ[ix_f] - θ[k])
-                                ) for k in neighbors[ix_f] if k != ix_f
-                            )
-                        #Jac: Active PF against same angle θ[ix_f]
-                        J[F_ix_f_r, X_ix_t_snd] =
-                            Vm[ix_f] * sum(
-                                Vm[k] * (
-                                    real(Yb[ix_f, k]) * -sin(θ[ix_f] - θ[k]) +
-                                    imag(Yb[ix_f, k]) * cos(θ[ix_f] - θ[k])
-                                ) for k in neighbors[ix_f] if k != ix_f
-                            )
-
-                        #Jac: Reactive PF against same voltage magnitude Vm[ix_f]
-                        J[F_ix_f_i, X_ix_t_fst] =
-                            -2 * imag(Yb[ix_f, ix_t]) * Vm[ix_f] + sum(
-                                Vm[k] * (
-                                    real(Yb[ix_f, k]) * sin(θ[ix_f] - θ[k]) -
-                                    imag(Yb[ix_f, k]) * cos(θ[ix_f] - θ[k])
-                                ) for k in neighbors[ix_f] if k != ix_f
-                            )
-                        #Jac: Reactive PF against same angle θ[ix_f]
-                        J[F_ix_f_i, X_ix_t_snd] =
-                            Vm[ix_f] * sum(
-                                Vm[k] * (
-                                    real(Yb[ix_f, k]) * cos(θ[ix_f] - θ[k]) -
-                                    imag(Yb[ix_f, k]) * -sin(θ[ix_f] - θ[k])
-                                ) for k in neighbors[ix_f] if k != ix_f
-                            )
-                    else
-                        g_ij = real(Yb[ix_f, ix_t])
-                        b_ij = imag(Yb[ix_f, ix_t])
-                        #Jac: Active PF w/r to different voltage magnitude Vm[ix_t]
-                        J[F_ix_f_r, X_ix_t_fst] =
-                            Vm[ix_f] *
-                            (g_ij * cos(θ[ix_f] - θ[ix_t]) + b_ij * sin(θ[ix_f] - θ[ix_t]))
-                        #Jac: Active PF w/r to different angle θ[ix_t]
-                        J[F_ix_f_r, X_ix_t_snd] =
-                            Vm[ix_f] *
-                            Vm[ix_t] *
-                            (g_ij * sin(θ[ix_f] - θ[ix_t]) + b_ij * -cos(θ[ix_f] - θ[ix_t]))
-
-                        #Jac: Reactive PF w/r to different voltage magnitude Vm[ix_t]
-                        J[F_ix_f_i, X_ix_t_fst] =
-                            Vm[ix_f] *
-                            (g_ij * sin(θ[ix_f] - θ[ix_t]) - b_ij * cos(θ[ix_f] - θ[ix_t]))
-                        #Jac: Reactive PF w/r to different angle θ[ix_t]
-                        J[F_ix_f_i, X_ix_t_snd] =
-                            Vm[ix_f] *
-                            Vm[ix_t] *
-                            (g_ij * -cos(θ[ix_f] - θ[ix_t]) - b_ij * sin(θ[ix_f] - θ[ix_t]))
-                    end
-                else
-                    error("Undefined Conditional")
+function _solve_powerflow!(
+    data::ACPowerFlowData,
+    check_reactive_power_limits;
+    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)
+                    return res
                 end
+            else
+                return res
             end
         end
+    else
+        return _nlsolve_powerflow(data; nlsolve_kwargs...)
     end
+end
 
-    res = similar(x0)
-    pf!(res, x0)
-    if sum(res) > 10 * (N_BUS * 2)
-        _, ix = findmax(res)
-        bus_no = PSY.get_number(buses[(ix + 1) ÷ 2])
-        @warn "Initial guess provided results in a large initial residual. Largest residual at bus $bus_no"
-    end
+function _nlsolve_powerflow(data::ACPowerFlowData; nlsolve_kwargs...)
+    pf = PolarPowerFlow(data)
+    J = PowerFlows.PolarPowerFlowJacobian(data, pf.x0)
 
-    if finite_diff
-        res = NLsolve.nlsolve(pf!, x0; kwargs...)
-    else
-        F0 = similar(x0)
-        df = NLsolve.OnceDifferentiable(pf!, jsp!, x0, F0, J0)
-        res = NLsolve.nlsolve(df, x0; kwargs...)
+    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/post_processing.jl b/src/post_processing.jl
index 250c8e8a..9e759f19 100644
--- a/src/post_processing.jl
+++ b/src/post_processing.jl
@@ -596,6 +596,7 @@ dictionary will therefore feature just one key linked to one DataFrame.
 function write_results(
     ::ACPowerFlow,
     sys::PSY.System,
+    data::ACPowerFlowData,
     result::Vector{Float64},
 )
     @info("Voltages are exported in pu. Powers are exported in MW/MVAr.")
@@ -603,7 +604,6 @@ function write_results(
     N_BUS = length(buses)
     bus_map = Dict(buses .=> 1:N_BUS)
     sys_basepower = PSY.get_base_power(sys)
-    sources = PSY.get_components(d -> !isa(d, PSY.ElectricLoad), PSY.StaticInjection, sys)
     Vm_vect = fill(0.0, N_BUS)
     θ_vect = fill(0.0, N_BUS)
     P_gen_vect = fill(0.0, N_BUS)
@@ -611,36 +611,27 @@ function write_results(
     P_load_vect = fill(0.0, N_BUS)
     Q_load_vect = fill(0.0, N_BUS)
 
-    for (ix, bus) in enumerate(buses)
-        P_load_vect[ix], Q_load_vect[ix] = _get_load_data(sys, bus) .* sys_basepower
-        P_admittance, Q_admittance = _get_fixed_admittance_power(sys, bus, result, ix)
+    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
+        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 bus.bustype == PSY.ACBusTypes.REF
-            Vm_vect[ix] = PSY.get_magnitude(bus)
-            θ_vect[ix] = PSY.get_angle(bus)
+        if bustype == PSY.ACBusTypes.REF
+            Vm_vect[ix] = data.bus_magnitude[ix]
+            θ_vect[ix] = data.bus_angles[ix]
             P_gen_vect[ix] = result[2 * ix - 1] * sys_basepower
             Q_gen_vect[ix] = result[2 * ix] * sys_basepower
-        elseif bus.bustype == PSY.ACBusTypes.PV
-            Vm_vect[ix] = PSY.get_magnitude(bus)
+        elseif bustype == PSY.ACBusTypes.PV
+            Vm_vect[ix] = data.bus_magnitude[ix]
             θ_vect[ix] = result[2 * ix]
-            for gen in sources
-                !PSY.get_available(gen) && continue
-                if gen.bus == bus
-                    P_gen_vect[ix] += PSY.get_active_power(gen) * sys_basepower
-                end
-            end
+            P_gen_vect[ix] = data.bus_activepower_injection[ix] * sys_basepower
             Q_gen_vect[ix] = result[2 * ix - 1] * sys_basepower
-        elseif bus.bustype == PSY.ACBusTypes.PQ
+        elseif bustype == PSY.ACBusTypes.PQ
             Vm_vect[ix] = result[2 * ix - 1]
             θ_vect[ix] = result[2 * ix]
-            for gen in sources
-                !PSY.get_available(gen) && continue
-                if gen.bus == bus
-                    P_gen_vect[ix] += PSY.get_active_power(gen) * sys_basepower
-                    Q_gen_vect[ix] += PSY.get_reactive_power(gen) * sys_basepower
-                end
-            end
+            P_gen_vect[ix] = data.bus_activepower_injection[ix] * sys_basepower
+            Q_gen_vect[ix] = data.bus_reactivepower_injection[ix] * sys_basepower
         end
     end
 
diff --git a/src/powerflow_types.jl b/src/powerflow_types.jl
index 6c75bcac..d8b6cffe 100644
--- a/src/powerflow_types.jl
+++ b/src/powerflow_types.jl
@@ -1,4 +1,7 @@
-struct ACPowerFlow end
+@kwdef struct ACPowerFlow
+    check_reactive_power_limits::Bool = false
+end
+
 struct DCPowerFlow end
 struct PTDFDCPowerFlow end
 struct vPTDFDCPowerFlow end
diff --git a/src/solve_dc_powerflow.jl b/src/solve_dc_powerflow.jl
index a18f7e4c..19e8f9fb 100644
--- a/src/solve_dc_powerflow.jl
+++ b/src/solve_dc_powerflow.jl
@@ -25,50 +25,6 @@ const ABAPowerFlowData = PowerFlowData{
         Tuple{Dict{Int64, Int64}, Dict{String, Int64}}},
 }
 
-# # method based on ABA and BA matrices
-# function _solve_powerflows!(
-#     data::ABAPowerFlowData,
-# )
-#     # get net injections
-#     power_injection = data.bus_activepower_injection - data.bus_activepower_withdrawals
-#     # save angles and power flows
-#     data.bus_angles[data.valid_ix, :] .=
-#         data.power_network_matrix.K \ @view power_injection[data.valid_ix, :]
-#     data.branch_flow_values .= data.aux_network_matrix.data' * data.bus_angles
-#     return
-# end
-
-# # method based on PTDF matrix
-# function _solve_powerflows!(
-#     data::PTDFPowerFlowData,
-# )
-#     # get net power injections
-#     power_injection = data.bus_activepower_injection .- data.bus_activepower_withdrawals
-#     # evaluate flows
-#     data.branch_flow_values .= data.power_network_matrix.data' * power_injection
-#     # evaluate bus angles
-#     p_inj = power_injection[data.valid_ix, :]
-#     data.bus_angles[data.valid_ix, :] .= data.aux_network_matrix.K \ p_inj
-#     return
-# end
-
-# # method based on Virtual PTDF
-# function _solve_powerflows!(
-#     data::vPTDFPowerFlowData,
-# )
-#     # get net power injections
-#     power_injection = data.bus_activepower_injection .- data.bus_activepower_withdrawals
-#     for i in axes(power_injection, 2)
-#         # evaluate flows (next line evaluates both both PTDF rows and line flows)
-#         data.branch_flow_values[:, i] .=
-#             my_mul_mt(data.power_network_matrix, power_injection[:, i])
-#         # evaluate bus angles
-#         p_inj = power_injection[data.valid_ix, i]
-#         data.bus_angles[data.valid_ix, i] .= data.aux_network_matrix.K \ p_inj
-#     end
-#     return
-# end
-
 """
 Evaluates the power flows on each system's branch and updates the PowerFlowData structure.
 
@@ -137,18 +93,18 @@ end
 # SINGLE PERIOD ##############################################################
 
 """
-Evaluates the power flows on each system's branch by means of the PTDF matrix. 
+Evaluates the power flows on each system's branch by means of the PTDF matrix.
 Updates the PowerFlowData structure and returns a dictionary containing a
 DataFrame for the single timestep considered.
-The DataFrame containts the flows and angles related to the information stored 
+The DataFrame containts the flows and angles related to the information stored
 in the PSY.System considered as input.
 
 # Arguments:
 - `::PTDFDCPowerFlow`:
-        use PTDFDCPowerFlow() to evaluate the power flows according to the 
+        use PTDFDCPowerFlow() to evaluate the power flows according to the
         method based on the PTDF matrix
 - `sys::PSY.System`:
-        container gathering the system data used for the evaluation of flows 
+        container gathering the system data used for the evaluation of flows
         and angles.
 """
 function solve_powerflow(
@@ -161,19 +117,19 @@ function solve_powerflow(
 end
 
 """
-Evaluates the power flows on each system's branch by means of the ABA and BA 
-matrices. 
+Evaluates the power flows on each system's branch by means of the ABA and BA
+matrices.
 Updates the PowerFlowData structure and returns a dictionary containing a
 DataFrame for the single timestep considered.
-The DataFrame containts the flows and angles related to the information stored 
+The DataFrame containts the flows and angles related to the information stored
 in the PSY.System considered as input.
 
 # Arguments:
 - `::DCPowerFlow`:
-        use DCPowerFlow() to evaluate the power flows according to the method 
+        use DCPowerFlow() to evaluate the power flows according to the method
         based on the ABA and BA matrices
 - `sys::PSY.System`:
-        container gathering the system data used for the evaluation of flows 
+        container gathering the system data used for the evaluation of flows
         and angles.
 """
 function solve_powerflow(
@@ -186,19 +142,19 @@ function solve_powerflow(
 end
 
 """
-Evaluates the power flows on each system's branch by means of the Virtual PTDF 
+Evaluates the power flows on each system's branch by means of the Virtual PTDF
 matrix.
 Updates the PowerFlowData structure "data" and returns a dictionary containing
 a number of DataFrames equal to the numeber of timestep considered in "data".
-The DataFrame containts the flows and angles related to the information stored 
+The DataFrame containts the flows and angles related to the information stored
 in the PSY.System considered as input.
 
 # Arguments:
 - `::vPTDFDCPowerFlow`:
-        use vPTDFDCPowerFlow() to evaluate the power flows according to the 
+        use vPTDFDCPowerFlow() to evaluate the power flows according to the
         method based on the Virtual PTDF matrix
 - `sys::PSY.System`:
-        container gathering the system data used for the evaluation of flows 
+        container gathering the system data used for the evaluation of flows
         and angles.
 """
 function solve_powerflow(
@@ -234,7 +190,7 @@ function solve_powerflow(
 end
 
 """
-Evaluates the power flows on each system's branch by means of the ABA and BA 
+Evaluates the power flows on each system's branch by means of the ABA and BA
 matrices.
 Updates the PowerFlowData structure "data" and returns a dictionary containing
 a number of DataFrames equal to the numeber of timestep considered in "data".
@@ -256,7 +212,7 @@ function solve_powerflow(
 end
 
 """
-Evaluates the power flows on each system's branch by means of Virtual PTDF 
+Evaluates the power flows on each system's branch by means of Virtual PTDF
 matrices.
 Updates the PowerFlowData structure "data" and returns a dictionary containing
 a number of DataFrames equal to the numeber of timestep considered in "data".
diff --git a/test/Project.toml b/test/Project.toml
index 79a25226..799f9cdd 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -11,7 +11,6 @@ PowerFlows = "94fada2c-fd9a-4e89-8d82-81405f5cb4f6"
 PowerNetworkMatrices = "bed98974-b02a-5e2f-9fe0-a103f5c450dd"
 PowerSystemCaseBuilder = "f00506e0-b84f-492a-93c2-c0a9afc4364e"
 PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"
-SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 638a0c9e..ae803eec 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -26,6 +26,8 @@ const TEST_FILES_DIR = test_file_dir
 const DIFF_INF_TOLERANCE = 1e-4
 const DIFF_L2_TOLERANCE = 1e-3
 
+MAIN_DIR = dirname(@__DIR__)
+
 include("test_utils/psse_results_compare.jl")
 
 LOG_FILE = "power-systems.log"
diff --git a/test/test_dc_powerflow.jl b/test/test_dc_powerflow.jl
index 76f8f6fd..5a487e9c 100644
--- a/test/test_dc_powerflow.jl
+++ b/test/test_dc_powerflow.jl
@@ -1,6 +1,3 @@
-# TODO:
-# - get_active_power(::StandardLoad) not working
-
 @testset "SINGLE PERIOD power flows evaluation: ABA, PTDF, VirtualPTDF" begin
     # get system
     sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
@@ -28,7 +25,7 @@
     ref_bus_angles[valid_ix] = matrix_data \ power_injection[valid_ix]
     ref_flow_values = transpose(aux_network_matrix.data) * ref_bus_angles
 
-    # CASE 1: ABA and BA matrices 
+    # CASE 1: ABA and BA matrices
     solved_data_ABA = solve_powerflow(DCPowerFlow(), sys)
     @test isapprox(
         solved_data_ABA["1"]["flow_results"].P_from_to +
diff --git a/test/test_multiperiod_dc_powerflow.jl b/test/test_multiperiod_dc_powerflow.jl
index 27417eb4..74633e4f 100644
--- a/test/test_multiperiod_dc_powerflow.jl
+++ b/test/test_multiperiod_dc_powerflow.jl
@@ -1,6 +1,4 @@
 @testset "MULTI-PERIOD power flows evaluation: ABA, PTDF, VirtualPTDF" begin
-    MAIN_DIR = dirname(@__DIR__)
-
     # get system
     sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
     injections = CSV.read(
@@ -28,9 +26,9 @@
 
     # create structure for multi-period case
     timesteps = 24
-    data_1 = PowerFlowData(DCPowerFlow(), sys, timesteps)
-    data_2 = PowerFlowData(PTDFDCPowerFlow(), sys, timesteps)
-    data_3 = PowerFlowData(vPTDFDCPowerFlow(), sys, timesteps)
+    data_1 = PowerFlowData(DCPowerFlow(), sys; timesteps = timesteps)
+    data_2 = PowerFlowData(PTDFDCPowerFlow(), sys; timesteps = timesteps)
+    data_3 = PowerFlowData(vPTDFDCPowerFlow(), sys; timesteps = timesteps)
 
     # allocate data from csv
     injs = Matrix(injections)
diff --git a/test/test_nlsolve_powerflow.jl b/test/test_nlsolve_powerflow.jl
index d676d044..7e33b99e 100644
--- a/test/test_nlsolve_powerflow.jl
+++ b/test/test_nlsolve_powerflow.jl
@@ -1,85 +1,64 @@
-result = [
-    2.32551,
-    -0.155293,
-    0.469214,
-    -0.0870457,
-    0.271364,
-    -0.222398,
-    1.01423,
-    -0.179009,
-    1.01724,
-    -0.152972,
-    0.216039,
-    -0.251637,
-    1.05034,
-    -0.231289,
-    0.245388,
-    -0.231289,
-    1.03371,
-    -0.258872,
-    1.03256,
-    -0.262519,
-    1.04748,
-    -0.259143,
-    1.0535,
-    -0.266484,
-    1.04711,
-    -0.267177,
-    1.02131,
-    -0.280381,
-]
-
-p_gen_matpower_3bus = [20.3512373930753, 100.0, 100.0]
-q_gen_matpower_3bus = [45.516916781567232, 10.453799727283879, -31.992561631394636]
-
-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")
-br = get_component(Line, pf_sys5_re, "6")
-PSY.set_x!(br, 20.0)
-PSY.set_r!(br, 2.0)
-
-@testset "NLsolve Power Flow testing" begin
-    # 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_powerflow!(pf_sys5_re; finite_diff = true)
-    )
+@testset "NLsolve Power Flow 14-Bus testing" begin
+    sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
+    result_14 = [
+        2.3255081760423684
+        -0.15529254415401786
+        0.4692141795968793
+        -0.08704571860678871
+        0.27136388547189727
+        -0.22239752522709744
+        1.014232467422514
+        -0.17900878100582837
+        1.0172382900002028
+        -0.15297197376986682
+        0.21603888720954267
+        -0.2516366411330649
+        1.0503438761049335
+        -0.23128945395825606
+        0.2453884476050094
+        -0.23128945395825604
+        1.0337109552890456
+        -0.2588720100456853
+        1.0325606430799592
+        -0.26251860149671896
+        1.0474837090281963
+        -0.25914254888894855
+        1.0535012072934
+        -0.26648433793564014
+        1.0471069158440354
+        -0.2671769024662062
+        1.0213119628726421
+        -0.2803812119374241
+    ]
+    data = PowerFlows.PowerFlowData(ACPowerFlow(), sys; check_connectivity = true)
     #Compare results between finite diff methods and Jacobian method
-    res_finite_diff = solve_powerflow(
-        ACPowerFlow(),
-        PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false);
-        finite_diff = true,
-    )
-    res_jacobian =
-        solve_powerflow(
-            ACPowerFlow(),
-            PSB.build_system(PSB.PSITestSystems,
-                "c_sys14"; add_forecasts = false),
-        )
-    @test LinearAlgebra.norm(
-        res_finite_diff["bus_results"].Vm - res_jacobian["bus_results"].Vm,
-    ) <= 1e-6
-    @test solve_powerflow!(
-        PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false),
-        finite_diff = true,
-        method = :newton,
-    )
+    res1 = PowerFlows._solve_powerflow!(data, false)
+    @test LinearAlgebra.norm(result_14 - res1.zero) <= 1e-6
+    @test solve_ac_powerflow!(sys; method = :newton)
 
+    # Test enforcing the reactive power Limits
+    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
+end
+
+@testset "NLsolve Power Flow 14-Bus Line Configurations" begin
     sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
+    base_res = solve_powerflow(ACPowerFlow(), sys)
     branch = first(PSY.get_components(Line, sys))
     dyn_branch = DynamicBranch(branch)
     add_component!(sys, dyn_branch)
-    @test dyn_pf = solve_powerflow!(sys)
+    @test dyn_pf = solve_ac_powerflow!(sys)
     dyn_pf = solve_powerflow(ACPowerFlow(), sys)
-    @test LinearAlgebra.norm(dyn_pf["bus_results"].Vm - res_jacobian["bus_results"].Vm) <=
+    @test LinearAlgebra.norm(dyn_pf["bus_results"].Vm - base_res["bus_results"].Vm) <=
           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_powerflow!(sys)
+    solve_ac_powerflow!(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)
@@ -90,11 +69,11 @@ PSY.set_r!(br, 2.0)
     res = solve_powerflow(ACPowerFlow(), sys)
     @test res["flow_results"].P_from_to[4] == 0.0
     @test res["flow_results"].P_to_from[4] == 0.0
+end
 
-    # Investigate why this test is failing
-    #sys = PSB.build_system(PSB.PSYTestSystems, "psse_240_parsing_sys")
-    #@test solve_powerflow!(sys)
-
+@testset "NLsolve Power Flow 3-Bus Fixed FixedAdmittance testing" begin
+    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)
@@ -104,6 +83,22 @@ PSY.set_r!(br, 2.0)
     @test isapprox(df["bus_results"].Q_gen, q_gen_matpower_3bus, atol = 1e-4)
 end
 
+@testset "NLsolve Power Flow convergence fail testing" begin
+    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")
+    br = get_component(Line, pf_sys5_re, "6")
+    PSY.set_x!(br, 20.0)
+    PSY.set_r!(br, 2.0)
+
+    # 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)
+    )
+end
+
 @testset "Test 240 Case PSS/e results" begin
     file = joinpath(
         TEST_FILES_DIR,
@@ -119,7 +114,7 @@ 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_powerflow!(system)
+    pf = solve_ac_powerflow!(system)
     @test pf
     pf_result_df = solve_powerflow(ACPowerFlow(), system)
 
@@ -136,3 +131,280 @@ end
     @test sum(q_diff) < DIFF_INF_TOLERANCE * base_power
     @test norm(q_diff, 2) / length(q_diff) < DIFF_L2_TOLERANCE
 end
+
+@testset "Multiple sources at ref" begin
+    sys = System(100.0)
+    b = ACBus(;
+        number = 1,
+        name = "01",
+        bustype = ACBusTypes.REF,
+        angle = 0.0,
+        magnitude = 1.1,
+        voltage_limits = (0.0, 2.0),
+        base_voltage = 230,
+    )
+    add_component!(sys, b)
+
+    #Test two sources with equal and opposite P and Q
+    s1 = Source(;
+        name = "source_1",
+        available = true,
+        bus = b,
+        active_power = 0.5,
+        reactive_power = 0.1,
+        R_th = 1e-5,
+        X_th = 1e-5,
+    )
+    add_component!(sys, s1)
+    s2 = Source(;
+        name = "source_2",
+        available = true,
+        bus = b,
+        active_power = -0.5,
+        reactive_power = -0.1,
+        R_th = 1e-5,
+        X_th = 1e-5,
+    )
+    add_component!(sys, s2)
+    @test solve_ac_powerflow!(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)
+end
+
+@testset "AC PowerFlow with Multiple sources at PV" begin
+    sys = System(100.0)
+    b1 = ACBus(;
+        number = 1,
+        name = "01",
+        bustype = ACBusTypes.REF,
+        angle = 0.0,
+        magnitude = 1.1,
+        voltage_limits = (0.0, 2.0),
+        base_voltage = 230,
+    )
+    add_component!(sys, b1)
+    b2 = ACBus(;
+        number = 2,
+        name = "02",
+        bustype = ACBusTypes.PV,
+        angle = 0.0,
+        magnitude = 1.1,
+        voltage_limits = (0.0, 2.0),
+        base_voltage = 230,
+    )
+    add_component!(sys, b2)
+    a = Arc(; from = b1, to = b2)
+    add_component!(sys, a)
+    l = Line(;
+        name = "l1",
+        available = true,
+        active_power_flow = 0.0,
+        reactive_power_flow = 0.0,
+        arc = a,
+        r = 1e-3,
+        x = 1e-3,
+        b = (from = 0.0, to = 0.0),
+        rate = 0.0,
+        angle_limits = (min = -pi / 2, max = pi / 2),
+    )
+    add_component!(sys, l)
+
+    #Test two sources with equal and opposite P and Q
+    s1 = Source(;
+        name = "source_1",
+        available = true,
+        bus = b1,
+        active_power = 0.5,
+        reactive_power = 0.1,
+        R_th = 1e-5,
+        X_th = 1e-5,
+    )
+    add_component!(sys, s1)
+    s2 = Source(;
+        name = "source_2",
+        available = true,
+        bus = b2,
+        active_power = 0.5,
+        reactive_power = 1.1,
+        R_th = 1e-5,
+        X_th = 1e-5,
+    )
+    add_component!(sys, s2)
+    s3 = Source(;
+        name = "source_3",
+        available = true,
+        bus = b2,
+        active_power = -0.5,
+        reactive_power = -1.1,
+        R_th = 1e-5,
+        X_th = 1e-5,
+    )
+    add_component!(sys, s3)
+
+    @test solve_ac_powerflow!(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!(
+        sys,
+    )
+end
+
+@testset "AC PowerFlow Source + non-source at Ref" begin
+    sys = System(100.0)
+    b = ACBus(;
+        number = 1,
+        name = "01",
+        bustype = ACBusTypes.REF,
+        angle = 0.0,
+        magnitude = 1.1,
+        voltage_limits = (0.0, 2.0),
+        base_voltage = 230,
+    )
+    add_component!(sys, b)
+
+    #Test two sources with equal and opposite P and Q
+    s1 = Source(;
+        name = "source_1",
+        available = true,
+        bus = b,
+        active_power = 0.5,
+        reactive_power = 0.1,
+        R_th = 1e-5,
+        X_th = 1e-5,
+    )
+    add_component!(sys, s1)
+    g1 = ThermalStandard(;
+        name = "init",
+        available = true,
+        status = false,
+        bus = b,
+        active_power = 0.1,
+        reactive_power = 0.1,
+        rating = 0.0,
+        active_power_limits = (min = 0.0, max = 0.0),
+        reactive_power_limits = nothing,
+        ramp_limits = nothing,
+        operation_cost = ThreePartCost(nothing),
+        base_power = 100.0,
+        time_limits = nothing,
+        prime_mover_type = PrimeMovers.OT,
+        fuel = ThermalFuels.OTHER,
+        services = Device[],
+        dynamic_injector = nothing,
+        ext = Dict{String, Any}(),
+        time_series_container = InfrastructureSystems.TimeSeriesContainer(),
+    )
+    add_component!(sys, g1)
+
+    @test solve_ac_powerflow!(sys)
+    @test isapprox(
+        get_active_power(get_component(Source, sys, "source_1")),
+        0.5;
+        atol = 0.001,
+    )
+    @test isapprox(
+        get_reactive_power(get_component(Source, sys, "source_1")),
+        0.1;
+        atol = 0.001,
+    )
+end
+
+@testset "AC PowerFlow Source + non-source at PV" begin
+    sys = System(100.0)
+    b1 = ACBus(;
+        number = 1,
+        name = "01",
+        bustype = ACBusTypes.REF,
+        angle = 0.0,
+        magnitude = 1.1,
+        voltage_limits = (0.0, 2.0),
+        base_voltage = 230,
+    )
+    add_component!(sys, b1)
+    b2 = ACBus(;
+        number = 2,
+        name = "02",
+        bustype = ACBusTypes.PV,
+        angle = 0.0,
+        magnitude = 1.1,
+        voltage_limits = (0.0, 2.0),
+        base_voltage = 230,
+    )
+    add_component!(sys, b2)
+    a = Arc(; from = b1, to = b2)
+    add_component!(sys, a)
+    l = Line(;
+        name = "l1",
+        available = true,
+        active_power_flow = 0.0,
+        reactive_power_flow = 0.0,
+        arc = a,
+        r = 1e-3,
+        x = 1e-3,
+        b = (from = 0.0, to = 0.0),
+        rate = 0.0,
+        angle_limits = (min = -pi / 2, max = pi / 2),
+    )
+    add_component!(sys, l)
+
+    #Test two sources with equal and opposite P and Q
+    s1 = Source(;
+        name = "source_1",
+        available = true,
+        bus = b1,
+        active_power = 0.5,
+        reactive_power = 0.1,
+        R_th = 1e-5,
+        X_th = 1e-5,
+    )
+    add_component!(sys, s1)
+    s2 = Source(;
+        name = "source_2",
+        available = true,
+        bus = b2,
+        active_power = 0.5,
+        reactive_power = 1.1,
+        R_th = 1e-5,
+        X_th = 1e-5,
+    )
+    add_component!(sys, s2)
+    g1 = ThermalStandard(;
+        name = "init",
+        available = true,
+        status = false,
+        bus = b2,
+        active_power = 0.1,
+        reactive_power = 0.1,
+        rating = 0.0,
+        active_power_limits = (min = 0.0, max = 0.0),
+        reactive_power_limits = nothing,
+        ramp_limits = nothing,
+        operation_cost = ThreePartCost(nothing),
+        base_power = 100.0,
+        time_limits = nothing,
+        prime_mover_type = PrimeMovers.OT,
+        fuel = ThermalFuels.OTHER,
+        services = Device[],
+        dynamic_injector = nothing,
+        ext = Dict{String, Any}(),
+        time_series_container = InfrastructureSystems.TimeSeriesContainer(),
+    )
+    add_component!(sys, g1)
+
+    @test solve_ac_powerflow!(sys)
+    @test isapprox(
+        get_active_power(get_component(Source, sys, "source_2")),
+        0.5;
+        atol = 0.001,
+    )
+    @test isapprox(
+        get_reactive_power(get_component(Source, sys, "source_2")),
+        1.1;
+        atol = 0.001,
+    )
+end
diff --git a/test/test_powerflow_data.jl b/test/test_powerflow_data.jl
index 729ff640..1f352156 100644
--- a/test/test_powerflow_data.jl
+++ b/test/test_powerflow_data.jl
@@ -10,7 +10,7 @@ end
     sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)
     timesteps = 24
     # TODO: "multiperiod AC still to implement"
-    PowerFlowData(DCPowerFlow(), sys, timesteps)
-    PowerFlowData(PTDFDCPowerFlow(), sys, timesteps)
-    PowerFlowData(vPTDFDCPowerFlow(), sys, timesteps)
+    PowerFlowData(DCPowerFlow(), sys; timesteps = timesteps)
+    PowerFlowData(PTDFDCPowerFlow(), sys; timesteps = timesteps)
+    PowerFlowData(vPTDFDCPowerFlow(), sys; timesteps = timesteps)
 end
diff --git a/test/test_powerflow_with_sources.jl b/test/test_powerflow_with_sources.jl
deleted file mode 100644
index 444d7f4d..00000000
--- a/test/test_powerflow_with_sources.jl
+++ /dev/null
@@ -1,277 +0,0 @@
-@testset "Multiple sources at ref" begin
-    sys = System(100.0)
-    b = ACBus(;
-        number = 1,
-        name = "01",
-        bustype = ACBusTypes.REF,
-        angle = 0.0,
-        magnitude = 1.1,
-        voltage_limits = (0.0, 2.0),
-        base_voltage = 230,
-    )
-    add_component!(sys, b)
-
-    #Test two sources with equal and opposite P and Q
-    s1 = Source(;
-        name = "source_1",
-        available = true,
-        bus = b,
-        active_power = 0.5,
-        reactive_power = 0.1,
-        R_th = 1e-5,
-        X_th = 1e-5,
-    )
-    add_component!(sys, s1)
-    s2 = Source(;
-        name = "source_2",
-        available = true,
-        bus = b,
-        active_power = -0.5,
-        reactive_power = -0.1,
-        R_th = 1e-5,
-        X_th = 1e-5,
-    )
-    add_component!(sys, s2)
-    @test solve_powerflow!(sys, finite_diff = true)
-
-    #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_powerflow!(sys, finite_diff = true)
-end
-
-@testset "Multiple sources at PV" begin
-    sys = System(100.0)
-    b1 = ACBus(;
-        number = 1,
-        name = "01",
-        bustype = ACBusTypes.REF,
-        angle = 0.0,
-        magnitude = 1.1,
-        voltage_limits = (0.0, 2.0),
-        base_voltage = 230,
-    )
-    add_component!(sys, b1)
-    b2 = ACBus(;
-        number = 2,
-        name = "02",
-        bustype = ACBusTypes.PV,
-        angle = 0.0,
-        magnitude = 1.1,
-        voltage_limits = (0.0, 2.0),
-        base_voltage = 230,
-    )
-    add_component!(sys, b2)
-    a = Arc(; from = b1, to = b2)
-    add_component!(sys, a)
-    l = Line(;
-        name = "l1",
-        available = true,
-        active_power_flow = 0.0,
-        reactive_power_flow = 0.0,
-        arc = a,
-        r = 1e-3,
-        x = 1e-3,
-        b = (from = 0.0, to = 0.0),
-        rate = 0.0,
-        angle_limits = (min = -pi / 2, max = pi / 2),
-    )
-    add_component!(sys, l)
-
-    #Test two sources with equal and opposite P and Q
-    s1 = Source(;
-        name = "source_1",
-        available = true,
-        bus = b1,
-        active_power = 0.5,
-        reactive_power = 0.1,
-        R_th = 1e-5,
-        X_th = 1e-5,
-    )
-    add_component!(sys, s1)
-    s2 = Source(;
-        name = "source_2",
-        available = true,
-        bus = b2,
-        active_power = 0.5,
-        reactive_power = 1.1,
-        R_th = 1e-5,
-        X_th = 1e-5,
-    )
-    add_component!(sys, s2)
-    s3 = Source(;
-        name = "source_3",
-        available = true,
-        bus = b2,
-        active_power = -0.5,
-        reactive_power = -1.1,
-        R_th = 1e-5,
-        X_th = 1e-5,
-    )
-    add_component!(sys, s3)
-
-    @test solve_powerflow!(sys, finite_diff = true)
-
-    #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_powerflow!(
-        sys,
-        finite_diff = true,
-    )
-end
-
-@testset "Source + non-source at Ref" begin
-    sys = System(100.0)
-    b = ACBus(;
-        number = 1,
-        name = "01",
-        bustype = ACBusTypes.REF,
-        angle = 0.0,
-        magnitude = 1.1,
-        voltage_limits = (0.0, 2.0),
-        base_voltage = 230,
-    )
-    add_component!(sys, b)
-
-    #Test two sources with equal and opposite P and Q
-    s1 = Source(;
-        name = "source_1",
-        available = true,
-        bus = b,
-        active_power = 0.5,
-        reactive_power = 0.1,
-        R_th = 1e-5,
-        X_th = 1e-5,
-    )
-    add_component!(sys, s1)
-    g1 = ThermalStandard(;
-        name = "init",
-        available = true,
-        status = false,
-        bus = b,
-        active_power = 0.1,
-        reactive_power = 0.1,
-        rating = 0.0,
-        active_power_limits = (min = 0.0, max = 0.0),
-        reactive_power_limits = nothing,
-        ramp_limits = nothing,
-        operation_cost = ThreePartCost(nothing),
-        base_power = 100.0,
-        time_limits = nothing,
-        prime_mover_type = PrimeMovers.OT,
-        fuel = ThermalFuels.OTHER,
-        services = Device[],
-        dynamic_injector = nothing,
-        ext = Dict{String, Any}(),
-        time_series_container = InfrastructureSystems.TimeSeriesContainer(),
-    )
-    add_component!(sys, g1)
-
-    @test solve_powerflow!(sys, finite_diff = true)
-    @test isapprox(
-        get_active_power(get_component(Source, sys, "source_1")),
-        0.5;
-        atol = 0.001,
-    )
-    @test isapprox(
-        get_reactive_power(get_component(Source, sys, "source_1")),
-        0.1;
-        atol = 0.001,
-    )
-end
-
-@testset "Source + non-source at PV" begin
-    sys = System(100.0)
-    b1 = ACBus(;
-        number = 1,
-        name = "01",
-        bustype = ACBusTypes.REF,
-        angle = 0.0,
-        magnitude = 1.1,
-        voltage_limits = (0.0, 2.0),
-        base_voltage = 230,
-    )
-    add_component!(sys, b1)
-    b2 = ACBus(;
-        number = 2,
-        name = "02",
-        bustype = ACBusTypes.PV,
-        angle = 0.0,
-        magnitude = 1.1,
-        voltage_limits = (0.0, 2.0),
-        base_voltage = 230,
-    )
-    add_component!(sys, b2)
-    a = Arc(; from = b1, to = b2)
-    add_component!(sys, a)
-    l = Line(;
-        name = "l1",
-        available = true,
-        active_power_flow = 0.0,
-        reactive_power_flow = 0.0,
-        arc = a,
-        r = 1e-3,
-        x = 1e-3,
-        b = (from = 0.0, to = 0.0),
-        rate = 0.0,
-        angle_limits = (min = -pi / 2, max = pi / 2),
-    )
-    add_component!(sys, l)
-
-    #Test two sources with equal and opposite P and Q
-    s1 = Source(;
-        name = "source_1",
-        available = true,
-        bus = b1,
-        active_power = 0.5,
-        reactive_power = 0.1,
-        R_th = 1e-5,
-        X_th = 1e-5,
-    )
-    add_component!(sys, s1)
-    s2 = Source(;
-        name = "source_2",
-        available = true,
-        bus = b2,
-        active_power = 0.5,
-        reactive_power = 1.1,
-        R_th = 1e-5,
-        X_th = 1e-5,
-    )
-    add_component!(sys, s2)
-    g1 = ThermalStandard(;
-        name = "init",
-        available = true,
-        status = false,
-        bus = b2,
-        active_power = 0.1,
-        reactive_power = 0.1,
-        rating = 0.0,
-        active_power_limits = (min = 0.0, max = 0.0),
-        reactive_power_limits = nothing,
-        ramp_limits = nothing,
-        operation_cost = ThreePartCost(nothing),
-        base_power = 100.0,
-        time_limits = nothing,
-        prime_mover_type = PrimeMovers.OT,
-        fuel = ThermalFuels.OTHER,
-        services = Device[],
-        dynamic_injector = nothing,
-        ext = Dict{String, Any}(),
-        time_series_container = InfrastructureSystems.TimeSeriesContainer(),
-    )
-    add_component!(sys, g1)
-
-    @test solve_powerflow!(sys, finite_diff = true)
-    @test isapprox(
-        get_active_power(get_component(Source, sys, "source_2")),
-        0.5;
-        atol = 0.001,
-    )
-    @test isapprox(
-        get_reactive_power(get_component(Source, sys, "source_2")),
-        1.1;
-        atol = 0.001,
-    )
-end