Add AVR ST6B

src/initialization/generator_components/init_avr.jl

@@ -571,3 +571,84 @@
     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

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

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.
 

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

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    /

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
+]

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
+    ],
+)