diff --git a/Project.toml b/Project.toml
index 33a058b8..63c022a2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "PowerFlows"
 uuid = "94fada2c-fd9a-4e89-8d82-81405f5cb4f6"
 authors = ["Jose Daniel Lara <jdlara@berkeley.edu>"]
-version = "0.4.0"
+version = "0.5.0"
 
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
@@ -17,6 +17,6 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 DataFrames = "1"
 InfrastructureSystems = "1"
 NLsolve = "4"
-PowerNetworkMatrices = "^0.8"
-PowerSystems = "2"
+PowerNetworkMatrices = "^0.9"
+PowerSystems = "^3"
 julia = "^1.6"
diff --git a/src/PowerFlowData.jl b/src/PowerFlowData.jl
index e751eb11..325ad17b 100644
--- a/src/PowerFlowData.jl
+++ b/src/PowerFlowData.jl
@@ -21,7 +21,7 @@ flows and angles, as well as these ones.
 - `bus_reactivepower_withdrawals::Matrix{Float64}`:
         "(b, t)" matrix containing the bus reactive power withdrawals. b: 
         number of buses, t: number of time period.
-- `bus_type::Vector{PSY.BusTypes}`:
+- `bus_type::Vector{PSY.ACBusTypes}`:
         vector containing type of buses present in the system, ordered 
         according to "bus_lookup".
 - `bus_magnitude::Matrix{Float64}`:
@@ -55,7 +55,7 @@ struct PowerFlowData{M <: PNM.PowerNetworkMatrix, N, S <: Union{String, Char}}
     bus_reactivepower_injection::Matrix{Float64}
     bus_activepower_withdrawals::Matrix{Float64}
     bus_reactivepower_withdrawals::Matrix{Float64}
-    bus_type::Vector{PSY.BusTypes}
+    bus_type::Vector{PSY.ACBusTypes}
     bus_magnitude::Matrix{Float64}
     bus_angles::Matrix{Float64}
     branch_flow_values::Matrix{Float64}
@@ -119,7 +119,7 @@ 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))
-    bus_type = Vector{PSY.BusTypes}(undef, n_buses)
+    bus_type = Vector{PSY.ACBusTypes}(undef, n_buses)
     bus_angles = zeros(Float64, n_buses)
     bus_magnitude = zeros(Float64, n_buses)
     temp_bus_map = Dict{Int, String}(
@@ -130,7 +130,7 @@ function PowerFlowData(
         bus_name = temp_bus_map[bus_no]
         bus = PSY.get_component(PSY.Bus, sys, bus_name)
         bus_type[ix] = PSY.get_bustype(bus)
-        if bus_type[ix] == PSY.BusTypes.REF
+        if bus_type[ix] == PSY.ACBusTypes.REF
             bus_angles[ix] = 0.0
         else
             bus_angles[ix] = PSY.get_angle(bus)
@@ -242,18 +242,18 @@ function PowerFlowData(
 
     bus_lookup = aux_network_matrix.lookup[1]
     branch_lookup = aux_network_matrix.lookup[2]
-    bus_type = Vector{PSY.BusTypes}(undef, n_buses)
+    bus_type = Vector{PSY.ACBusTypes}(undef, n_buses)
     bus_angles = zeros(Float64, n_buses)
     bus_magnitude = zeros(Float64, n_buses)
     temp_bus_map = Dict{Int, String}(
-        PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.Bus, sys)
+        PSY.get_number(b) => PSY.get_name(b) for b in PSY.get_components(PSY.ACBus, sys)
     )
 
     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)
-        if bus_type[ix] == PSY.BusTypes.REF
+        if bus_type[ix] == PSY.ACBusTypes.REF
             bus_angles[ix] = 0.0
         else
             bus_angles[ix] = PSY.get_angle(bus)
@@ -293,6 +293,8 @@ 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
@@ -368,7 +370,7 @@ function PowerFlowData(
 
     bus_lookup = power_network_matrix.lookup[1]
     branch_lookup = power_network_matrix.lookup[2]
-    bus_type = Vector{PSY.BusTypes}(undef, n_buses)
+    bus_type = Vector{PSY.ACBusTypes}(undef, n_buses)
     bus_angles = zeros(Float64, n_buses)
     bus_magnitude = zeros(Float64, n_buses)
     temp_bus_map = Dict{Int, String}(
@@ -379,7 +381,7 @@ function PowerFlowData(
         bus_name = temp_bus_map[bus_no]
         bus = PSY.get_component(PSY.Bus, sys, bus_name)
         bus_type[ix] = PSY.get_bustype(bus)
-        if bus_type[ix] == PSY.BusTypes.REF
+        if bus_type[ix] == PSY.ACBusTypes.REF
             bus_angles[ix] = 0.0
         else
             bus_angles[ix] = PSY.get_angle(bus)
@@ -494,7 +496,7 @@ function PowerFlowData(
 
     bus_lookup = power_network_matrix.lookup[2]
     branch_lookup = power_network_matrix.lookup[1]
-    bus_type = Vector{PSY.BusTypes}(undef, n_buses)
+    bus_type = Vector{PSY.ACBusTypes}(undef, n_buses)
     bus_angles = zeros(Float64, n_buses)
     bus_magnitude = zeros(Float64, n_buses)
     temp_bus_map = Dict{Int, String}(
@@ -505,7 +507,7 @@ function PowerFlowData(
         bus_name = temp_bus_map[bus_no]
         bus = PSY.get_component(PSY.Bus, sys, bus_name)
         bus_type[ix] = PSY.get_bustype(bus)
-        if bus_type[ix] == PSY.BusTypes.REF
+        if bus_type[ix] == PSY.ACBusTypes.REF
             bus_angles[ix] = 0.0
         else
             bus_angles[ix] = PSY.get_angle(bus)
diff --git a/src/nlsolve_ac_powerflow.jl b/src/nlsolve_ac_powerflow.jl
index ecdcbd09..747864ce 100644
--- a/src/nlsolve_ac_powerflow.jl
+++ b/src/nlsolve_ac_powerflow.jl
@@ -122,7 +122,7 @@ function _solve_powerflow(system::PSY.System, finite_diff::Bool; kwargs...)
             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.BusTypes.REF
+            if b == PSY.ACBusTypes.REF
                 if ix_f == ix_t
                     #Active PF w/r Local Active Power
                     push!(J0_I, F_ix_f_r)
@@ -133,7 +133,7 @@ function _solve_powerflow(system::PSY.System, finite_diff::Bool; kwargs...)
                     push!(J0_J, X_ix_t_snd)
                     push!(J0_V, 0.0)
                 end
-            elseif b == PSY.BusTypes.PV
+            elseif b == PSY.ACBusTypes.PV
                 #Active PF w/r Angle
                 push!(J0_I, F_ix_f_r)
                 push!(J0_J, X_ix_t_snd)
@@ -148,7 +148,7 @@ function _solve_powerflow(system::PSY.System, finite_diff::Bool; kwargs...)
                     push!(J0_J, X_ix_t_fst)
                     push!(J0_V, 0.0)
                 end
-            elseif b == PSY.BusTypes.PQ
+            elseif b == PSY.ACBusTypes.PQ
                 #Active PF w/r VoltageMag
                 push!(J0_I, F_ix_f_r)
                 push!(J0_J, X_ix_t_fst)
@@ -200,7 +200,7 @@ function _solve_powerflow(system::PSY.System, finite_diff::Bool; kwargs...)
         end
 
         P_LOAD_BUS[ix], Q_LOAD_BUS[ix] = _get_load_data(system, b)
-        if b.bustype == PSY.BusTypes.REF
+        if b.bustype == PSY.ACBusTypes.REF
             injection_components = PSY.get_components(
                 d -> PSY.get_available(d) && PSY.get_bus(d) == b,
                 PSY.StaticInjection,
@@ -215,11 +215,11 @@ function _solve_powerflow(system::PSY.System, finite_diff::Bool; kwargs...)
             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.BusTypes.PV
+        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.BusTypes.PQ
+        elseif b.bustype == PSY.ACBusTypes.PQ
             x0[state_variable_count] = bus_voltage_magnitude
             x0[state_variable_count + 1] = bus_angle
             state_variable_count += 2
@@ -232,16 +232,16 @@ function _solve_powerflow(system::PSY.System, finite_diff::Bool; kwargs...)
     bus_types = PSY.get_bustype.(buses)
     function pf!(F::Vector{Float64}, X::Vector{Float64})
         for (ix, b) in enumerate(bus_types)
-            if b == PSY.BusTypes.REF
+            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.BusTypes.PV
+            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.BusTypes.PQ
+            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]
@@ -286,7 +286,7 @@ function _solve_powerflow(system::PSY.System, finite_diff::Bool; kwargs...)
                 X_ix_t_snd = 2 * ix_t
                 b = bus_types[ix_t]
 
-                if b == PSY.BusTypes.REF
+                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]
@@ -295,7 +295,7 @@ function _solve_powerflow(system::PSY.System, finite_diff::Bool; kwargs...)
                         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.BusTypes.PV
+                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]
@@ -334,7 +334,7 @@ function _solve_powerflow(system::PSY.System, finite_diff::Bool; kwargs...)
                             Vm[ix_t] *
                             (g_ij * -cos(θ[ix_f] - θ[ix_t]) - b_ij * sin(θ[ix_f] - θ[ix_t]))
                     end
-                elseif b == PSY.BusTypes.PQ
+                elseif b == PSY.ACBusTypes.PQ
                     # State variables are Voltage Magnitude and Voltage Angle
                     # Everything appears in everything
                     if ix_f == ix_t
diff --git a/src/post_processing.jl b/src/post_processing.jl
index 8def3034..250c8e8a 100644
--- a/src/post_processing.jl
+++ b/src/post_processing.jl
@@ -146,7 +146,11 @@ function _get_fixed_admittance_power(
         !PSY.get_available(l) && continue
         if (l.bus == b)
             Vm_squared =
-                b.bustype == PSY.BusTypes.PQ ? result[2 * ix - 1]^2 : PSY.get_magnitude(b)^2
+                if b.bustype == PSY.ACBusTypes.PQ
+                    result[2 * ix - 1]^2
+                else
+                    PSY.get_magnitude(b)^2
+                end
             active_power += Vm_squared * real(PSY.get_Y(l))
             reactive_power -= Vm_squared * imag(PSY.get_Y(l))
         end
@@ -433,15 +437,15 @@ function write_powerflow_solution!(sys::PSY.System, result::Vector{Float64})
     )
 
     for (ix, bus) in buses
-        if bus.bustype == PSY.BusTypes.REF
+        if bus.bustype == PSY.ACBusTypes.REF
             P_gen = result[2 * ix - 1]
             Q_gen = result[2 * ix]
             _power_redistribution_ref(sys, P_gen, Q_gen, bus)
-        elseif bus.bustype == PSY.BusTypes.PV
+        elseif bus.bustype == PSY.ACBusTypes.PV
             Q_gen = result[2 * ix - 1]
             bus.angle = result[2 * ix]
             _reactive_power_redistribution_pv(sys, Q_gen, bus)
-        elseif bus.bustype == PSY.BusTypes.PQ
+        elseif bus.bustype == PSY.ACBusTypes.PQ
             Vm = result[2 * ix - 1]
             θ = result[2 * ix]
             PSY.set_magnitude!(bus, Vm)
@@ -460,6 +464,8 @@ end
 
 # returns list of branches names and buses numbers: PTDF and virtual PTDF case
 function _get_branches_buses(data::Union{PTDFPowerFlowData, vPTDFPowerFlowData})
+    PNM.get_branch_ax(data.power_network_matrix)
+    PNM.get_bus_ax(data.power_network_matrix)
     return PNM.get_branch_ax(data.power_network_matrix),
     PNM.get_bus_ax(data.power_network_matrix)
 end
@@ -610,12 +616,12 @@ function write_results(
         P_admittance, Q_admittance = _get_fixed_admittance_power(sys, bus, result, ix)
         P_load_vect[ix] += P_admittance * sys_basepower
         Q_load_vect[ix] += Q_admittance * sys_basepower
-        if bus.bustype == PSY.BusTypes.REF
+        if bus.bustype == PSY.ACBusTypes.REF
             Vm_vect[ix] = PSY.get_magnitude(bus)
             θ_vect[ix] = PSY.get_angle(bus)
             P_gen_vect[ix] = result[2 * ix - 1] * sys_basepower
             Q_gen_vect[ix] = result[2 * ix] * sys_basepower
-        elseif bus.bustype == PSY.BusTypes.PV
+        elseif bus.bustype == PSY.ACBusTypes.PV
             Vm_vect[ix] = PSY.get_magnitude(bus)
             θ_vect[ix] = result[2 * ix]
             for gen in sources
@@ -625,7 +631,7 @@ function write_results(
                 end
             end
             Q_gen_vect[ix] = result[2 * ix - 1] * sys_basepower
-        elseif bus.bustype == PSY.BusTypes.PQ
+        elseif bus.bustype == PSY.ACBusTypes.PQ
             Vm_vect[ix] = result[2 * ix - 1]
             θ_vect[ix] = result[2 * ix]
             for gen in sources
diff --git a/test/test_multiperiod_dc_powerflow.jl b/test/test_multiperiod_dc_powerflow.jl
index d667a9a9..27417eb4 100644
--- a/test/test_multiperiod_dc_powerflow.jl
+++ b/test/test_multiperiod_dc_powerflow.jl
@@ -1,11 +1,28 @@
 @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("test_data/c_sys14_injections.csv", DataFrame; header = 0)
-    withdrawals = CSV.read("test_data/c_sys14_withdrawals.csv", DataFrame; header = 0)
-    flows = CSV.read("test_data/c_sys14_flows.csv", DataFrame; header = 0)
-    angles = CSV.read("test_data/c_sys14_angles.csv", DataFrame; header = 0)
+    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,
+    )
 
     ##############################################################################
 
diff --git a/test/test_powerflow_with_sources.jl b/test/test_powerflow_with_sources.jl
index 8e351448..444d7f4d 100644
--- a/test/test_powerflow_with_sources.jl
+++ b/test/test_powerflow_with_sources.jl
@@ -1,9 +1,9 @@
 @testset "Multiple sources at ref" begin
     sys = System(100.0)
-    b = Bus(;
+    b = ACBus(;
         number = 1,
         name = "01",
-        bustype = BusTypes.REF,
+        bustype = ACBusTypes.REF,
         angle = 0.0,
         magnitude = 1.1,
         voltage_limits = (0.0, 2.0),
@@ -43,20 +43,20 @@ end
 
 @testset "Multiple sources at PV" begin
     sys = System(100.0)
-    b1 = Bus(;
+    b1 = ACBus(;
         number = 1,
         name = "01",
-        bustype = BusTypes.REF,
+        bustype = ACBusTypes.REF,
         angle = 0.0,
         magnitude = 1.1,
         voltage_limits = (0.0, 2.0),
         base_voltage = 230,
     )
     add_component!(sys, b1)
-    b2 = Bus(;
+    b2 = ACBus(;
         number = 2,
         name = "02",
-        bustype = BusTypes.PV,
+        bustype = ACBusTypes.PV,
         angle = 0.0,
         magnitude = 1.1,
         voltage_limits = (0.0, 2.0),
@@ -123,10 +123,10 @@ end
 
 @testset "Source + non-source at Ref" begin
     sys = System(100.0)
-    b = Bus(;
+    b = ACBus(;
         number = 1,
         name = "01",
-        bustype = BusTypes.REF,
+        bustype = ACBusTypes.REF,
         angle = 0.0,
         magnitude = 1.1,
         voltage_limits = (0.0, 2.0),
@@ -159,7 +159,7 @@ end
         operation_cost = ThreePartCost(nothing),
         base_power = 100.0,
         time_limits = nothing,
-        prime_mover = PrimeMovers.OT,
+        prime_mover_type = PrimeMovers.OT,
         fuel = ThermalFuels.OTHER,
         services = Device[],
         dynamic_injector = nothing,
@@ -183,20 +183,20 @@ end
 
 @testset "Source + non-source at PV" begin
     sys = System(100.0)
-    b1 = Bus(;
+    b1 = ACBus(;
         number = 1,
         name = "01",
-        bustype = BusTypes.REF,
+        bustype = ACBusTypes.REF,
         angle = 0.0,
         magnitude = 1.1,
         voltage_limits = (0.0, 2.0),
         base_voltage = 230,
     )
     add_component!(sys, b1)
-    b2 = Bus(;
+    b2 = ACBus(;
         number = 2,
         name = "02",
-        bustype = BusTypes.PV,
+        bustype = ACBusTypes.PV,
         angle = 0.0,
         magnitude = 1.1,
         voltage_limits = (0.0, 2.0),
@@ -254,7 +254,7 @@ end
         operation_cost = ThreePartCost(nothing),
         base_power = 100.0,
         time_limits = nothing,
-        prime_mover = PrimeMovers.OT,
+        prime_mover_type = PrimeMovers.OT,
         fuel = ThermalFuels.OTHER,
         services = Device[],
         dynamic_injector = nothing,