diff --git a/test/results/results_eigenvalues.jl b/test/results/results_eigenvalues.jl
index b37b79b4c..2e4edb5ea 100644
--- a/test/results/results_eigenvalues.jl
+++ b/test/results/results_eigenvalues.jl
@@ -971,3 +971,21 @@ test58_eigvals_GateFlag = [
     -0.19849913626541832 + 0.17939599481963514im
     -0.04238523317744444 + 0.0im
 ]
+
+test60_eigvals = [
+    -106.13231871102334 + 0.0im
+    -91.13024819676329 + 0.0im
+    -63.8688081304463 + 0.0im
+    -40.03565106938418 + 0.0im
+    -38.75323547447537 - 8.125818822049656im
+    -38.75323547447537 + 8.125818822049656im
+    -6.104826779980879 + 0.0im
+    -1.395127240875047 - 5.966607970662552im
+    -1.395127240875047 + 5.966607970662552im
+    -0.9449531063097173 + 0.0im
+    -0.48179149983179287 + 0.0im
+    -0.1941504953444676 - 0.18213498194497163im
+    -0.1941504953444676 + 0.18213498194497163im
+    -0.03057660237077197 - 0.30118043163774183im
+    -0.03057660237077197 + 0.30118043163774183im
+]
diff --git a/test/results/results_initial_conditions.jl b/test/results/results_initial_conditions.jl
index d54376ab5..aca059dff 100644
--- a/test/results/results_initial_conditions.jl
+++ b/test/results/results_initial_conditions.jl
@@ -1752,7 +1752,7 @@ test58_x0_init_PowerFlag = Dict(
     ],
 )
 
-test58_x0_init_GateFlag = Dict(
+test60_x0_init = Dict(
     "V_R" => [
         1.05
         1.0197944718502572
@@ -1766,18 +1766,18 @@ test58_x0_init_GateFlag = Dict(
     "generator-102-1" => [
         0.7538967192769663
         0.5555487379562241
-        0.7042452148052647
-        0.7246287886385533
+        0.7042452148052648
+        0.7246287886385532
         0.9158416020737494
         1.0
         1.4986692863524897
         0.044960078577949814
         0.0
-        0.95
-        0.95
+        3.871458709170383e-12
+        3.871569731472846e-12
         0.0
-        0.95
-        0.95
+        3.871569731472846e-12
+        0.7417441860465115
         2.099999999999999
     ],
 )
diff --git a/test/test_case60_wpidhy.jl b/test/test_case60_wpidhy.jl
new file mode 100644
index 000000000..458a303ec
--- /dev/null
+++ b/test/test_case60_wpidhy.jl
@@ -0,0 +1,208 @@
+"""
+Validation PSSE/WPIDHY:
+This case study defines a three bus system with an infinite bus, GENROU+SEXS+PIDGOV and a load.
+The fault drop the line connecting the infinite bus and GENROU
+"""
+
+##################################################
+############### SOLVE PROBLEM ####################
+##################################################
+
+# Define dyr files
+
+names = ["WPIDHY"]
+
+dyr_files = [
+    joinpath(TEST_FILES_DIR, "benchmarks/psse/WPIDHY/ThreeBus_WPIDHY.dyr"),
+]
+
+csv_files = [
+    joinpath(TEST_FILES_DIR, "benchmarks/psse/WPIDHY/WPIDHY_results.csv"),
+]
+
+init_conditions = [test60_x0_init]
+
+eigs_values = [test60_eigvals]
+
+raw_file_dir = joinpath(TEST_FILES_DIR, "benchmarks/psse/WPIDHY/ThreeBusMulti.raw")
+tspan = (0.0, 20.0)
+
+function test_wpidhy_implicit(dyr_file, csv_file, init_cond, eigs_value)
+    path = (joinpath(pwd(), "test-psse-wpidhy"))
+    !isdir(path) && mkdir(path)
+    try
+        sys = System(raw_file_dir, dyr_file)
+        for l in get_components(PSY.StandardLoad, sys)
+            transform_load_to_constant_impedance(l)
+        end
+
+        # Define Simulation Problem
+        sim = Simulation!(
+            ResidualModel,
+            sys, #system
+            path,
+            tspan, #time span
+            BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
+        ) #Type of Fault
+
+        # Test Initial Condition
+        diff_val = [0.0]
+        res = get_init_values_for_comparison(sim)
+        for (k, v) in init_cond
+            diff_val[1] += LinearAlgebra.norm(res[k] - v)
+        end
+
+        @test (diff_val[1] < 1e-3)
+
+        # Obtain small signal results for initial conditions. Testing the simulation reset
+        small_sig = small_signal_analysis(sim)
+        eigs = small_sig.eigenvalues
+        @test small_sig.stable
+
+        # Test Eigenvalues
+        @test LinearAlgebra.norm(eigs - eigs_value) < 1e-3
+
+        # Solve problem
+        @test execute!(sim, IDA(); dtmax = 0.005, saveat = 0.005) ==
+              PSID.SIMULATION_FINALIZED
+        results = read_results(sim)
+
+        # Obtain data for voltage magnitude at bus 102
+        series = get_voltage_magnitude_series(results, 102)
+        t = series[1]
+        V = series[2]
+        # Test Vf, τm and branch series flows with PSSE
+        _, P101_103 = get_activepower_branch_flow(results, "BUS 1-BUS 3-i_1", :from)
+        _, Q101_103 = get_reactivepower_branch_flow(results, "BUS 1-BUS 3-i_1", :from)
+        _, P103_101 = get_activepower_branch_flow(results, "BUS 1-BUS 3-i_1", :to)
+        _, Q103_101 = get_reactivepower_branch_flow(results, "BUS 1-BUS 3-i_1", :to)
+        _, Vf = get_field_voltage_series(results, "generator-102-1")
+        _, τm = get_mechanical_torque_series(results, "generator-102-1")
+
+        # TODO: Get PSSE CSV files and enable tests
+        #M = get_csv_data(csv_file)
+        #t_psse = M[:, 1]
+        #V_psse = M[:, 2]
+        #P101_103_psse = M[:, 3] ./ 100.0 # convert to pu
+        #Q101_103_psse = M[:, 4] ./ 100.0 # convert to pu
+        #P103_101_psse = M[:, 5] ./ 100.0 # convert to pu
+        #Q103_101_psse = M[:, 6] ./ 100.0 # convert to pu
+        #Vf_psse = M[:, 7]
+        #τm_psse = M[:, 8]
+
+        # Test Transient Simulation Results
+        #@test LinearAlgebra.norm(V - V_psse, Inf) <= 1e-2
+        #@test LinearAlgebra.norm(P101_103 - P101_103_psse, Inf) <= 1e-2
+        #@test LinearAlgebra.norm(Q101_103 - Q101_103_psse, Inf) <= 1e-2
+        #@test LinearAlgebra.norm(P103_101 - P103_101_psse, Inf) <= 1e-2
+        #@test LinearAlgebra.norm(Q103_101 - Q103_101_psse, Inf) <= 1e-2
+        #@test LinearAlgebra.norm(Vf - Vf_psse, Inf) <= 1e-2
+        #@test LinearAlgebra.norm(τm - τm_psse, Inf) <= 1e-2
+        #@test LinearAlgebra.norm(t - round.(t_psse, digits = 3)) == 0.0
+    finally
+        @info("removing test files")
+        rm(path; force = true, recursive = true)
+    end
+end
+
+function test_wpidhy_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
+    path = (joinpath(pwd(), "test-psse-wpidhy"))
+    !isdir(path) && mkdir(path)
+    try
+        sys = System(raw_file_dir, dyr_file)
+        for l in get_components(PSY.StandardLoad, sys)
+            transform_load_to_constant_impedance(l)
+        end
+
+        # Define Simulation Problem
+        sim = Simulation!(
+            MassMatrixModel,
+            sys, #system
+            path,
+            tspan, #time span
+            BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
+        ) #Type of Fault
+
+        # Test Initial Condition
+        diff_val = [0.0]
+        res = get_init_values_for_comparison(sim)
+        for (k, v) in init_cond
+            diff_val[1] += LinearAlgebra.norm(res[k] - v)
+        end
+
+        @test (diff_val[1] < 1e-3)
+
+        # Obtain small signal results for initial conditions. Testing the simulation reset
+        small_sig = small_signal_analysis(sim)
+        eigs = small_sig.eigenvalues
+        @test small_sig.stable
+
+        # Test Eigenvalues
+        @test LinearAlgebra.norm(eigs - eigs_value) < 1e-3
+
+        # Solve problem
+        @test execute!(sim, Rodas4(); dtmax = 0.005, saveat = 0.005) ==
+              PSID.SIMULATION_FINALIZED
+        results = read_results(sim)
+
+        # Obtain data for voltage magnitude at bus 102
+        series = get_voltage_magnitude_series(results, 102)
+        t = series[1]
+        V = series[2]
+        # Test Vf, τm and branch series flows with PSSE
+        _, P101_103 = get_activepower_branch_flow(results, "BUS 1-BUS 3-i_1", :from)
+        _, Q101_103 = get_reactivepower_branch_flow(results, "BUS 1-BUS 3-i_1", :from)
+        _, P103_101 = get_activepower_branch_flow(results, "BUS 1-BUS 3-i_1", :to)
+        _, Q103_101 = get_reactivepower_branch_flow(results, "BUS 1-BUS 3-i_1", :to)
+        _, Vf = get_field_voltage_series(results, "generator-102-1")
+        _, τm = get_mechanical_torque_series(results, "generator-102-1")
+
+        # TODO: Get PSSE CSV files and enable tests
+        #M = get_csv_data(csv_file)
+        #t_psse = M[:, 1]
+        #V_psse = M[:, 2]
+        #P101_103_psse = M[:, 3] ./ 100.0 # convert to pu
+        #Q101_103_psse = M[:, 4] ./ 100.0 # convert to pu
+        #P103_101_psse = M[:, 5] ./ 100.0 # convert to pu
+        #Q103_101_psse = M[:, 6] ./ 100.0 # convert to pu
+        #Vf_psse = M[:, 7]
+        #τm_psse = M[:, 8]
+
+        # Test Transient Simulation Results
+        #@test LinearAlgebra.norm(V - V_psse, Inf) <= 1e-2
+        #@test LinearAlgebra.norm(P101_103 - P101_103_psse, Inf) <= 1e-2
+        #@test LinearAlgebra.norm(Q101_103 - Q101_103_psse, Inf) <= 1e-2
+        #@test LinearAlgebra.norm(P103_101 - P103_101_psse, Inf) <= 1e-2
+        #@test LinearAlgebra.norm(Q103_101 - Q103_101_psse, Inf) <= 1e-2
+        #@test LinearAlgebra.norm(Vf - Vf_psse, Inf) <= 1e-2
+        #@test LinearAlgebra.norm(τm - τm_psse, Inf) <= 1e-2
+        #@test LinearAlgebra.norm(t - round.(t_psse, digits = 3)) == 0.0
+    finally
+        @info("removing test files")
+        rm(path; force = true, recursive = true)
+    end
+end
+
+@testset "Test 60 WPIDHY ResidualModel" begin
+    for (ix, name) in enumerate(names)
+        @testset "$(name)" begin
+            dyr_file = dyr_files[ix]
+            csv_file = csv_files[ix]
+            init_cond = init_conditions[ix]
+            eigs_value = eigs_values[ix]
+            test_wpidhy_implicit(dyr_file, csv_file, init_cond, eigs_value)
+        end
+    end
+end
+
+@testset "Test 60 WPIDHY MassMatrixModel" begin
+    for (ix, name) in enumerate(names)
+        @testset "$(name)" begin
+            dyr_file = dyr_files[ix]
+            csv_file = csv_files[ix]
+            init_cond = init_conditions[ix]
+            eigs_value = eigs_values[ix]
+            test_wpidhy_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
+        end
+    end
+end