ESST1A model and test (#348)

* Add ESST1A and test

* Fix format and error

* Fix wrong variable name

src/PowerSimulationsDynamics.jl

@@ -48,6 +48,7 @@ export get_source_real_current_series
 export get_source_imaginary_current_series
 export get_field_current_series
 export get_field_voltage_series
+export get_pss_output_series
 export get_mechanical_torque_series
 export get_real_current_branch_flow
 export get_imaginary_current_branch_flow

src/base/simulation_results.jl

@@ -206,6 +206,23 @@ function post_proc_field_voltage_series(
     return ts, Vf
 end
 
+"""
+Function to compute the pss output time series of a Dynamic Injection series out of the DAE Solution. It receives the solution and the string
+name of the Dynamic Injection device.
+
+"""
+function post_proc_pss_output_series(
+    res::SimulationResults,
+    name::String,
+    dt::Union{Nothing, Float64},
+)
+    system = get_system(res)
+    device = PSY.get_component(PSY.StaticInjection, system, name)
+    dyn_device = PSY.get_dynamic_injector(device)
+    ts, Vs = compute_pss_output(res, dyn_device, dt)
+    return ts, Vs
+end
+
 """
 Function to compute the mechanical torque output time series of a Dynamic Injection series out of the DAE Solution. It receives the solution and the
 string name of the Dynamic Injection device.
@@ -430,6 +447,23 @@ function get_field_voltage_series(res::SimulationResults, name::String; dt = not
     return post_proc_field_voltage_series(res, name, dt)
 end
 
+"""
+    get_pss_output_series(
+            res::SimulationResults,
+            name::String,
+    )
+
+Function to obtain the pss output time series of a Dynamic Generator out of the DAE Solution.
+
+# Arguments
+
+- `res::SimulationResults` : Simulation Results object that contains the solution
+- `name::String` : Name to identify the specified device
+"""
+function get_pss_output_series(res::SimulationResults, name::String; dt = nothing)
+    return post_proc_pss_output_series(res, name, dt)
+end
+
 """
     get_mechanical_torque_series(
             res::SimulationResults,

src/initialization/generator_components/init_avr.jl

@@ -512,3 +512,62 @@ function initialize_avr!(
     avr_states[5] = Vr3 #Vr3
     return
 end
+
+function initialize_avr!(
+    device_states,
+    static::PSY.StaticInjection,
+    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, PSY.ESST1A, 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
+    Vt = sqrt(inner_vars[VR_gen_var]^2 + inner_vars[VI_gen_var]^2)
+    #Obtain field winding current 
+    Ifd = inner_vars[Xad_Ifd_var] # machine's field current in exciter base (for the available generator models)
+
+    #Get parameters
+    avr = PSY.get_avr(dynamic_device)
+    Tc = PSY.get_Tc(avr)
+    Tb = PSY.get_Tb(avr)
+    Tc1 = PSY.get_Tc1(avr)
+    Tb1 = PSY.get_Tb1(avr)
+    Ka = PSY.get_Ka(avr)
+    Vr_min, Vr_max = PSY.get_Vr_lim(avr)
+    Kc = PSY.get_Kc(avr)
+    Kf = PSY.get_Kf(avr)
+    Tf = PSY.get_Tf(avr)
+    K_lr = PSY.get_K_lr(avr)
+    I_lr = PSY.get_I_lr(avr)
+
+    # Check limits to field voltage 
+    if (Vt * Vr_min > Vf0) || (Vf0 > Vt * Vr_max - Kc * Ifd)
+        @error(
+            "Field Voltage for AVR in $(PSY.get_name(dynamic_device)) is $(Vf0) pu, which is outside its limits.  Consider updating the operating point."
+        )
+    end
+
+    #Compute auxiliary parameters
+    Itemp = K_lr * (Ifd - I_lr)
+    Iresult = Itemp > 0 ? Itemp : 0
+
+    Va = Vf0 + Iresult
+
+    #Update V_ref
+    Vref0 = Vt + Va / Ka
+
+    PSY.set_V_ref!(avr, Vref0)
+    set_V_ref(dynamic_device, Vref0)
+
+    #States of ESST1A_PTI are Vm, Vr1, Vr2, Va, Vr3
+
+    #Update AVR states
+    avr_ix = get_local_state_ix(dynamic_device, PSY.ESST1A)
+    avr_states = @view device_states[avr_ix]
+    avr_states[1] = Vt
+    avr_states[2] = (1 - Tc / Tb) * Va / Ka
+    avr_states[3] = (1 - Tc1 / Tb1) * Va / Ka
+    avr_states[4] = Va
+    avr_states[5] = -Kf / Tf * Vf0
+    return
+end

src/models/generator_models/avr_models.jl

@@ -49,6 +49,18 @@ function mass_matrix_avr_entries!(
     return
 end
 
+function mass_matrix_avr_entries!(
+    mass_matrix,
+    avr::PSY.ESST1A,
+    global_index::Base.ImmutableDict{Symbol, Int64},
+)
+    mass_matrix[global_index[:Vm], global_index[:Vm]] = PSY.get_Tr(avr)
+    mass_matrix[global_index[:Vr1], global_index[:Vr1]] = PSY.get_Tb(avr)
+    mass_matrix[global_index[:Vr2], global_index[:Vr2]] = PSY.get_Tb1(avr)
+    mass_matrix[global_index[:Va], global_index[:Va]] = PSY.get_Ta(avr)
+    return
+end
+
 function mdl_avr_ode!(
     ::AbstractArray{<:ACCEPTED_REAL_TYPES},
     ::AbstractArray{<:ACCEPTED_REAL_TYPES},
@@ -523,3 +535,83 @@ 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.ESST1A, TG, P}},
+) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}
+
+    #Obtain references
+    V0_ref = get_V_ref(dynamic_device)
+
+    #Obtain indices for component w/r to device
+    local_ix = get_local_state_ix(dynamic_device, PSY.ESST1A)
+
+    #Define inner states for component
+    internal_states = @view device_states[local_ix]
+    Vm = internal_states[1]
+    Vr1 = internal_states[2]
+    Vr2 = internal_states[3]
+    Va = internal_states[4]
+    Vr3 = internal_states[5]
+
+    #Define external states for device
+    Vt = 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 
+    Ifd = inner_vars[Xad_Ifd_var] # machine's field current in exciter base 
+
+    #Get parameters
+    avr = PSY.get_avr(dynamic_device)
+    UEL = PSY.get_UEL_flags(avr)
+    VOS = PSY.get_PSS_flags(avr)
+    Tr = PSY.get_Tr(avr)
+    Vi_min, Vi_max = PSY.get_Vi_lim(avr)
+    Tc = PSY.get_Tc(avr)
+    Tb = PSY.get_Tb(avr)
+    Tc1 = PSY.get_Tc1(avr)
+    Tb1 = PSY.get_Tb1(avr)
+    Ka = PSY.get_Ka(avr)
+    Ta = PSY.get_Ta(avr)
+    Va_min, Va_max = PSY.get_Va_lim(avr)
+    Vr_min, Vr_max = PSY.get_Vr_lim(avr)
+    Kc = PSY.get_Kc(avr)
+    Kf = PSY.get_Kf(avr)
+    Tf = PSY.get_Tf(avr)
+    K_lr = PSY.get_K_lr(avr)
+    I_lr = PSY.get_I_lr(avr)
+
+    #Compute auxiliary parameters
+    Itemp = K_lr * (Ifd - I_lr)
+    Iresult = Itemp > 0.0 ? Itemp : 0.0
+
+    if VOS == 1
+        V_ref = V0_ref + Vs
+        Va_sum = Va - Iresult
+    elseif VOS == 2
+        V_ref = V0_ref
+        Va_sum = Va - Iresult + Vs
+    end
+
+    # Compute block derivatives
+    _, dVm_dt = low_pass_mass_matrix(Vt, Vm, 1.0, Tr)
+    y_hp, dVr3_dt = high_pass(Va_sum, Vr3, Kf, Tf)
+    y_ll1, dVr1_dt =
+        lead_lag_mass_matrix(clamp(V_ref - Vm - y_hp, Vi_min, Vi_max), Vr1, 1.0, Tc, Tb)
+    y_ll2, dVr2_dt =
+        lead_lag_mass_matrix(y_ll1, Vr2, 1.0, Tc1, Tb1)
+    _, dVa_dt = low_pass_nonwindup_mass_matrix(y_ll2, Va, Ka, Ta, Va_min, Va_max)
+
+    #Compute 5 States AVR ODE:
+    output_ode[local_ix[1]] = dVm_dt
+    output_ode[local_ix[2]] = dVr1_dt
+    output_ode[local_ix[3]] = dVr2_dt
+    output_ode[local_ix[4]] = dVa_dt
+    output_ode[local_ix[5]] = dVr3_dt
+
+    #Update inner_vars
+    Vf = clamp(Va_sum, Vt * Vr_min, Vt * Vr_max - Kc * Ifd)
+    inner_vars[Vf_var] = Vf
+    return
+end

src/post_processing/post_proc_generator.jl

@@ -179,6 +179,22 @@ function compute_field_voltage(
     return _field_voltage(avr, PSY.get_name(dynamic_device), res, dt)
 end
 
+"""
+Function to obtain the pss output time series of a Dynamic Generator model out of the DAE Solution. It receives the simulation inputs,
+the dynamic device and bus voltage. It is dispatched for device type to compute the specific output.
+
+"""
+function compute_pss_output(
+    res::SimulationResults,
+    dynamic_device::G,
+    dt::Union{Nothing, Float64},
+) where {G <: PSY.DynamicGenerator}
+
+    #Get PSS
+    pss = PSY.get_pss(dynamic_device)
+    return _pss_output(pss, PSY.get_name(dynamic_device), res, dt)
+end
+
 """
 Function to obtain the mechanical torque time series of a Dynamic Generator model out of the DAE Solution. It receives the simulation inputs,
 the dynamic device and bus voltage. It is dispatched for device type to compute the specific torque.
@@ -753,6 +769,149 @@ 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.ESST1A,
+    name::String,
+    res::SimulationResults,
+    dt::Union{Nothing, Float64},
+)
+    # Obtain state Va
+    ts, Va = post_proc_state_series(res, (name, :Va), dt)
+
+    # Obtain field current
+    _, Ifd = post_proc_field_current_series(res, name, dt)
+
+    # Obtain PSS output
+    _, Vs = post_proc_pss_output_series(res, name, dt)
+
+    # Get parameters
+    VOS = PSY.get_PSS_flags(avr)
+    Vr_min, Vr_max = PSY.get_Vr_lim(avr)
+    Kc = PSY.get_Kc(avr)
+    K_lr = PSY.get_K_lr(avr)
+    I_lr = PSY.get_I_lr(avr)
+
+    # Obtain machine's terminal voltage
+    bus_str = split(name, "-")[2]
+    bus_num = parse(Int, bus_str)
+    _, Vt = get_voltage_magnitude_series(res, bus_num; dt = dt)
+
+    # Compute auxiliary arrays
+    Itemp = similar(Ifd)
+    Iresult = similar(Ifd)
+
+    for (ix, Ifd_ix) in enumerate(Ifd)
+        Itemp[ix] = K_lr * (Ifd_ix - I_lr)
+        Iresult[ix] = Itemp[ix] > 0.0 ? Itemp[ix] : 0.0
+    end
+
+    # Compute Va_sum
+    Va_sum = similar(Va)
+
+    if VOS == 1
+        for (ix, Va_ix) in enumerate(Va)
+            Va_sum[ix] = Va_ix - Iresult[ix]
+        end
+    elseif VOS == 2
+        for (ix, Va_ix) in enumerate(Va)
+            Va_sum[ix] = Va_ix - Iresult[ix] + Vs[ix]
+        end
+    end
+
+    # Compute field voltage
+    Vf = similar(Va_sum)
+    for (ix, Va_sum_ix) in enumerate(Va_sum)
+        Vf[ix] = clamp(Va_sum_ix, Vt[ix] * Vr_min, Vt[ix] * Vr_max - Kc * Ifd[ix])
+    end
+
+    return ts, Vf
+end
+
+"""
+Function to obtain the pss output time series of a Dynamic Generator with pss PSSFixed.
+
+"""
+function _pss_output(
+    pss::PSY.PSSFixed,
+    name::String,
+    res::SimulationResults,
+    dt::Union{Nothing, Float64},
+)
+    ts, _ = post_proc_state_series(res, (name, :δ), dt)
+    return ts, zeros(length(ts))
+end
+
+"""
+Function to obtain the pss output time series of a Dynamic Generator with pss IEEEST.
+
+"""
+function _pss_output(
+    pss::PSY.IEEEST,
+    name::String,
+    res::SimulationResults,
+    dt::Union{Nothing, Float64},
+)
+    # Obtain states
+    ts, x_p2 = post_proc_state_series(res, (name, :x_p2), dt)
+    _, x_p3 = post_proc_state_series(res, (name, :x_p3), dt)
+    _, x_p4 = post_proc_state_series(res, (name, :x_p4), dt)
+    _, x_p5 = post_proc_state_series(res, (name, :x_p5), dt)
+    _, x_p6 = post_proc_state_series(res, (name, :x_p6), dt)
+    _, x_p7 = post_proc_state_series(res, (name, :x_p7), dt)
+
+    # Get Parameters
+    A1 = PSY.get_A1(pss)
+    A2 = PSY.get_A2(pss)
+    A5 = PSY.get_A5(pss)
+    A6 = PSY.get_A6(pss)
+    T1 = PSY.get_T1(pss)
+    T2 = PSY.get_T2(pss)
+    T3 = PSY.get_T3(pss)
+    T4 = PSY.get_T4(pss)
+    T5 = PSY.get_T5(pss)
+    T6 = PSY.get_T6(pss)
+    Ks = PSY.get_Ks(pss)
+    Ls_min, Ls_max = PSY.get_Ls_lim(pss)
+    V_cu = PSY.get_Vcu(pss)
+    V_cl = PSY.get_Vcl(pss)
+
+    # Compute intermediate variables
+    y_f = similar(x_p2)
+    y_LL1 = similar(x_p2)
+    y_LL2 = similar(x_p2)
+    y_out = similar(x_p2)
+
+    for (ix, _) in enumerate(x_p2)
+        y_f[ix], _, _ =
+            lead_lag_2nd_mass_matrix(x_p2[ix], x_p3[ix], x_p4[ix], A1, A2, A5, A6)
+        y_LL1[ix], _ = lead_lag_mass_matrix(y_f[ix], x_p5[ix], 1.0, T1, T2)
+        y_LL2[ix], _ = lead_lag_mass_matrix(y_LL1[ix], x_p6[ix], 1.0, T3, T4)
+        y_out[ix], _ = high_pass_mass_matrix(y_LL2[ix], x_p7[ix], Ks * T5, T6)
+    end
+
+    # Compute V_ss
+    V_ss = clamp.(y_out, Ls_min, Ls_max)
+
+    # Obtain machine's terminal voltage
+    bus_str = split(name, "-")[2]
+    bus_num = parse(Int, bus_str)
+    # TODO: Figure out how to compensate terminal voltage
+    _, V_ct = get_voltage_magnitude_series(res, bus_num; dt = dt)
+
+    # Compute PSS output signal
+    V_pss = similar(V_ss)
+
+    for (ix, V_ss_ix) in enumerate(V_ss)
+        V_pss[ix] = output_pss_limiter(V_ss_ix, V_ct[ix], V_cl, V_cu)
+    end
+
+    return ts, V_pss
+end
+
 """
 Function to obtain the mechanical torque time series of a Dynamic Generator with TGFixed Turbine Governor.