Add wpidhy (#391)

* Add WPIDHY governor model

* Fix typo at the input of the regulator

* Add initialization for WPIDHY Turbine Governor

* Fix typo at initialization

* Add files for testing in 3 machine system. Sign of the reg parameter was corrected in .dyr file, it must be negative (See PSSE documentation)

* Formatter

* Add post processing

* Add Test

* fix rebase issue

* formatter

* change type

---------

Co-authored-by: Rodrigo Henríquez-Auba <rodrigomha@gmail.com>

Author: SebastianManriqueM

src/initialization/generator_components/init_tg.jl

@@ -406,3 +406,100 @@ function initialize_tg!(
     end
     return
 end
+
+function initialize_tg!(
+    device_states,
+    static::PSY.StaticInjection,
+    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, A, PSY.WPIDHY, 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]
+    τe0 = inner_vars[τe_var]
+    P0 = PSY.get_active_power(static)
+
+    #Get parameters
+    tg = PSY.get_prime_mover(dynamic_device)
+    reg = PSY.get_reg(tg)
+    Kp = PSY.get_Kp(tg)
+    Ki = PSY.get_Ki(tg)
+    Kd = PSY.get_Kd(tg)
+    Ta = PSY.get_Ta(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)
+    P_min, P_max = PSY.get_P_lim(tg)
+    Tw = PSY.get_Tw(tg)
+
+    #Compute controller parameters for equivalent TF
+    Kp_prime = (-Ta * Ki) + Kp
+    Kd_prime = (Ta^2 * Ki) - (Ta * Kp) + Kd
+    Ki_prime = Ki
+
+    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]
+
+        pid_input = x_g1
+        pi_out, dxg2_dt = pi_block(pid_input, x_g2, Kp_prime, Ki_prime)
+        pd_out, dxg4_dt = high_pass(pid_input, x_g4, Kd_prime, Ta)
+
+        power_at_gate =
+            three_level_gate_to_power_map(x_g6, gate_openings, power_gate_openings)
+
+        y_LL_out, dxg7_dt = lead_lag(power_at_gate, x_g7, 1.0, -Tw, Tw / 2.0)
+
+        out[1] = y_LL_out - τm0
+        out[2] = (P0 - P_ref) * reg - x_g1
+        out[3] = dxg2_dt
+        out[4] = pi_out + pd_out - x_g3
+        out[5] = dxg4_dt
+        out[6] = x_g3 - x_g5
+        out[7] = x_g5
+        out[8] = dxg7_dt
+    end
+    gate0 = three_level_power_to_gate_map(τm0, gate_openings, power_gate_openings)
+    P_ref_guess = P0
+    xg1_guess = τe0 - P_ref_guess
+
+    #x0 = [P_ref_guess, xg1_guess, 0.0, 0.0, 0.0, 0.0, gate0, 3.0 * τm0]
+    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(
+                "WPIDHY Turbine Governor $(PSY.get_name(static)) $(sol_x0[7]) outside its gate limits $G_min, $G_max. Consider updating the operating point.",
+            )
+        elseif τm0 > P_max || τm0 < P_min
+            @error(
+                "WPIDHY Turbine Governor $(PSY.get_name(static)) $(sol_x0[8]) outside its power limits $P_min, $P_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

src/models/common_controls.jl

@@ -1606,7 +1606,7 @@ function three_level_gate_to_power_map(
     p0 = 0.0
     g3 = 1.0
     if input <= g0
-        return zero(T)
+        return zero(Z)
     elseif g0 < input <= g1
         return ((p1 - p0) / (g1 - g0)) * (input - g0) + p0
     elseif g1 < input <= g2

src/models/generator_models/tg_models.jl

@@ -43,6 +43,15 @@ function mass_matrix_tg_entries!(
     return
 end
 
+function mass_matrix_tg_entries!(
+    mass_matrix,
+    tg::PSY.WPIDHY,
+    global_index::Base.ImmutableDict{Symbol, Int64},
+)
+    mass_matrix[global_index[:x_g1], global_index[:x_g1]] = PSY.get_T_reg(tg)
+    return
+end
+
 ##################################
 ##### Differential Equations #####
 ##################################
@@ -409,6 +418,102 @@ 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.WPIDHY, 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.WPIDHY)
+
+    #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 PID 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.WPIDHY)
+    ω = @view device_states[external_ix]
+
+    # Read Inner Vars
+    τ_e = inner_vars[τe_var]
+
+    #Get Parameters
+    tg = PSY.get_prime_mover(device)
+    T_reg = PSY.get_T_reg(tg)
+    reg = PSY.get_reg(tg)
+    Kp = PSY.get_Kp(tg)
+    Ki = PSY.get_Ki(tg)
+    Kd = PSY.get_Kd(tg)#Kd>0
+    Ta = PSY.get_Ta(tg)
+    Tb = PSY.get_Tb(tg)
+    V_min, V_max = PSY.get_V_lim(tg)
+    G_min, G_max = PSY.get_G_lim(tg)
+    P_min, P_max = PSY.get_P_lim(tg)
+    Tw = PSY.get_Tw(tg)
+    D = PSY.get_D(tg)
+    gate_openings = PSY.get_gate_openings(tg)
+    power_gate_openings = PSY.get_power_gate_openings(tg)
+
+    #Compute controller parameters for equivalent TF
+    Kp_prime = (-Ta * Ki) + Kp
+    Kd_prime = (Ta^2 * Ki) - (Ta * Kp) + Kd
+    Ki_prime = Ki
+
+    x_in = τ_e - P_ref
+    #Compute block derivatives
+    _, dxg1_dt = low_pass_mass_matrix(x_in, x_g1, reg, T_reg)
+    pid_input = x_g1 - (ω[1] - ω_sys)
+    pi_out, dxg2_dt = pi_block(pid_input, x_g2, Kp_prime, Ki_prime)
+    pd_out, dxg4_dt = high_pass(pid_input, x_g4, Kd_prime, Ta)
+    pid_out = pi_out + pd_out
+    _, dxg3_dt = low_pass(pid_out, x_g3, 1.0, Ta)
+    _, dxg5_dt = low_pass(x_g3, x_g5, 1.0, Tb)
+
+    #Set clamping for G_vel.
+    G_vel_sat = clamp(x_g5, V_min, V_max)
+
+    # Compute integrator
+    xg6_sat, dxg6_dt = integrator_windup(G_vel_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
+    ll_out, dxg7_dt = lead_lag(power_at_gate, x_g7, 1.0, -Tw, Tw / 2.0)
+    Power_sat = clamp(ll_out, P_min, P_max)
+
+    #Compute output torque
+    P_m = Power_sat - (D * (ω[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
+
 function mdl_tg_ode!(
     device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
     output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},

src/post_processing/post_proc_generator.jl

@@ -1179,3 +1179,34 @@ function _mechanical_torque(
     τm = Pe ./ ω
     return ts, τm
 end
+
+function _mechanical_torque(
+    tg::PSY.WPIDHY,
+    name::String,
+    res::SimulationResults,
+    dt::Union{Nothing, Float64, Vector{Float64}},
+)
+    # Get params
+    D = PSY.get_D(tg)
+    gate_openings = PSY.get_gate_openings(tg)
+    power_gate_openings = PSY.get_power_gate_openings(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)
+    Pm = 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)
+        ll_out, _ = lead_lag(power_at_gate, x7, 1.0, -Tw, Tw / 2.0)
+        Pm[ix] = ll_out - D * (ω[ix] - ω_ref)
+    end
+    τm = Pm ./ ω
+    return ts, τm
+end

test/benchmarks/psse/WPIDHY/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/WPIDHY/ThreeBus_WPIDHY.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 'WPIDHY'   1    0.10000E-01  -0.50000E-01   3.0000      0.50000
+          3.0000      0.20000E-01  0.10000E-01  0.25000E-01 -0.25000E-01
+          1.0000      0.20000       3.2500       52001.       0.0000
+          0.0000      0.12000      0.54500      0.40000      0.82700
+         0.83000       1.0000       /

test/results/results_eigenvalues.jl

@@ -983,3 +983,21 @@ test59_eigvals = [
     -0.36823401395037775 - 0.5649228138733879im
     -0.36823401395037775 + 0.5649228138733879im
 ]
+
+test60_eigvals = [
+    -106.13231871102334 + 0.0im
+    -91.13024819676329 + 0.0im
+    -63.8688081304463 + 0.0im
+    -40.03565106938418 + 0.0im
+    -38.75323547447537 - 8.125818822049656im
+    -38.75323547447537 + 8.125818822049656im
+    -6.104826779980879 + 0.0im
+    -1.395127240875047 - 5.966607970662552im
+    -1.395127240875047 + 5.966607970662552im
+    -0.9449531063097173 + 0.0im
+    -0.48179149983179287 + 0.0im
+    -0.1941504953444676 - 0.18213498194497163im
+    -0.1941504953444676 + 0.18213498194497163im
+    -0.03057660237077197 - 0.30118043163774183im
+    -0.03057660237077197 + 0.30118043163774183im
+]

test/results/results_initial_conditions.jl

@@ -1806,3 +1806,33 @@ 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.7538967192769663
+        0.5555487379562241
+        0.7042452148052648
+        0.7246287886385532
+        0.9158416020737494
+        1.0
+        1.4986692863524897
+        0.044960078577949814
+        0.0
+        3.871458709170383e-12
+        3.871569731472846e-12
+        0.0
+        3.871569731472846e-12
+        0.7417441860465115
+        2.099999999999999
+    ],
+)