From 3f959bc3dd605a83be3677a216b72dbedec6fa5e Mon Sep 17 00:00:00 2001 From: rodrigomha Date: Wed, 21 Aug 2024 16:39:13 -0700 Subject: [PATCH 1/9] add gate to power function for PIDGOV --- src/models/common_controls.jl | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/src/models/common_controls.jl b/src/models/common_controls.jl index c320f843b..df26482d7 100644 --- a/src/models/common_controls.jl +++ b/src/models/common_controls.jl @@ -1591,3 +1591,27 @@ function integrator_nonwindup( y_sat, dydt_scaled = integrator_nonwindup_mass_matrix(u, y, K, T, y_min, y_max) return y_sat, (1.0 / T) * dydt_scaled end + +""" +Three level gate to power transformation for PIDGOV +(G0, 0), (G1, P1), (G2, P2), (1, P3) are (x,y) coordinates of Power vs Gate functions +""" +function three_level_gate_to_power_map( + input::Z, + gates::Tuple{Float64, Float64, Float64}, + powers::Tuple{Float64, Float64, Float64}, +) where {Z <: ACCEPTED_REAL_TYPES} + g0, g1, g2 = gates + p1, p2, p3 = powers + p0 = 0.0 + g3 = 1.0 + if input <= g0 + return zero(T) + elseif g0 < input <= g1 + return ((p1 - p0) / (g1 - g0)) * (input - g0) + p0 + elseif g1 < input <= g2 + return ((p2 - p1) / (g2 - g1)) * (input - g1) + p1 + elseif g2 < input + return ((p3 - p2) / (g3 - g2)) * (input - g2) + p2 + end +end From b5f4f39e08d122a76bf1bba5f197bbc14bfd8cd8 Mon Sep 17 00:00:00 2001 From: rodrigomha Date: Wed, 21 Aug 2024 16:39:28 -0700 Subject: [PATCH 2/9] update inner vars to update electric torque --- src/models/device.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/models/device.jl b/src/models/device.jl index 5e5a1c078..4d1fd9549 100644 --- a/src/models/device.jl +++ b/src/models/device.jl @@ -333,10 +333,14 @@ function _update_inner_vars!( I_d = (1.0 / (R^2 + Xq_pp * Xd_pp)) * (-R * (V_dq[1] - ψq_pp) + Xq_pp * (-V_dq[2] + ψd_pp)) + I_q = + (1.0 / (R^2 + Xq_pp * Xd_pp)) * (Xd_pp * (V_dq[1] - ψq_pp) + R * (-V_dq[2] + ψd_pp)) Se = saturation_function(machine, ψ_pp) Xad_Ifd = eq_p + (Xd - Xd_p) * (γ_d1 * I_d - γ_d2 * ψ_kd + γ_d2 * eq_p) + Se * ψd_pp + τ_e = I_d * (V_dq[1] + I_d * R) + I_q * (V_dq[2] + I_q * R) #Update Xad_Ifd + inner_vars[τe_var] = τ_e inner_vars[Xad_Ifd_var] = Xad_Ifd return end @@ -386,11 +390,14 @@ function _update_inner_vars!( #Currents I_d = (1.0 / (R^2 + Xd_pp^2)) * (-R * (V_d + ψq_pp) + Xd_pp * (ψd_pp - V_q)) + I_q = (1.0 / (R^2 + Xd_pp^2)) * (Xd_pp * (V_d + ψq_pp) + R * (ψd_pp - V_q)) Se = saturation_function(machine, eq_p) Xad_Ifd = eq_p + Se * eq_p + (Xd - Xd_p) * (I_d + γ_d2 * (eq_p - ψ_kd - (Xd_p - Xl) * I_d)) + τ_e = I_d * (V_d + I_d * R) + I_q * (V_q + I_q * R) #Update Xad_Ifd + inner_vars[τe_var] = τ_e inner_vars[Xad_Ifd_var] = Xad_Ifd return end @@ -442,11 +449,14 @@ function _update_inner_vars!( #Currents I_d = (1.0 / (R^2 + Xd_pp^2)) * (-R * (V_d - ψq_pp) + Xq_pp * (-V_q + ψd_pp)) + I_q = (1.0 / (R^2 + Xd_pp^2)) * (Xd_pp * (V_d - ψq_pp) + R * (-V_q + ψd_pp)) Se = saturation_function(machine, ψ_pp) Xad_Ifd = eq_p + Se * ψd_pp + (Xd - Xd_p) * (I_d + γ_d2 * (eq_p - ψ_kd - (Xd_p - Xl) * I_d)) + τ_e = I_d * (V_d + I_d * R) + I_q * (V_q + I_q * R) #Update Xad_Ifd + inner_vars[τe_var] = τ_e inner_vars[Xad_Ifd_var] = Xad_Ifd return end From 1997c56d89ece85dd74538499e9324a6ae9e6dd1 Mon Sep 17 00:00:00 2001 From: rodrigomha Date: Wed, 21 Aug 2024 16:39:47 -0700 Subject: [PATCH 3/9] add PIDGOV model --- src/models/generator_models/tg_models.jl | 117 +++++++++++++++++++++++ 1 file changed, 117 insertions(+) diff --git a/src/models/generator_models/tg_models.jl b/src/models/generator_models/tg_models.jl index 67755c0f6..d473c63b6 100644 --- a/src/models/generator_models/tg_models.jl +++ b/src/models/generator_models/tg_models.jl @@ -1,3 +1,7 @@ +################################## +###### Mass Matrix Entries ####### +################################## + function mass_matrix_tg_entries!( mass_matrix, tg::TG, @@ -27,6 +31,22 @@ function mass_matrix_tg_entries!( return end +function mass_matrix_tg_entries!( + mass_matrix, + tg::PSY.PIDGOV, + global_index::Base.ImmutableDict{Symbol, Int64}, +) + mass_matrix[global_index[:x_g1], global_index[:x_g1]] = PSY.get_T_reg(tg) + mass_matrix[global_index[:x_g3], global_index[:x_g3]] = PSY.get_Ta(tg) + mass_matrix[global_index[:x_g4], global_index[:x_g4]] = PSY.get_Ta(tg) + mass_matrix[global_index[:x_g5], global_index[:x_g5]] = PSY.get_Ta(tg) + return +end + +################################## +##### Differential Equations ##### +################################## + function mdl_tg_ode!( device_states::AbstractArray{<:ACCEPTED_REAL_TYPES}, ::AbstractArray{<:ACCEPTED_REAL_TYPES}, @@ -388,3 +408,100 @@ 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.PIDGOV, 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.PIDGOV) + + #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 PI 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.PIDGOV) + ω = @view device_states[external_ix] + + # Read Inner Vars + τ_e = inner_vars[τe_var] + + #Get Parameters + tg = PSY.get_prime_mover(device) + feedback_flag = PSY.get_feedback_flag(tg) + Rperm = PSY.get_Rperm(tg) + T_reg = PSY.get_T_reg(tg) + Kp = PSY.get_Kp(tg) + Ki = PSY.get_Ki(tg) + Kd = PSY.get_Kd(tg) + Ta = PSY.get_Ta(tg) + Tb = PSY.get_Tb(tg) + D_turb = PSY.get_D_turb(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) + A_tw = PSY.get_A_tw(tg) + Tw = PSY.get_Tw(tg) + V_min, V_max = PSY.get_V_lim(tg) + + #Compute auxiliary parameters + if feedback_flag == 0 + x_in = P_ref - τ_e + else + x_in = P_ref - x_g6 + end + + #Compute block derivatives + _, dxg1_dt = low_pass_mass_matrix(x_in, x_g1, Rperm, T_reg) + pid_input = x_g1 - ω[1] + pi_out, dxg2_dt = pi_block(pid_input, x_g2, Kp, Ki) + _, dxg3_dt = low_pass_mass_matrix(pi_out, x_g3, 1.0, Ta) + pd_out, dxg4_dt = high_pass_mass_matrix(pid_input, x_g4, Kd, Ta) + _, dxg5_dt = low_pass_mass_matrix(x_g3 + pd_out, x_g5, 1.0, Ta) + + # Compute integrator + integrator_input = (1.0 / Tb) * (x_g5 - x_g6) + integrator_input_sat = clamp(integrator_input, V_min, V_max) + xg6_sat, dxg6_dt = + integrator_nonwindup(integrator_input_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 + Tz = A_tw * Tw + ll_out, dxg7_dt = lead_lag(power_at_gate, x_g7, 1.0, -Tz, Tz / 2.0) + + #Compute output torque + P_m = ll_out - D_turb * (ω[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 From e99fd56da7567f5f1bd1f681ab0ba9354ba52d09 Mon Sep 17 00:00:00 2001 From: rodrigomha Date: Wed, 21 Aug 2024 16:40:03 -0700 Subject: [PATCH 4/9] [WIP] add test PIDGOV --- test/benchmarks/psse/PIDGOV/ThreeBusMulti.raw | 32 +++++++++++++++++++ .../psse/PIDGOV/ThreeBus_PIDGOV_GateFlag.dyr | 10 ++++++ .../psse/PIDGOV/ThreeBus_PIDGOV_PowerFlag.dyr | 10 ++++++ test/test_case58_pidgov.jl | 12 +++++++ 4 files changed, 64 insertions(+) create mode 100644 test/benchmarks/psse/PIDGOV/ThreeBusMulti.raw create mode 100644 test/benchmarks/psse/PIDGOV/ThreeBus_PIDGOV_GateFlag.dyr create mode 100644 test/benchmarks/psse/PIDGOV/ThreeBus_PIDGOV_PowerFlag.dyr create mode 100644 test/test_case58_pidgov.jl diff --git a/test/benchmarks/psse/PIDGOV/ThreeBusMulti.raw b/test/benchmarks/psse/PIDGOV/ThreeBusMulti.raw new file mode 100644 index 000000000..7ab73f141 --- /dev/null +++ b/test/benchmarks/psse/PIDGOV/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, 250.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 ', 153.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 ', 100.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/PIDGOV/ThreeBus_PIDGOV_GateFlag.dyr b/test/benchmarks/psse/PIDGOV/ThreeBus_PIDGOV_GateFlag.dyr new file mode 100644 index 000000000..bb2ec7b09 --- /dev/null +++ b/test/benchmarks/psse/PIDGOV/ThreeBus_PIDGOV_GateFlag.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 'PIDGOV' 1 1 0.50000E-01 0.10000E-01 5.0000 + 1.0000 0.10000 0.50000 0.99000 0.10000E-02 + 0.10000 0.70000 0.30000 0.90000 0.50000 + 0.90000 1.0000 0.50000E-01 1.0000 0.60000 + 0.40000E-01 -0.40000E-01 / diff --git a/test/benchmarks/psse/PIDGOV/ThreeBus_PIDGOV_PowerFlag.dyr b/test/benchmarks/psse/PIDGOV/ThreeBus_PIDGOV_PowerFlag.dyr new file mode 100644 index 000000000..aaa0597c7 --- /dev/null +++ b/test/benchmarks/psse/PIDGOV/ThreeBus_PIDGOV_PowerFlag.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 'PIDGOV' 1 0 0.50000E-01 0.10000E-01 5.0000 + 1.0000 0.10000 0.50000 0.99000 0.10000E-02 + 0.10000 0.70000 0.30000 0.90000 0.50000 + 0.90000 1.0000 0.50000E-01 1.0000 0.60000 + 0.40000E-01 -0.40000E-01 / diff --git a/test/test_case58_pidgov.jl b/test/test_case58_pidgov.jl new file mode 100644 index 000000000..22bbb88d7 --- /dev/null +++ b/test/test_case58_pidgov.jl @@ -0,0 +1,12 @@ +""" +Validation PSSE/PIDGOV: +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 #################### +################################################## + +raw_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/PIDGOV/ThreeBusMulti.raw") +dyr_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/PIDGOV/ThreeBus_PIDGOV_PowerFlag.dyr") From 06bc65de04a69e4e271c63055ae1deb9b9e6e05e Mon Sep 17 00:00:00 2001 From: rodrigomha Date: Fri, 23 Aug 2024 16:06:10 -0700 Subject: [PATCH 5/9] add initialization routine for pidgov --- .../generator_components/init_tg.jl | 95 +++++++++++++++++++ 1 file changed, 95 insertions(+) diff --git a/src/initialization/generator_components/init_tg.jl b/src/initialization/generator_components/init_tg.jl index 7b91d492a..178834af5 100644 --- a/src/initialization/generator_components/init_tg.jl +++ b/src/initialization/generator_components/init_tg.jl @@ -311,3 +311,98 @@ function initialize_tg!( end return end + +function initialize_tg!( + device_states, + static::PSY.StaticInjection, + dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, A, PSY.PIDGOV, 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] + P0 = PSY.get_active_power(static) + Δω = 0.0 + #Get parameters + tg = PSY.get_prime_mover(dynamic_device) + feedback_flag = PSY.get_feedback_flag(tg) + Rperm = PSY.get_Rperm(tg) + Kp = PSY.get_Kp(tg) + Ki = PSY.get_Ki(tg) + Kd = PSY.get_Kd(tg) + Ta = PSY.get_Ta(tg) + Tb = PSY.get_Tb(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) + A_tw = PSY.get_A_tw(tg) + Tw = PSY.get_Tw(tg) + + 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] + + if feedback_flag == 0 + P_in = P_ref - P0 + else + P_in = P_ref - x_g6 + end + pid_input = x_g1 + pi_out, dxg2_dt = pi_block(pid_input, x_g2, Kp, Ki) + pd_out, dxg4_dt = high_pass_mass_matrix(pid_input, x_g4, Kd, Ta) + + integrator_input = (1.0 / Tb) * (x_g5 - x_g6) + power_at_gate = + three_level_gate_to_power_map(x_g6, gate_openings, power_gate_openings) + Tz = A_tw * Tw + y_LL_out, dxg7_dt = lead_lag(power_at_gate, x_g7, 1.0, -Tz, Tz / 2.0) + + out[1] = y_LL_out - τm0 + out[2] = (Rperm * P_in - x_g1) + out[3] = dxg2_dt + out[4] = pi_out - x_g3 + out[5] = dxg4_dt + out[6] = x_g3 + pd_out - x_g5 + out[7] = integrator_input + out[8] = dxg7_dt + end + gate0 = three_level_power_to_gate_map(τm0, gate_openings, power_gate_openings) + if feedback_flag == 0 + P_ref_guess = P0 + else + P_ref_guess = gate0 + end + 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( + "PIDGOV Turbine Governor $(PSY.get_name(static)) $(sol_x0[4]) outside its gate limits $G_min, $G_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 From e21638a440f17df580df2221c456ab4ddfd0b7a0 Mon Sep 17 00:00:00 2001 From: rodrigomha Date: Fri, 23 Aug 2024 16:06:23 -0700 Subject: [PATCH 6/9] add inverse map for gate to power --- src/models/common_controls.jl | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/src/models/common_controls.jl b/src/models/common_controls.jl index df26482d7..e8e55be5b 100644 --- a/src/models/common_controls.jl +++ b/src/models/common_controls.jl @@ -1615,3 +1615,27 @@ function three_level_gate_to_power_map( return ((p3 - p2) / (g3 - g2)) * (input - g2) + p2 end end + +""" +Inverse Three level gate to power transformation for PIDGOV +(G0, 0), (G1, P1), (G2, P2), (1, P3) are (x,y) coordinates of Power vs Gate functions +""" +function three_level_power_to_gate_map( + power_input::Z, + gates::Tuple{Float64, Float64, Float64}, + powers::Tuple{Float64, Float64, Float64}, +) where {Z <: ACCEPTED_REAL_TYPES} + g0, g1, g2 = gates + p1, p2, p3 = powers + p0 = 0.0 + g3 = 1.0 + if power_input <= p0 + return g0 + elseif p0 < power_input <= p1 + return ((g1 - g0) / (p1 - p0)) * (power_input - p0) + g0 + elseif p1 < power_input <= p2 + return ((g2 - g1) / (p2 - p1)) * (power_input - p1) + g1 + elseif p2 < power_input + return ((g3 - g2) / (p3 - p2)) * (power_input - p2) + g2 + end +end From 081299985995ad76b126656a99763671e16a3af4 Mon Sep 17 00:00:00 2001 From: rodrigomha Date: Fri, 23 Aug 2024 16:06:36 -0700 Subject: [PATCH 7/9] fix error for frequency difference --- src/models/generator_models/tg_models.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/models/generator_models/tg_models.jl b/src/models/generator_models/tg_models.jl index d473c63b6..445ff0f6a 100644 --- a/src/models/generator_models/tg_models.jl +++ b/src/models/generator_models/tg_models.jl @@ -469,7 +469,7 @@ function mdl_tg_ode!( #Compute block derivatives _, dxg1_dt = low_pass_mass_matrix(x_in, x_g1, Rperm, T_reg) - pid_input = x_g1 - ω[1] + pid_input = x_g1 - (ω[1] - ω_sys) pi_out, dxg2_dt = pi_block(pid_input, x_g2, Kp, Ki) _, dxg3_dt = low_pass_mass_matrix(pi_out, x_g3, 1.0, Ta) pd_out, dxg4_dt = high_pass_mass_matrix(pid_input, x_g4, Kd, Ta) From a5a4e1f1bea53889f7e418ddfd142de6589985b5 Mon Sep 17 00:00:00 2001 From: rodrigomha Date: Fri, 23 Aug 2024 16:06:48 -0700 Subject: [PATCH 8/9] add post calculation for mechanical torque for pidgov --- src/post_processing/post_proc_generator.jl | 35 ++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/src/post_processing/post_proc_generator.jl b/src/post_processing/post_proc_generator.jl index 1f25258dd..e9a63b74f 100644 --- a/src/post_processing/post_proc_generator.jl +++ b/src/post_processing/post_proc_generator.jl @@ -1081,3 +1081,38 @@ function _mechanical_torque( τm = ((x_g4 .- q_nl) .* h * At - D_T * Δω .* x_g3) ./ ω return ts, τm end + +""" +Function to obtain the mechanical torque time series of a Dynamic Generator with HydroTurbineGov (HYGOV) Turbine Governor. + +""" +function _mechanical_torque( + tg::PSY.PIDGOV, + name::String, + res::SimulationResults, + dt::Union{Nothing, Float64, Vector{Float64}}, +) + # Get params + D_turb = PSY.get_D_turb(tg) + gate_openings = PSY.get_gate_openings(tg) + power_gate_openings = PSY.get_power_gate_openings(tg) + A_tw = PSY.get_A_tw(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) + Pe = 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) + Tz = A_tw * Tw + ll_out, _ = lead_lag(power_at_gate, x7, 1.0, -Tz, Tz / 2.0) + Pe[ix] = ll_out - D_turb * (ω[ix] - ω_ref) + end + τm = Pe ./ ω + return ts, τm +end From 077891743b3b7b2af8cf13372d2492c12734eb78 Mon Sep 17 00:00:00 2001 From: rodrigomha Date: Fri, 23 Aug 2024 16:06:59 -0700 Subject: [PATCH 9/9] add test58 --- test/benchmarks/psse/PIDGOV/ThreeBusMulti.raw | 6 +- test/results/results_eigenvalues.jl | 36 ++++ test/results/results_initial_conditions.jl | 60 ++++++ test/test_case58_pidgov.jl | 202 +++++++++++++++++- 4 files changed, 299 insertions(+), 5 deletions(-) diff --git a/test/benchmarks/psse/PIDGOV/ThreeBusMulti.raw b/test/benchmarks/psse/PIDGOV/ThreeBusMulti.raw index 7ab73f141..318f19291 100644 --- a/test/benchmarks/psse/PIDGOV/ThreeBusMulti.raw +++ b/test/benchmarks/psse/PIDGOV/ThreeBusMulti.raw @@ -5,11 +5,11 @@ 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, 250.000, 30.000, 0.000, 0.000, 0.000, 0.000, 1,1,0 + 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 ', 153.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 ', 100.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 + 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 diff --git a/test/results/results_eigenvalues.jl b/test/results/results_eigenvalues.jl index e8dd1a497..b37b79b4c 100644 --- a/test/results/results_eigenvalues.jl +++ b/test/results/results_eigenvalues.jl @@ -935,3 +935,39 @@ test_57_eigvals = [ 0.0 + 0.0im, 0.0 + 0.0im, ] + +test58_eigvals_PowerFlag = [ + -99.99999037034084 + 0.0im + -41.17481961789709 + 0.0im + -38.789855128930164 + 0.0im + -5.647427439841257 + 0.0im + -3.8982943228132214 - 1.5754284098308229im + -3.8982943228132214 + 1.5754284098308229im + -2.000000000000001 + 0.0im + -0.9492553729947647 + 0.0im + -0.8686675700247483 - 8.496111510550772im + -0.8686675700247483 + 8.496111510550772im + -0.28282148961179976 - 1.1097759323927108im + -0.28282148961179976 + 1.1097759323927108im + -0.1975067616652887 - 0.17956822187766im + -0.1975067616652887 + 0.17956822187766im + -0.11697960134622151 + 0.0im +] + +test58_eigvals_GateFlag = [ + -100.00010646060846 + 0.0im + -41.17482850153586 + 0.0im + -38.789870976223455 + 0.0im + -5.718238174509167 + 0.0im + -3.7935658880265133 + 0.0im + -2.0454156913133885 + 0.0im + -1.9999999999999976 + 0.0im + -1.2228555790477964 - 0.7200047656459937im + -1.2228555790477964 + 0.7200047656459937im + -0.9655206559688869 + 0.0im + -0.9001334037955144 - 8.447836359367825im + -0.9001334037955144 + 8.447836359367825im + -0.19849913626541832 - 0.17939599481963514im + -0.19849913626541832 + 0.17939599481963514im + -0.04238523317744444 + 0.0im +] diff --git a/test/results/results_initial_conditions.jl b/test/results/results_initial_conditions.jl index df2b268d3..d54376ab5 100644 --- a/test/results/results_initial_conditions.jl +++ b/test/results/results_initial_conditions.jl @@ -1721,3 +1721,63 @@ test_57_x0_init = Dict( 1.0 ], ) + +test58_x0_init_PowerFlag = 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.7042452148052647 + 0.7246287886385533 + 0.9158416020737494 + 1.0 + 1.4986692863524897 + 0.044960078577949814 + 0.0 + 0.95 + 0.95 + 0.0 + 0.95 + 0.95 + 2.099999999999999 + ], +) + +test58_x0_init_GateFlag = 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.7042452148052647 + 0.7246287886385533 + 0.9158416020737494 + 1.0 + 1.4986692863524897 + 0.044960078577949814 + 0.0 + 0.95 + 0.95 + 0.0 + 0.95 + 0.95 + 2.099999999999999 + ], +) diff --git a/test/test_case58_pidgov.jl b/test/test_case58_pidgov.jl index 22bbb88d7..f75823e25 100644 --- a/test/test_case58_pidgov.jl +++ b/test/test_case58_pidgov.jl @@ -8,5 +8,203 @@ The fault drop the line connecting the infinite bus and GENROU ############### SOLVE PROBLEM #################### ################################################## -raw_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/PIDGOV/ThreeBusMulti.raw") -dyr_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/PIDGOV/ThreeBus_PIDGOV_PowerFlag.dyr") +# Define dyr files + +names = ["PIDGOV Power Flag", "PIDGOV Gate Flag"] + +dyr_files = [ + joinpath(TEST_FILES_DIR, "benchmarks/psse/PIDGOV/ThreeBus_PIDGOV_PowerFlag.dyr"), + joinpath(TEST_FILES_DIR, "benchmarks/psse/PIDGOV/ThreeBus_PIDGOV_GateFlag.dyr"), +] + +csv_files = [ + joinpath(TEST_FILES_DIR, "benchmarks/psse/PIDGOV/PIDGOV_results_PowerFlag.csv"), + joinpath(TEST_FILES_DIR, "benchmarks/psse/PIDGOV/PIDGOV_results_GateFlag.csv"), +] + +init_conditions = [test58_x0_init_PowerFlag, test58_x0_init_GateFlag] + +eigs_values = [test58_eigvals_PowerFlag, test58_eigvals_GateFlag] + +raw_file_dir = joinpath(TEST_FILES_DIR, "benchmarks/psse/PIDGOV/ThreeBusMulti.raw") +tspan = (0.0, 20.0) + +function test_pidgov_implicit(dyr_file, csv_file, init_cond, eigs_value) + path = (joinpath(pwd(), "test-psse-pidgov")) + !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_pidgov_mass_matrix(dyr_file, csv_file, init_cond, eigs_value) + path = (joinpath(pwd(), "test-psse-pidgov")) + !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 58 PIDGOV 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_pidgov_implicit(dyr_file, csv_file, init_cond, eigs_value) + end + end +end + +@testset "Test 58 PIDGOV 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_pidgov_mass_matrix(dyr_file, csv_file, init_cond, eigs_value) + end + end +end