From c97be587766b519617fd7083f1c30faafda90966 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebasti=C3=A1n=20M?= <104728742+SebastianManriqueM@users.noreply.github.com> Date: Wed, 18 Sep 2024 19:49:45 -0600 Subject: [PATCH] Add wpidhy (#391) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Add WPIDHY governor model * Fix typo at the input of the regulator * Add initialization for WPIDHY Turbine Governor * Fix typo at initialization * Add files for testing in 3 machine system. Sign of the reg parameter was corrected in .dyr file, it must be negative (See PSSE documentation) * Formatter * Add post processing * Add Test * fix rebase issue * formatter * change type --------- Co-authored-by: Rodrigo Henríquez-Auba --- .../generator_components/init_tg.jl | 97 ++++++++ src/models/common_controls.jl | 2 +- src/models/generator_models/tg_models.jl | 105 +++++++++ src/post_processing/post_proc_generator.jl | 31 +++ test/benchmarks/psse/WPIDHY/ThreeBusMulti.raw | 32 +++ .../psse/WPIDHY/ThreeBus_WPIDHY.dyr | 10 + test/results/results_eigenvalues.jl | 18 ++ test/results/results_initial_conditions.jl | 30 +++ test/test_case60_wpidhy.jl | 208 ++++++++++++++++++ 9 files changed, 532 insertions(+), 1 deletion(-) create mode 100644 test/benchmarks/psse/WPIDHY/ThreeBusMulti.raw create mode 100644 test/benchmarks/psse/WPIDHY/ThreeBus_WPIDHY.dyr create mode 100644 test/test_case60_wpidhy.jl diff --git a/src/initialization/generator_components/init_tg.jl b/src/initialization/generator_components/init_tg.jl index 178834af5..e0c4b1028 100644 --- a/src/initialization/generator_components/init_tg.jl +++ b/src/initialization/generator_components/init_tg.jl @@ -406,3 +406,100 @@ function initialize_tg!( end return end + +function initialize_tg!( + device_states, + static::PSY.StaticInjection, + dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, A, PSY.WPIDHY, P}}, + inner_vars::AbstractVector, +) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, P <: PSY.PSS} + + #Get mechanical torque to SyncMach + τm0 = inner_vars[τm_var] + τe0 = inner_vars[τe_var] + P0 = PSY.get_active_power(static) + + #Get parameters + tg = PSY.get_prime_mover(dynamic_device) + reg = PSY.get_reg(tg) + Kp = PSY.get_Kp(tg) + Ki = PSY.get_Ki(tg) + Kd = PSY.get_Kd(tg) + Ta = PSY.get_Ta(tg) + gate_openings = PSY.get_gate_openings(tg) + power_gate_openings = PSY.get_power_gate_openings(tg) + G_min, G_max = PSY.get_G_lim(tg) + P_min, P_max = PSY.get_P_lim(tg) + Tw = PSY.get_Tw(tg) + + #Compute controller parameters for equivalent TF + Kp_prime = (-Ta * Ki) + Kp + Kd_prime = (Ta^2 * Ki) - (Ta * Kp) + Kd + Ki_prime = Ki + + function f!(out, x) + P_ref = x[1] + x_g1 = x[2] + x_g2 = x[3] + x_g3 = x[4] + x_g4 = x[5] + x_g5 = x[6] + x_g6 = x[7] + x_g7 = x[8] + + pid_input = x_g1 + pi_out, dxg2_dt = pi_block(pid_input, x_g2, Kp_prime, Ki_prime) + pd_out, dxg4_dt = high_pass(pid_input, x_g4, Kd_prime, Ta) + + power_at_gate = + three_level_gate_to_power_map(x_g6, gate_openings, power_gate_openings) + + y_LL_out, dxg7_dt = lead_lag(power_at_gate, x_g7, 1.0, -Tw, Tw / 2.0) + + out[1] = y_LL_out - τm0 + out[2] = (P0 - P_ref) * reg - x_g1 + out[3] = dxg2_dt + out[4] = pi_out + pd_out - x_g3 + out[5] = dxg4_dt + out[6] = x_g3 - x_g5 + out[7] = x_g5 + out[8] = dxg7_dt + end + gate0 = three_level_power_to_gate_map(τm0, gate_openings, power_gate_openings) + P_ref_guess = P0 + xg1_guess = τe0 - P_ref_guess + + #x0 = [P_ref_guess, xg1_guess, 0.0, 0.0, 0.0, 0.0, gate0, 3.0 * τm0] + x0 = [P_ref_guess, 0.0, gate0, gate0, 0.0, gate0, gate0, 3.0 * τm0] + sol = NLsolve.nlsolve(f!, x0; ftol = STRICT_NLSOLVE_F_TOLERANCE) + if !NLsolve.converged(sol) + @warn("Initialization of Turbine Governor $(PSY.get_name(static)) failed") + else + sol_x0 = sol.zero + #Error if x_g3 is outside PI limits + if sol_x0[7] > G_max || sol_x0[7] < G_min + @error( + "WPIDHY Turbine Governor $(PSY.get_name(static)) $(sol_x0[7]) outside its gate limits $G_min, $G_max. Consider updating the operating point.", + ) + elseif τm0 > P_max || τm0 < P_min + @error( + "WPIDHY Turbine Governor $(PSY.get_name(static)) $(sol_x0[8]) outside its power limits $P_min, $P_max. Consider updating the operating point.", + ) + end + + #Update Control Refs + PSY.set_P_ref!(tg, sol_x0[1]) + set_P_ref(dynamic_device, sol_x0[1]) + #Update states + tg_ix = get_local_state_ix(dynamic_device, typeof(tg)) + tg_states = @view device_states[tg_ix] + tg_states[1] = sol_x0[2] + tg_states[2] = sol_x0[3] + tg_states[3] = sol_x0[4] + tg_states[4] = sol_x0[5] + tg_states[5] = sol_x0[6] + tg_states[6] = sol_x0[7] + tg_states[7] = sol_x0[8] + end + return +end diff --git a/src/models/common_controls.jl b/src/models/common_controls.jl index e8e55be5b..d420f3205 100644 --- a/src/models/common_controls.jl +++ b/src/models/common_controls.jl @@ -1606,7 +1606,7 @@ function three_level_gate_to_power_map( p0 = 0.0 g3 = 1.0 if input <= g0 - return zero(T) + return zero(Z) elseif g0 < input <= g1 return ((p1 - p0) / (g1 - g0)) * (input - g0) + p0 elseif g1 < input <= g2 diff --git a/src/models/generator_models/tg_models.jl b/src/models/generator_models/tg_models.jl index 445ff0f6a..2cf8baad4 100644 --- a/src/models/generator_models/tg_models.jl +++ b/src/models/generator_models/tg_models.jl @@ -43,6 +43,15 @@ function mass_matrix_tg_entries!( return end +function mass_matrix_tg_entries!( + mass_matrix, + tg::PSY.WPIDHY, + global_index::Base.ImmutableDict{Symbol, Int64}, +) + mass_matrix[global_index[:x_g1], global_index[:x_g1]] = PSY.get_T_reg(tg) + return +end + ################################## ##### Differential Equations ##### ################################## @@ -409,6 +418,102 @@ function mdl_tg_ode!( return end +function mdl_tg_ode!( + device_states::AbstractArray{<:ACCEPTED_REAL_TYPES}, + output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES}, + inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES}, + ω_sys::ACCEPTED_REAL_TYPES, + device::DynamicWrapper{PSY.DynamicGenerator{M, S, A, PSY.WPIDHY, P}}, + h, + t, +) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, P <: PSY.PSS} + + #Obtain references + P_ref = get_P_ref(device) + + #Obtain indices for component w/r to device + local_ix = get_local_state_ix(device, PSY.WPIDHY) + + #Define internal states for component + internal_states = @view device_states[local_ix] + x_g1 = internal_states[1] # Filter Input + x_g2 = internal_states[2] # PI Block + x_g3 = internal_states[3] # Regulator After PID Block + x_g4 = internal_states[4] # Derivative (High Pass) Block + x_g5 = internal_states[5] # Second Regulator Block + x_g6 = internal_states[6] # Gate State + x_g7 = internal_states[7] # Water Inertia State + + #Obtain external states inputs for component + external_ix = get_input_port_ix(device, PSY.WPIDHY) + ω = @view device_states[external_ix] + + # Read Inner Vars + τ_e = inner_vars[τe_var] + + #Get Parameters + tg = PSY.get_prime_mover(device) + T_reg = PSY.get_T_reg(tg) + reg = PSY.get_reg(tg) + Kp = PSY.get_Kp(tg) + Ki = PSY.get_Ki(tg) + Kd = PSY.get_Kd(tg)#Kd>0 + Ta = PSY.get_Ta(tg) + Tb = PSY.get_Tb(tg) + V_min, V_max = PSY.get_V_lim(tg) + G_min, G_max = PSY.get_G_lim(tg) + P_min, P_max = PSY.get_P_lim(tg) + Tw = PSY.get_Tw(tg) + D = PSY.get_D(tg) + gate_openings = PSY.get_gate_openings(tg) + power_gate_openings = PSY.get_power_gate_openings(tg) + + #Compute controller parameters for equivalent TF + Kp_prime = (-Ta * Ki) + Kp + Kd_prime = (Ta^2 * Ki) - (Ta * Kp) + Kd + Ki_prime = Ki + + x_in = τ_e - P_ref + #Compute block derivatives + _, dxg1_dt = low_pass_mass_matrix(x_in, x_g1, reg, T_reg) + pid_input = x_g1 - (ω[1] - ω_sys) + pi_out, dxg2_dt = pi_block(pid_input, x_g2, Kp_prime, Ki_prime) + pd_out, dxg4_dt = high_pass(pid_input, x_g4, Kd_prime, Ta) + pid_out = pi_out + pd_out + _, dxg3_dt = low_pass(pid_out, x_g3, 1.0, Ta) + _, dxg5_dt = low_pass(x_g3, x_g5, 1.0, Tb) + + #Set clamping for G_vel. + G_vel_sat = clamp(x_g5, V_min, V_max) + + # Compute integrator + xg6_sat, dxg6_dt = integrator_windup(G_vel_sat, x_g6, 1.0, 1.0, G_min, G_max) + + power_at_gate = + three_level_gate_to_power_map(xg6_sat, gate_openings, power_gate_openings) + + # Compute Lead-Lag Block + ll_out, dxg7_dt = lead_lag(power_at_gate, x_g7, 1.0, -Tw, Tw / 2.0) + Power_sat = clamp(ll_out, P_min, P_max) + + #Compute output torque + P_m = Power_sat - (D * (ω[1] - ω_sys)) + + #Compute 1 State TG ODE: + output_ode[local_ix[1]] = dxg1_dt + output_ode[local_ix[2]] = dxg2_dt + output_ode[local_ix[3]] = dxg3_dt + output_ode[local_ix[4]] = dxg4_dt + output_ode[local_ix[5]] = dxg5_dt + output_ode[local_ix[6]] = dxg6_dt + output_ode[local_ix[7]] = dxg7_dt + + #Update mechanical torque + inner_vars[τm_var] = P_m / ω[1] + + return +end + function mdl_tg_ode!( device_states::AbstractArray{<:ACCEPTED_REAL_TYPES}, output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES}, diff --git a/src/post_processing/post_proc_generator.jl b/src/post_processing/post_proc_generator.jl index db8d715fb..d3da8ca96 100644 --- a/src/post_processing/post_proc_generator.jl +++ b/src/post_processing/post_proc_generator.jl @@ -1179,3 +1179,34 @@ function _mechanical_torque( τm = Pe ./ ω return ts, τm end + +function _mechanical_torque( + tg::PSY.WPIDHY, + name::String, + res::SimulationResults, + dt::Union{Nothing, Float64, Vector{Float64}}, +) + # Get params + D = PSY.get_D(tg) + gate_openings = PSY.get_gate_openings(tg) + power_gate_openings = PSY.get_power_gate_openings(tg) + Tw = PSY.get_Tw(tg) + setpoints = get_setpoints(res) + ω_ref = setpoints[name]["ω_ref"] + + # Get state results + ts, x_g7 = post_proc_state_series(res, (name, :x_g7), dt) + _, x_g6 = post_proc_state_series(res, (name, :x_g6), dt) + _, ω = post_proc_state_series(res, (name, :ω), dt) + Pm = similar(x_g7) + + for (ix, x7) in enumerate(x_g7) + x6 = x_g6[ix] + power_at_gate = + three_level_gate_to_power_map(x6, gate_openings, power_gate_openings) + ll_out, _ = lead_lag(power_at_gate, x7, 1.0, -Tw, Tw / 2.0) + Pm[ix] = ll_out - D * (ω[ix] - ω_ref) + end + τm = Pm ./ ω + return ts, τm +end diff --git a/test/benchmarks/psse/WPIDHY/ThreeBusMulti.raw b/test/benchmarks/psse/WPIDHY/ThreeBusMulti.raw new file mode 100644 index 000000000..6c04ff717 --- /dev/null +++ b/test/benchmarks/psse/WPIDHY/ThreeBusMulti.raw @@ -0,0 +1,32 @@ +0, 100.00, 33, 0, 0, 60.00 / PSS(R)E 33 RAW created by rawd33 TUE, JUL 21 2020 17:55 + + + 101,'BUS 1', 138.0000,3, 1, 1, 1,1.05000, 0.0000,1.10000,0.90000,1.10000,0.90000 + 102,'BUS 2', 138.0000,2, 1, 1, 1,1.02000, -0.9440,1.10000,0.90000,1.10000,0.90000 + 103,'BUS 3', 138.0000,1, 1, 1, 1,0.99341, -8.7697,1.10000,0.90000,1.10000,0.90000 + 0 /End of Bus data, Begin Load data + 103,'1 ',1, 1, 1, 200.000, 30.000, 0.000, 0.000, 0.000, 0.000, 1,1,0 + 0 /End of Load data, Begin Fixed shunt data + 0 /End of Fixed shunt data, Begin Generator data + 101,'1 ', 133.335, 73.271, 100.000, -100.000,1.05000, 0, 100.000, 0.00000E+0, 1.00000E-5, 0.00000E+0, 0.00000E+0,1.00000,1, 100.0, 318.000, 0.000, 1,1.0000 + 102,'1 ', 70.000, -3.247, 100.000, -100.000,1.02000, 0, 100.000, 0.00000E+0, 2.500E-1, 0.00000E+0, 0.00000E+0,1.00000,1, 100.0, 318.000, 0.000, 1,1.0000 + 0 /End of Generator data, Begin Branch data + 101, 102,'1 ', 1.00000E-2, 1.20000E-1, 0.00000, 250.00, 250.00, 250.00, 0.00000, 0.00000, 0.00000, 0.00000,1,1, 0.00, 1,1.0000 + 101, 103,'1 ', 1.00000E-2, 1.20000E-1, 0.00000, 250.00, 250.00, 250.00, 0.00000, 0.00000, 0.00000, 0.00000,1,1, 0.00, 1,1.0000 + 102, 103,'1 ', 1.00000E-2, 1.20000E-1, 0.00000, 250.00, 250.00, 250.00, 0.00000, 0.00000, 0.00000, 0.00000,1,1, 0.00, 1,1.0000 + 0 /End of Branch data, Begin Transformer data + 0 /End of Transformer data, Begin Area interchange data + 0 /End of Area interchange data, Begin Two-terminal dc line data + 0 /End of Two-terminal dc line data, Begin VSC dc line data + 0 /End of VSC dc line data, Begin Impedance correction table data + 0 /End of Impedance correction table data, Begin Multi-terminal dc line data + 0 /End of Multi-terminal dc line data, Begin Multi-section line data + 0 /End of Multi-section line data, Begin Zone data + 0 /End of Zone data, Begin Inter-area transfer data + 0 /End of Inter-area transfer data, Begin Owner data + 0 /End of Owner data, Begin FACTS device data + 0 /End of FACTS device data, Begin Switched shunt data + 0 /End of Switched shunt data, Begin GNE device data + 0 /End of GNE device data, Begin Induction machine data + 0 /End of Induction machine data +Q \ No newline at end of file diff --git a/test/benchmarks/psse/WPIDHY/ThreeBus_WPIDHY.dyr b/test/benchmarks/psse/WPIDHY/ThreeBus_WPIDHY.dyr new file mode 100644 index 000000000..d71618c71 --- /dev/null +++ b/test/benchmarks/psse/WPIDHY/ThreeBus_WPIDHY.dyr @@ -0,0 +1,10 @@ + 101 'GENCLS' 1 0.0000 0.0000 / + 102 'GENROU' 1 8.0000 0.30000E-01 0.40000 0.50000E-01 + 6.1750 0.50000E-01 1.8000 1.7000 0.30000 + 0.55000 0.25000 0.20000 0.0000 1.0000 / + 102 'SEXS' 1 0.4 5.0 20.0 1.0 -50.0 50.0 / + 102 'WPIDHY' 1 0.10000E-01 -0.50000E-01 3.0000 0.50000 + 3.0000 0.20000E-01 0.10000E-01 0.25000E-01 -0.25000E-01 + 1.0000 0.20000 3.2500 52001. 0.0000 + 0.0000 0.12000 0.54500 0.40000 0.82700 + 0.83000 1.0000 / diff --git a/test/results/results_eigenvalues.jl b/test/results/results_eigenvalues.jl index c7cb82601..2d3c2df98 100644 --- a/test/results/results_eigenvalues.jl +++ b/test/results/results_eigenvalues.jl @@ -983,3 +983,21 @@ test59_eigvals = [ -0.36823401395037775 - 0.5649228138733879im -0.36823401395037775 + 0.5649228138733879im ] + +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 cb38572b8..d62797126 100644 --- a/test/results/results_initial_conditions.jl +++ b/test/results/results_initial_conditions.jl @@ -1806,3 +1806,33 @@ test59_x0_init = Dict( 2.153115531256307 ], ) + +test60_x0_init = Dict( + "V_R" => [ + 1.05 + 1.0197944718502572 + 0.9923907751848658 + ], + "V_I" => [ + 0.0 + -0.020475233421259967 + -0.1243212484458825 + ], + "generator-102-1" => [ + 0.7538967192769663 + 0.5555487379562241 + 0.7042452148052648 + 0.7246287886385532 + 0.9158416020737494 + 1.0 + 1.4986692863524897 + 0.044960078577949814 + 0.0 + 3.871458709170383e-12 + 3.871569731472846e-12 + 0.0 + 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