From 0c375d2d25901a3347727c3e505c30d54c252c48 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebasti=C3=A1n=20M?= <104728742+SebastianManriqueM@users.noreply.github.com> Date: Thu, 5 Sep 2024 15:05:28 -0600 Subject: [PATCH] Add AVR ST6B (#388) * add mass matrix entries for ST6B * Add Comment * Add ST6B model * Add initialization for ST6B * add mass matrix entries for ST6B * Add Test 59 * Fix Typo --------- Co-authored-by: Machado --- .../generator_components/init_avr.jl | 81 +++++++++ src/models/generator_models/avr_models.jl | 90 ++++++++++ src/post_processing/post_proc_generator.jl | 63 +++++++ test/benchmarks/psse/ST6B/ThreeBusMulti.raw | 32 ++++ test/benchmarks/psse/ST6B/ThreeBus_ST6B.dyr | 10 ++ test/results/results_eigenvalues.jl | 12 ++ test/results/results_initial_conditions.jl | 25 +++ test/test_case59_st6b.jl | 157 ++++++++++++++++++ 8 files changed, 470 insertions(+) create mode 100644 test/benchmarks/psse/ST6B/ThreeBusMulti.raw create mode 100644 test/benchmarks/psse/ST6B/ThreeBus_ST6B.dyr create mode 100644 test/test_case59_st6b.jl diff --git a/src/initialization/generator_components/init_avr.jl b/src/initialization/generator_components/init_avr.jl index 4d18c35a5..87965eadd 100644 --- a/src/initialization/generator_components/init_avr.jl +++ b/src/initialization/generator_components/init_avr.jl @@ -571,3 +571,84 @@ function initialize_avr!( avr_states[5] = -Kf / Tf * Vf0 return end + +function initialize_avr!( + device_states, + static::PSY.StaticInjection, + dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, PSY.ST6B, TG, P}}, + inner_vars::AbstractVector, +) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS} + #Obtain Vf0 solved from Machine + Vf0 = inner_vars[Vf_var] + #Obtain measured terminal voltage + Vt0 = sqrt(inner_vars[VR_gen_var]^2 + inner_vars[VI_gen_var]^2) + #Obtain field winding current + Ifd0 = inner_vars[Xad_Ifd_var] + + #Get parameters + avr = PSY.get_avr(dynamic_device) + Tr = PSY.get_Tr(avr) + K_pa = PSY.get_K_pa(avr) #k_pa>0 + K_ia = PSY.get_K_ia(avr) + K_da = PSY.get_K_da(avr) + T_da = PSY.get_T_da(avr) + Va_min, Va_max = PSY.get_Va_lim(avr) + K_ff = PSY.get_K_ff(avr) + K_m = PSY.get_K_m(avr) + K_ci = PSY.get_K_ci(avr) #K_cl in pss + K_lr = PSY.get_K_lr(avr) + I_lr = PSY.get_I_lr(avr) + Vr_min, Vr_max = PSY.get_Vr_lim(avr) + Kg = PSY.get_Kg(avr) + Tg = PSY.get_Tg(avr) #T_g>0 + + #Compute auxiliary parameters + Vr = Vf0 / Vt0 + Vg = Kg * Vf0 + + Vr2 = max(Vr_min, ((I_lr * K_ci) - Ifd0) * K_lr) + if (Vr - Vr2 > eps()) # or using STRICT_NLSOLVE_F_TOLERANCE? + + #Check saturation + + Va = (Vr + K_m * Vg) / (K_ff + K_m) + + #Check ss error according to control parameters + if K_ia < eps() + Vref0 = Va / K_pa + Vt0 + error( + "AVR in $(PSY.get_name(dynamic_device)) has non positive integrator gain. Please update it.", + ) + else + Vref0 = Vt0 + end + + if !((Vr < Vr_max) && (Vr > Vr_min)) + @error( + "AVR controller in $(PSY.get_name(dynamic_device)) is saturated. Consider updating the operating point." + ) + end + else + error( + "Current limiter of AVR in $(PSY.get_name(dynamic_device)) is activated. Consider updating the operating point.", + ) + end + + PSY.set_V_ref!(avr, Vref0) + set_V_ref(dynamic_device, Vref0) + @warn( + "I_LR parameter was updated from $(I_lr) to $(Ifd0) of $(PSY.get_name(dynamic_device)) to have zero current field error." + ) + PSY.set_I_lr!(avr, Ifd0 * 10.0) + + #States of ST6B [:Vm, :x_i, :x_d, :Vg] + #Update AVR states + avr_ix = get_local_state_ix(dynamic_device, PSY.ST6B) + avr_states = @view device_states[avr_ix] + avr_states[1] = Vt0 + avr_states[2] = Va / K_ia + avr_states[3] = 0 + avr_states[4] = Vg + + return +end diff --git a/src/models/generator_models/avr_models.jl b/src/models/generator_models/avr_models.jl index 220fc2e4a..748e38362 100644 --- a/src/models/generator_models/avr_models.jl +++ b/src/models/generator_models/avr_models.jl @@ -1,3 +1,7 @@ +################################## +###### Mass Matrix Entries ####### +################################## + function mass_matrix_avr_entries!( mass_matrix, avr::T, @@ -61,6 +65,20 @@ function mass_matrix_avr_entries!( return end +function mass_matrix_avr_entries!( + mass_matrix, + avr::PSY.ST6B, + global_index::Base.ImmutableDict{Symbol, Int64}, +) + mass_matrix[global_index[:Vm], global_index[:Vm]] = PSY.get_Tr(avr) + mass_matrix[global_index[:x_d], global_index[:x_d]] = PSY.get_T_da(avr) + return +end + +################################## +##### Differential Equations ##### +################################## + function mdl_avr_ode!( ::AbstractArray{<:ACCEPTED_REAL_TYPES}, ::AbstractArray{<:ACCEPTED_REAL_TYPES}, @@ -634,3 +652,75 @@ function mdl_avr_ode!( inner_vars[Vf_var] = Vf return end + +###################################################################### +function mdl_avr_ode!( + device_states::AbstractArray, + output_ode::AbstractArray, + inner_vars::AbstractArray, + dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, PSY.ST6B, TG, P}}, + h, + t, +) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS} + + #Obtain references + V_ref = get_V_ref(dynamic_device) + + #Obtain avr + avr = PSY.get_avr(dynamic_device) + + #Obtain indices for component w/r to device + local_ix = get_local_state_ix(dynamic_device, typeof(avr)) + + #Define inner states for component + internal_states = @view device_states[local_ix] + Vm = internal_states[1] + x_i = internal_states[2] + x_d = internal_states[3] + Vg = internal_states[4] + + #Define external states for device + V_th = sqrt(inner_vars[VR_gen_var]^2 + inner_vars[VI_gen_var]^2) # machine's terminal voltage + Vs = inner_vars[V_pss_var] # PSS output + Xad_Ifd = inner_vars[Xad_Ifd_var] # machine's field current in exciter base + + #Get parameters + Tr = PSY.get_Tr(avr) + K_pa = PSY.get_K_pa(avr) #k_pa>0 + K_ia = PSY.get_K_ia(avr) + K_da = PSY.get_K_da(avr) + T_da = PSY.get_T_da(avr) + Va_min, Va_max = PSY.get_Va_lim(avr) + K_ff = PSY.get_K_ff(avr) + K_m = PSY.get_K_m(avr) + K_ci = PSY.get_K_ci(avr) #K_cl in pss + K_lr = PSY.get_K_lr(avr) + I_lr = PSY.get_I_lr(avr) + Vr_min, Vr_max = PSY.get_Vr_lim(avr) + Kg = PSY.get_Kg(avr) + Tg = PSY.get_Tg(avr) #T_g>0 + + #Compute block derivatives + _, dVm_dt = low_pass_mass_matrix(V_th, Vm, 1.0, Tr) + pid_input = V_ref + Vs - Vm + pi_out, dx_i = pi_block_nonwindup(pid_input, x_i, K_pa, K_ia, Va_min, Va_max) + pd_out, dx_d = high_pass_mass_matrix(pid_input, x_d, K_da, T_da) + Va = pi_out + pd_out + + ff_out = ((Va - Vg) * K_m) + (K_ff * Va) + V_r1 = max(((I_lr * K_ci) - Xad_Ifd) * K_lr, Vr_min) + V_r2 = clamp(ff_out, Vr_min, Vr_max) + V_r = min(V_r1, V_r2) + E_fd = V_r * Vm + _, dVg = low_pass(E_fd, Vg, Kg, Tg) + + #Compute 4 States AVR ODE: + output_ode[local_ix[1]] = dVm_dt + output_ode[local_ix[2]] = dx_i + output_ode[local_ix[3]] = dx_d + output_ode[local_ix[4]] = dVg + + #Update inner_vars + inner_vars[Vf_var] = E_fd + return +end diff --git a/src/post_processing/post_proc_generator.jl b/src/post_processing/post_proc_generator.jl index e9a63b74f..db8d715fb 100644 --- a/src/post_processing/post_proc_generator.jl +++ b/src/post_processing/post_proc_generator.jl @@ -831,6 +831,69 @@ function _field_voltage( return ts, Vf end +""" +Function to obtain the field voltage time series of a Dynamic Generator with avr ESST1A. + +""" +function _field_voltage( + avr::PSY.ST6B, + name::String, + res::SimulationResults, + dt::Union{Nothing, Float64, Vector{Float64}}, +) + #ASSUMPTION THAT PSS IS NOT ACTIVATE + # Obtain state Vm + ts, Vm = post_proc_state_series(res, (name, :Vm), dt) + + # Obtain state x_i + ts, x_i = post_proc_state_series(res, (name, :x_i), dt) + + # Obtain state x_d + ts, x_d = post_proc_state_series(res, (name, :x_d), dt) + + # Obtain state x_d + ts, Vg = post_proc_state_series(res, (name, :Vg), dt) + + # Obtain state Xad_Ifd + ts, Xad_Ifd = post_proc_field_current_series(res, name, dt) + + # Obtain PSS output + _, Vs = post_proc_pss_output_series(res, name, dt) + + V_ref = PSY.get_V_ref(avr) + + #Get parameters + Tr = PSY.get_Tr(avr) + K_pa = PSY.get_K_pa(avr) #k_pa>0 + K_ia = PSY.get_K_ia(avr) + K_da = PSY.get_K_da(avr) + T_da = PSY.get_T_da(avr) + Va_min, Va_max = PSY.get_Va_lim(avr) + K_ff = PSY.get_K_ff(avr) + K_m = PSY.get_K_m(avr) + K_ci = PSY.get_K_ci(avr) #K_cl in pss + K_lr = PSY.get_K_lr(avr) + I_lr = PSY.get_I_lr(avr) + Vr_min, Vr_max = PSY.get_Vr_lim(avr) + Kg = PSY.get_Kg(avr) + Tg = PSY.get_Tg(avr) #T_g>0 + + Efd = zeros(length(ts)) + for (ix, t) in enumerate(ts) + pid_input = V_ref + Vs[ix] - Vm[ix] + pi_out, dx_i = pi_block_nonwindup(pid_input, x_i[ix], K_pa, K_ia, Va_min, Va_max) + pd_out, dx_d = high_pass_mass_matrix(pid_input, x_d[ix], K_da, T_da) + Va = pi_out + pd_out + ff_out = ((Va - Vg[ix]) * K_m) + (K_ff * Va) + V_r1 = max(((I_lr * K_ci) - Xad_Ifd[ix]) * K_lr, Vr_min) + V_r2 = clamp(ff_out, Vr_min, Vr_max) + V_r = min(V_r1, V_r2) + Efd[ix] = V_r * Vm[ix] + end + + return ts, Efd +end + """ Function to obtain the pss output time series of a Dynamic Generator with pss PSSFixed. diff --git a/test/benchmarks/psse/ST6B/ThreeBusMulti.raw b/test/benchmarks/psse/ST6B/ThreeBusMulti.raw new file mode 100644 index 000000000..7ec76d042 --- /dev/null +++ b/test/benchmarks/psse/ST6B/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 diff --git a/test/benchmarks/psse/ST6B/ThreeBus_ST6B.dyr b/test/benchmarks/psse/ST6B/ThreeBus_ST6B.dyr new file mode 100644 index 000000000..850a4e2af --- /dev/null +++ b/test/benchmarks/psse/ST6B/ThreeBus_ST6B.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.10000 0.80000 / + 102 'ST6B' 1 0 + 0.12000E-01 50.038 45.094 + 0.0000 0.0000 4.8100 + -3.8500 1.0000 1.0000 + 1.0577 17.330 2.1500 + 4.8100 -3.8500 1.0000 0.20000E-01 / diff --git a/test/results/results_eigenvalues.jl b/test/results/results_eigenvalues.jl index b37b79b4c..c7cb82601 100644 --- a/test/results/results_eigenvalues.jl +++ b/test/results/results_eigenvalues.jl @@ -971,3 +971,15 @@ test58_eigvals_GateFlag = [ -0.19849913626541832 + 0.17939599481963514im -0.04238523317744444 + 0.0im ] + +test59_eigvals = [ + -99.69916689593907 + 0.0im + -84.05122388240402 + 0.0im + -39.79684261365684 + 0.0im + -38.3444969887499 + 0.0im + -8.550153616541023 + 0.0im + -0.6772738854081002 - 8.637736379730404im + -0.6772738854081002 + 8.637736379730404im + -0.36823401395037775 - 0.5649228138733879im + -0.36823401395037775 + 0.5649228138733879im +] diff --git a/test/results/results_initial_conditions.jl b/test/results/results_initial_conditions.jl index d54376ab5..cb38572b8 100644 --- a/test/results/results_initial_conditions.jl +++ b/test/results/results_initial_conditions.jl @@ -1781,3 +1781,28 @@ test58_x0_init_GateFlag = Dict( 2.099999999999999 ], ) + +test59_x0_init = Dict( + "V_R" => [ + 1.05 + 1.0198615749401696 + 0.9817964461177487 + ], + "V_I" => [ + 0.0 + -0.01680380791834505 + -0.15145839832533242 + ], + "generator-102-1" => [ + 0.8084146853601003 + 0.5302614530689579 + 0.7288775579091268 + 0.7311890654084042 + 0.9615873122791483 + 1.0 + 1.0199999999999998 + 0.047279162537359226 + 0.0 + 2.153115531256307 + ], +) diff --git a/test/test_case59_st6b.jl b/test/test_case59_st6b.jl new file mode 100644 index 000000000..5ccb57311 --- /dev/null +++ b/test/test_case59_st6b.jl @@ -0,0 +1,157 @@ +""" +Validation PSSE/ST6B: +This case study defines a three bus system with an infinite bus, GENROU and a load. +The GENROU machine has connected an ST6B Excitation System. +The fault drop the line connecting the infinite bus and GENROU. +""" +################################################## +############### LOAD DATA ######################## +################################################## + +raw_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/ST6B/ThreeBusMulti.raw") +dyr_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/ST6B/ThreeBus_ST6B.dyr") +#csv_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/ST6B/results_PSSe.csv") + +@testset "Test 59 ST6B ResidualModel" begin + path = joinpath(pwd(), "test-psse-ST6B") + !isdir(path) && mkdir(path) + try + # Define system + sys = System(raw_file, 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, + (0.0, 20.0), #time span + BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault + ) + + # Test Initial Condition + diff_val = [0.0] + res = get_init_values_for_comparison(sim) + for (k, v) in test59_x0_init + diff_val[1] += LinearAlgebra.norm(res[k] - v) + end + + @test diff_val[1] < 1e-3 + + # Obtain small signal results for initial conditions + small_sig = small_signal_analysis(sim) + eigs = small_sig.eigenvalues + @test small_sig.stable + + # Test Eigenvalues + @test LinearAlgebra.norm(eigs - test59_eigvals) < 1e-3 + + # Solve problem + @test execute!(sim, IDA(); dtmax = 0.005, saveat = 0.005) == + PSID.SIMULATION_FINALIZED + results = read_results(sim) + + # Obtain results + + t_psid, v2_psid = get_voltage_magnitude_series(results, 102) + _, v3_psid = get_voltage_magnitude_series(results, 103) + _, ω_psid = get_state_series(results, ("generator-102-1", :ω)) + _, Vf = get_field_voltage_series(results, "generator-102-1") + + #= + # Obtain PSSE results + M = get_csv_data(csv_file) + t_psse = M[:, 1] + v1_psse = M[:, 2] + v2_psse = M[:, 3] + v3_psse = M[:, 4] + v4_psse = M[:, 5] + ω_psse = M[:, 6] .+ 1.0 + efd_psse = M[:, 7] + + # Test Transient Simulation Results + + @test LinearAlgebra.norm(t_psid - round.(t_psse, digits = 3)) == 0.0 + @test LinearAlgebra.norm(v2_psid - v2_psse, Inf) <= 1e-3 + @test LinearAlgebra.norm(v3_psid - v3_psse, Inf) <= 1e-3 + @test LinearAlgebra.norm(ω_psid - ω_psse, Inf) <= 1e-3 + =# + finally + @info("removing test files") + rm(path; force = true, recursive = true) + end +end + +@testset "Test 59 ST6B MassMatrixModel" begin + path = joinpath(pwd(), "test-psse-ST6B") + !isdir(path) && mkdir(path) + try + # Define system + sys = System(raw_file, 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, + (0.0, 20.0), #time span + BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault + ) + + # Test Initial Condition + diff_val = [0.0] + res = get_init_values_for_comparison(sim) + for (k, v) in test59_x0_init + diff_val[1] += LinearAlgebra.norm(res[k] - v) + end + + @test diff_val[1] < 1e-3 + + # Obtain small signal results for initial conditions + small_sig = small_signal_analysis(sim) + eigs = small_sig.eigenvalues + @test small_sig.stable + + # Test Eigenvalues + @test LinearAlgebra.norm(eigs - test59_eigvals) < 1e-3 + + # Solve problem + @test execute!(sim, Rodas4(); dtmax = 0.005, saveat = 0.005) == + PSID.SIMULATION_FINALIZED + results = read_results(sim) + + # Obtain results + + t_psid, v2_psid = get_voltage_magnitude_series(results, 102) + _, v3_psid = get_voltage_magnitude_series(results, 103) + _, ω_psid = get_state_series(results, ("generator-102-1", :ω)) + _, Vf = get_field_voltage_series(results, "generator-102-1") + + #= + # Obtain PSSE results + M = get_csv_data(csv_file) + t_psse = M[:, 1] + v1_psse = M[:, 2] + v2_psse = M[:, 3] + v3_psse = M[:, 4] + v4_psse = M[:, 5] + ω_psse = M[:, 6] .+ 1.0 + efd_psse = M[:, 7] + + # Test Transient Simulation Results + + @test LinearAlgebra.norm(t_psid - round.(t_psse, digits = 3)) == 0.0 + @test LinearAlgebra.norm(v2_psid - v2_psse, Inf) <= 1e-3 + @test LinearAlgebra.norm(v3_psid - v3_psse, Inf) <= 1e-3 + @test LinearAlgebra.norm(ω_psid - ω_psse, Inf) <= 1e-3 + =# + finally + @info("removing test files") + rm(path; force = true, recursive = true) + end +end