Add ST8C, DEGOV1 and TGSimple models. Also, remove Windows PR Testing temporarily

.github/workflows/main-tests.yml

@@ -15,7 +15,7 @@ jobs:
       matrix:
         julia-version: ['1', 'nightly']
         julia-arch: [x64]
-        os: [ubuntu-latest, windows-latest]
+        os: [ubuntu-latest, macOS-latest]
 
     steps:
       - uses: actions/checkout@v2

.github/workflows/pr_testing.yml

@@ -12,7 +12,7 @@ jobs:
       matrix:
         julia-version: ['1']
         julia-arch: [x64]
-        os: [ubuntu-latest, windows-latest, macOS-latest]
+        os: [ubuntu-latest, macOS-latest]
 
     steps:
       - uses: actions/checkout@v2

src/initialization/generator_components/init_avr.jl

@@ -652,3 +652,116 @@
 
     return
 end
+
+function initialize_avr!(
+    device_states,
+    static::PSY.StaticInjection,
+    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, PSY.ST8C, 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 Ifd limiter
+    Ifd = inner_vars[Xad_Ifd_var] # read Lad Ifd (field current times Lad)
+    #Obtain measured terminal voltage
+    Vt0 = sqrt(inner_vars[VR_gen_var]^2 + inner_vars[VI_gen_var]^2)
+
+    #Get parameters
+    avr = PSY.get_avr(dynamic_device)
+    SW1_Flag = PSY.get_SW1_Flag(avr)
+    if SW1_Flag == 1
+        error("Source from generator terminal voltage not supported.")
+    end
+    K_pr = PSY.get_K_pr(avr)
+    K_ir = PSY.get_K_ir(avr)
+    Vpi_min, Vpi_max = PSY.get_Vpi_lim(avr)
+    K_pa = PSY.get_K_pa(avr)
+    K_ia = PSY.get_K_ia(avr)
+    Va_min, Va_max = PSY.get_Va_lim(avr)
+    K_a = PSY.get_K_a(avr)
+    T_a = PSY.get_T_a(avr)
+    Vr_min, Vr_max = PSY.get_Vr_lim(avr)
+    K_f = PSY.get_K_f(avr)
+    T_f = PSY.get_T_f(avr)
+    K_c1 = PSY.get_K_c1(avr)
+    K_p = PSY.get_K_p(avr)
+    K_i2 = PSY.get_K_i2(avr)
+
+    if K_i2 != 0.0
+        error("Feedforward Current for AVR ST8C not implemented yet.")
+    end
+    if K_ia == 0.0
+        error(
+            "Integrator gain cannot be zero for AVR ST8C of $(PSY.get_name(dynamic_device))",
+        )
+    end
+    #To solve V_ref, Vr
+    function f!(out, x)
+        V_ref = x[1]
+        Vm = x[2] # Sensed Voltage
+        x_a1 = x[3] # Regulator Integrator
+        x_a2 = x[4] # Field Regulator
+        x_a3 = x[5] # Controller Integrator
+        x_a4 = x[6] # Regulator Feedback
+
+        #TODO: Implement Terminal Current FF for AVR
+        V_b2 = 0.0
+        #TODO: Implement Voltage Compensation if needed
+        V_e = K_p
+
+        # Compute block derivatives
+        Vs = 0.0
+        V_pi_in = V_ref + Vs - Vm
+        Ifd_ref, dxa1_dt = pi_block_nonwindup(V_pi_in, x_a1, K_pr, K_ir, Vpi_min, Vpi_max)
+        Ifd_diff = Ifd_ref - x_a4
+        pi_out, dxa2_dt = pi_block_nonwindup(Ifd_diff, x_a2, K_pa, K_ia, Va_min, Va_max)
+        _, dxa3_dt = low_pass_nonwindup_mass_matrix(pi_out, x_a3, K_a, T_a, Vr_min, Vr_max)
+        _, dxa4_dt = low_pass_mass_matrix(Ifd, x_a4, K_f, T_f)
+
+        # Compute V_b1
+        I_n1 = K_c1 * Ifd / V_e
+        F_ex = rectifier_function(I_n1)
+        V_b1 = F_ex * V_e
+
+        Efd = V_b1 * x_a3 + V_b2
+
+        #V_ll output first block 
+        out[1] = Efd - Vf0 # we are asking for Vf0 
+        out[2] = Vm - Vt0 # Low Pass Sensor
+        out[3] = dxa1_dt
+        out[4] = dxa2_dt
+        out[5] = dxa3_dt
+        out[6] = dxa4_dt
+    end
+    # Initial Guess
+    V_e0 = K_p
+    I_n10 = K_c1 * Ifd / V_e0
+    F_ex0 = rectifier_function(I_n10)
+    V_b10 = F_ex0 * V_e0
+    x_a30 = Vf0 / V_b10
+    x_a20 = x_a30 / (K_a * K_ia)
+    x0 = [Vt0, Vt0, K_f * Ifd / K_ir, x_a20, x_a30, Ifd * K_f]
+    sol = NLsolve.nlsolve(f!, x0; ftol = STRICT_NLSOLVE_F_TOLERANCE)
+    if !NLsolve.converged(sol)
+        @warn("Initialization of AVR in $(PSY.get_name(static)) failed")
+    else # if converge
+        sol_x0 = sol.zero
+        Vreg = sol_x0[6]
+        if (Vreg > Vr_max) || (Vreg < Vr_min)
+            @error(
+                "Regulator Voltage V_R = $(Vreg) outside the limits for AVR of $(PSY.get_name(dynamic_device)). Consider updating its limits."
+            )
+        end
+        #Update V_ref
+        PSY.set_V_ref!(avr, sol_x0[1])
+        set_V_ref(dynamic_device, sol_x0[1])
+        #Update AVR states
+        avr_ix = get_local_state_ix(dynamic_device, PSY.ST8C)
+        avr_states = @view device_states[avr_ix]
+        avr_states[1] = sol_x0[2]
+        avr_states[2] = sol_x0[3]
+        avr_states[3] = sol_x0[4]
+        avr_states[4] = sol_x0[5]
+        avr_states[5] = sol_x0[6]
+    end
+end

src/models/generator_models/avr_models.jl

@@ -75,6 +75,17 @@
     return
 end
 
+function mass_matrix_avr_entries!(
+    mass_matrix,
+    avr::PSY.ST8C,
+    global_index::Base.ImmutableDict{Symbol, Int64},
+)
+    mass_matrix[global_index[:Vm], global_index[:Vm]] = PSY.get_Tr(avr)
+    mass_matrix[global_index[:x_a3], global_index[:x_a3]] = PSY.get_T_a(avr)
+    mass_matrix[global_index[:x_a4], global_index[:x_a4]] = PSY.get_T_f(avr)
+    return
+end
+
 ##################################
 ##### Differential Equations #####
 ##################################
@@ -724,3 +735,89 @@
     inner_vars[Vf_var] = E_fd
     return
 end
+
+function mdl_avr_ode!(
+    device_states::AbstractArray,
+    output_ode::AbstractArray,
+    inner_vars::AbstractArray,
+    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, PSY.ST8C, 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 indices for component w/r to device
+    local_ix = get_local_state_ix(dynamic_device, PSY.ST8C)
+
+    #Define inner states for component
+    internal_states = @view device_states[local_ix]
+    Vm = internal_states[1] # Sensed Voltage
+    x_a1 = internal_states[2] # Regulator Integrator
+    x_a2 = internal_states[3] # Field Regulator
+    x_a3 = internal_states[4] # Controller Integrator
+    x_a4 = internal_states[5] # Regulator Feedback
+
+    #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)
+    SW1_Flag = PSY.get_SW1_Flag(avr)
+    if SW1_Flag == 1
+        error("Source from generator terminal voltage not supported.")
+    end
+    Tr = PSY.get_Tr(avr)
+    K_pr = PSY.get_K_pr(avr)
+    K_ir = PSY.get_K_ir(avr)
+    Vpi_min, Vpi_max = PSY.get_Vpi_lim(avr)
+    K_pa = PSY.get_K_pa(avr)
+    K_ia = PSY.get_K_ia(avr)
+    Va_min, Va_max = PSY.get_Va_lim(avr)
+    K_a = PSY.get_K_a(avr)
+    T_a = PSY.get_T_a(avr)
+    Vr_min, Vr_max = PSY.get_Vr_lim(avr)
+    K_f = PSY.get_K_f(avr)
+    T_f = PSY.get_T_f(avr)
+    K_c1 = PSY.get_K_c1(avr)
+    K_p = PSY.get_K_p(avr)
+    K_i2 = PSY.get_K_i2(avr)
+    VB1_max = PSY.get_VB1_max(avr)
+
+    if K_i2 != 0.0
+        error("Feedforward Current for AVR ST8C not implemented yet.")
+    end
+    #TODO: Implement Terminal Current FF for AVR
+    V_b2 = 0.0
+    #TODO: Implement Voltage Compensation if needed
+    V_e = K_p
+
+    # Compute block derivatives
+    _, dVm_dt = low_pass_mass_matrix(Vt, Vm, 1.0, Tr)
+    V_pi_in = V_ref + Vs - Vm
+    Ifd_ref, dxa1_dt = pi_block_nonwindup(V_pi_in, x_a1, K_pr, K_ir, Vpi_min, Vpi_max)
+    Ifd_diff = Ifd_ref - x_a4
+    pi_out, dxa2_dt = pi_block_nonwindup(Ifd_diff, x_a2, K_pa, K_ia, Va_min, Va_max)
+    _, dxa3_dt = low_pass_nonwindup_mass_matrix(pi_out, x_a3, K_a, T_a, Vr_min, Vr_max)
+    _, dxa4_dt = low_pass_mass_matrix(Ifd, x_a4, K_f, T_f)
+
+    # Compute V_b1
+    I_n1 = K_c1 * Ifd / V_e
+    F_ex = rectifier_function(I_n1)
+    V_b1 = min(F_ex * V_e, VB1_max)
+
+    #Compute 5 States AVR ODE:
+    output_ode[local_ix[1]] = dVm_dt
+    output_ode[local_ix[2]] = dxa1_dt
+    output_ode[local_ix[3]] = dxa2_dt
+    output_ode[local_ix[4]] = dxa3_dt
+    output_ode[local_ix[5]] = dxa4_dt
+
+    #Update inner_vars
+    Efd = V_b1 * x_a3 + V_b2
+    inner_vars[Vf_var] = Efd
+    return
+end

src/post_processing/post_proc_generator.jl

@@ -832,7 +832,7 @@ function _field_voltage(
 end
 
 """
-Function to obtain the field voltage time series of a Dynamic Generator with avr ESST1A.
+Function to obtain the field voltage time series of a Dynamic Generator with avr ST6B.
 
 """
 function _field_voltage(
@@ -894,6 +894,42 @@ function _field_voltage(
     return ts, Efd
 end
 
+"""
+Function to obtain the field voltage time series of a Dynamic Generator with avr ST8C.
+
+"""
+function _field_voltage(
+    avr::PSY.ST8C,
+    name::String,
+    res::SimulationResults,
+    dt::Union{Nothing, Float64, Vector{Float64}},
+)
+    # Obtain state Vm
+    ts, x_a3 = post_proc_state_series(res, (name, :x_a3), dt)
+
+    # Obtain state Xad_Ifd
+    ts, Xad_Ifd = post_proc_field_current_series(res, name, dt)
+
+    #Get parameters
+    K_c1 = PSY.get_K_c1(avr)
+    K_p = PSY.get_K_p(avr)
+    VB1_max = PSY.get_VB1_max(avr)
+
+    Efd = zeros(length(ts))
+    for (ix, t) in enumerate(ts)
+        #TODO: Get proper V_e and V_b2
+        V_e = K_p
+        V_b2 = 0.0
+        Ifd = Xad_Ifd[ix]
+        I_N1 = K_c1 * Ifd / V_e
+        F_ex = rectifier_function(I_N1)
+        V_b1 = min(F_ex * V_e, VB1_max)
+        Efd[ix] = V_b1 * x_a3[ix] + V_b2
+    end
+
+    return ts, Efd
+end
+
 """
 Function to obtain the pss output time series of a Dynamic Generator with pss PSSFixed.
 

test/benchmarks/psse/ST8C/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

test/benchmarks/psse/ST8C/ThreeBus_ST8C.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 'ST8C'   1       0      0      0      0
+         0.10000E-01    5.0000        5.0000        99.9900        0.0000
+          4.0000        0.1000       99.99620      -99.86600        1.3500
+         0.30000E-02    99.3449       -99.1691       0.50000       0.50000E-02
+         0.10000        4.0000        0.0000        0.0000        0.0000
+          99.000        0.0000        0.0000        99.000        /

test/results/results_eigenvalues.jl

@@ -983,3 +983,17 @@ test59_eigvals = [
     -0.36823401395037775 - 0.5649228138733879im
     -0.36823401395037775 + 0.5649228138733879im
 ]
+
+test60_eigvals = [
+    -351.6684302703973 + 0.0im
+    -163.4858530295313 + 0.0im
+    -98.32270507099068 + 0.0im
+    -49.907783868123786 + 0.0im
+    -39.63477584570041 + 0.0im
+    -10.878948461924733 + 0.0im
+    -2.752172823341148 - 8.057381046225618im
+    -2.752172823341148 + 8.057381046225618im
+    -1.5898057550991769 + 0.0im
+    -0.3644442941887857 + 0.0im
+    -0.025119358661025975 + 0.0im
+]

test/results/results_initial_conditions.jl

@@ -1806,3 +1806,29 @@ 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.7884763261345278
+        0.5085259590582243
+        0.7414384926786595
+        0.6865254295904037
+        0.8631405314160565
+        1.0
+        1.02
+        0.15882529455273622
+        3.0101738941140193
+        0.4063734757053926
+        0.7941264727636811
+    ],
+)