Add PIDGOV model (#387)

* add gate to power function for PIDGOV

* update inner vars to update electric torque

* add PIDGOV model

* [WIP] add test PIDGOV

* add initialization routine for pidgov

* add inverse map for gate to power

* fix error for frequency difference

* add post calculation for mechanical torque for pidgov

* add test58

src/initialization/generator_components/init_tg.jl

@@ -311,3 +311,98 @@ function initialize_tg!(
     end
     return
 end
+
+function initialize_tg!(
+    device_states,
+    static::PSY.StaticInjection,
+    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, A, PSY.PIDGOV, 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]
+    P0 = PSY.get_active_power(static)
+    Δω = 0.0
+    #Get parameters
+    tg = PSY.get_prime_mover(dynamic_device)
+    feedback_flag = PSY.get_feedback_flag(tg)
+    Rperm = PSY.get_Rperm(tg)
+    Kp = PSY.get_Kp(tg)
+    Ki = PSY.get_Ki(tg)
+    Kd = PSY.get_Kd(tg)
+    Ta = PSY.get_Ta(tg)
+    Tb = PSY.get_Tb(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)
+    A_tw = PSY.get_A_tw(tg)
+    Tw = PSY.get_Tw(tg)
+
+    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]
+
+        if feedback_flag == 0
+            P_in = P_ref - P0
+        else
+            P_in = P_ref - x_g6
+        end
+        pid_input = x_g1
+        pi_out, dxg2_dt = pi_block(pid_input, x_g2, Kp, Ki)
+        pd_out, dxg4_dt = high_pass_mass_matrix(pid_input, x_g4, Kd, Ta)
+
+        integrator_input = (1.0 / Tb) * (x_g5 - x_g6)
+        power_at_gate =
+            three_level_gate_to_power_map(x_g6, gate_openings, power_gate_openings)
+        Tz = A_tw * Tw
+        y_LL_out, dxg7_dt = lead_lag(power_at_gate, x_g7, 1.0, -Tz, Tz / 2.0)
+
+        out[1] = y_LL_out - τm0
+        out[2] = (Rperm * P_in - x_g1)
+        out[3] = dxg2_dt
+        out[4] = pi_out - x_g3
+        out[5] = dxg4_dt
+        out[6] = x_g3 + pd_out - x_g5
+        out[7] = integrator_input
+        out[8] = dxg7_dt
+    end
+    gate0 = three_level_power_to_gate_map(τm0, gate_openings, power_gate_openings)
+    if feedback_flag == 0
+        P_ref_guess = P0
+    else
+        P_ref_guess = gate0
+    end
+    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(
+                "PIDGOV Turbine Governor $(PSY.get_name(static)) $(sol_x0[4]) outside its gate limits $G_min, $G_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

@@ -1591,3 +1591,51 @@ function integrator_nonwindup(
     y_sat, dydt_scaled = integrator_nonwindup_mass_matrix(u, y, K, T, y_min, y_max)
     return y_sat, (1.0 / T) * dydt_scaled
 end
+
+"""
+Three level gate to power transformation for PIDGOV    
+(G0, 0), (G1, P1), (G2, P2), (1, P3) are (x,y) coordinates of Power vs Gate functions
+"""
+function three_level_gate_to_power_map(
+    input::Z,
+    gates::Tuple{Float64, Float64, Float64},
+    powers::Tuple{Float64, Float64, Float64},
+) where {Z <: ACCEPTED_REAL_TYPES}
+    g0, g1, g2 = gates
+    p1, p2, p3 = powers
+    p0 = 0.0
+    g3 = 1.0
+    if input <= g0
+        return zero(T)
+    elseif g0 < input <= g1
+        return ((p1 - p0) / (g1 - g0)) * (input - g0) + p0
+    elseif g1 < input <= g2
+        return ((p2 - p1) / (g2 - g1)) * (input - g1) + p1
+    elseif g2 < input
+        return ((p3 - p2) / (g3 - g2)) * (input - g2) + p2
+    end
+end
+
+"""
+Inverse Three level gate to power transformation for PIDGOV    
+(G0, 0), (G1, P1), (G2, P2), (1, P3) are (x,y) coordinates of Power vs Gate functions
+"""
+function three_level_power_to_gate_map(
+    power_input::Z,
+    gates::Tuple{Float64, Float64, Float64},
+    powers::Tuple{Float64, Float64, Float64},
+) where {Z <: ACCEPTED_REAL_TYPES}
+    g0, g1, g2 = gates
+    p1, p2, p3 = powers
+    p0 = 0.0
+    g3 = 1.0
+    if power_input <= p0
+        return g0
+    elseif p0 < power_input <= p1
+        return ((g1 - g0) / (p1 - p0)) * (power_input - p0) + g0
+    elseif p1 < power_input <= p2
+        return ((g2 - g1) / (p2 - p1)) * (power_input - p1) + g1
+    elseif p2 < power_input
+        return ((g3 - g2) / (p3 - p2)) * (power_input - p2) + g2
+    end
+end

src/models/device.jl

@@ -333,10 +333,14 @@ function _update_inner_vars!(
     I_d =
         (1.0 / (R^2 + Xq_pp * Xd_pp)) *
         (-R * (V_dq[1] - ψq_pp) + Xq_pp * (-V_dq[2] + ψd_pp))
+    I_q =
+        (1.0 / (R^2 + Xq_pp * Xd_pp)) * (Xd_pp * (V_dq[1] - ψq_pp) + R * (-V_dq[2] + ψd_pp))
     Se = saturation_function(machine, ψ_pp)
     Xad_Ifd = eq_p + (Xd - Xd_p) * (γ_d1 * I_d - γ_d2 * ψ_kd + γ_d2 * eq_p) + Se * ψd_pp
+    τ_e = I_d * (V_dq[1] + I_d * R) + I_q * (V_dq[2] + I_q * R)
 
     #Update Xad_Ifd
+    inner_vars[τe_var] = τ_e
     inner_vars[Xad_Ifd_var] = Xad_Ifd
     return
 end
@@ -386,11 +390,14 @@ function _update_inner_vars!(
 
     #Currents
     I_d = (1.0 / (R^2 + Xd_pp^2)) * (-R * (V_d + ψq_pp) + Xd_pp * (ψd_pp - V_q))
+    I_q = (1.0 / (R^2 + Xd_pp^2)) * (Xd_pp * (V_d + ψq_pp) + R * (ψd_pp - V_q))
     Se = saturation_function(machine, eq_p)
     Xad_Ifd =
         eq_p + Se * eq_p + (Xd - Xd_p) * (I_d + γ_d2 * (eq_p - ψ_kd - (Xd_p - Xl) * I_d))
+    τ_e = I_d * (V_d + I_d * R) + I_q * (V_q + I_q * R)
 
     #Update Xad_Ifd
+    inner_vars[τe_var] = τ_e
     inner_vars[Xad_Ifd_var] = Xad_Ifd
     return
 end
@@ -442,11 +449,14 @@ function _update_inner_vars!(
 
     #Currents
     I_d = (1.0 / (R^2 + Xd_pp^2)) * (-R * (V_d - ψq_pp) + Xq_pp * (-V_q + ψd_pp))
+    I_q = (1.0 / (R^2 + Xd_pp^2)) * (Xd_pp * (V_d - ψq_pp) + R * (-V_q + ψd_pp))
     Se = saturation_function(machine, ψ_pp)
     Xad_Ifd =
         eq_p + Se * ψd_pp + (Xd - Xd_p) * (I_d + γ_d2 * (eq_p - ψ_kd - (Xd_p - Xl) * I_d))
+    τ_e = I_d * (V_d + I_d * R) + I_q * (V_q + I_q * R)
 
     #Update Xad_Ifd
+    inner_vars[τe_var] = τ_e
     inner_vars[Xad_Ifd_var] = Xad_Ifd
     return
 end

src/models/generator_models/tg_models.jl

@@ -1,3 +1,7 @@
+##################################
+###### Mass Matrix Entries #######
+##################################
+
 function mass_matrix_tg_entries!(
     mass_matrix,
     tg::TG,
@@ -27,6 +31,22 @@ function mass_matrix_tg_entries!(
     return
 end
 
+function mass_matrix_tg_entries!(
+    mass_matrix,
+    tg::PSY.PIDGOV,
+    global_index::Base.ImmutableDict{Symbol, Int64},
+)
+    mass_matrix[global_index[:x_g1], global_index[:x_g1]] = PSY.get_T_reg(tg)
+    mass_matrix[global_index[:x_g3], global_index[:x_g3]] = PSY.get_Ta(tg)
+    mass_matrix[global_index[:x_g4], global_index[:x_g4]] = PSY.get_Ta(tg)
+    mass_matrix[global_index[:x_g5], global_index[:x_g5]] = PSY.get_Ta(tg)
+    return
+end
+
+##################################
+##### Differential Equations #####
+##################################
+
 function mdl_tg_ode!(
     device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
     ::AbstractArray{<:ACCEPTED_REAL_TYPES},
@@ -388,3 +408,100 @@ 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.PIDGOV, 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.PIDGOV)
+
+    #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 PI 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.PIDGOV)
+    ω = @view device_states[external_ix]
+
+    # Read Inner Vars
+    τ_e = inner_vars[τe_var]
+
+    #Get Parameters
+    tg = PSY.get_prime_mover(device)
+    feedback_flag = PSY.get_feedback_flag(tg)
+    Rperm = PSY.get_Rperm(tg)
+    T_reg = PSY.get_T_reg(tg)
+    Kp = PSY.get_Kp(tg)
+    Ki = PSY.get_Ki(tg)
+    Kd = PSY.get_Kd(tg)
+    Ta = PSY.get_Ta(tg)
+    Tb = PSY.get_Tb(tg)
+    D_turb = PSY.get_D_turb(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)
+    A_tw = PSY.get_A_tw(tg)
+    Tw = PSY.get_Tw(tg)
+    V_min, V_max = PSY.get_V_lim(tg)
+
+    #Compute auxiliary parameters
+    if feedback_flag == 0
+        x_in = P_ref - τ_e
+    else
+        x_in = P_ref - x_g6
+    end
+
+    #Compute block derivatives
+    _, dxg1_dt = low_pass_mass_matrix(x_in, x_g1, Rperm, T_reg)
+    pid_input = x_g1 - (ω[1] - ω_sys)
+    pi_out, dxg2_dt = pi_block(pid_input, x_g2, Kp, Ki)
+    _, dxg3_dt = low_pass_mass_matrix(pi_out, x_g3, 1.0, Ta)
+    pd_out, dxg4_dt = high_pass_mass_matrix(pid_input, x_g4, Kd, Ta)
+    _, dxg5_dt = low_pass_mass_matrix(x_g3 + pd_out, x_g5, 1.0, Ta)
+
+    # Compute integrator
+    integrator_input = (1.0 / Tb) * (x_g5 - x_g6)
+    integrator_input_sat = clamp(integrator_input, V_min, V_max)
+    xg6_sat, dxg6_dt =
+        integrator_nonwindup(integrator_input_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
+    Tz = A_tw * Tw
+    ll_out, dxg7_dt = lead_lag(power_at_gate, x_g7, 1.0, -Tz, Tz / 2.0)
+
+    #Compute output torque
+    P_m = ll_out - D_turb * (ω[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

src/post_processing/post_proc_generator.jl

@@ -1081,3 +1081,38 @@ function _mechanical_torque(
     τm = ((x_g4 .- q_nl) .* h * At - D_T * Δω .* x_g3) ./ ω
     return ts, τm
 end
+
+"""
+Function to obtain the mechanical torque time series of a Dynamic Generator with HydroTurbineGov (HYGOV) Turbine Governor.
+
+"""
+function _mechanical_torque(
+    tg::PSY.PIDGOV,
+    name::String,
+    res::SimulationResults,
+    dt::Union{Nothing, Float64, Vector{Float64}},
+)
+    # Get params
+    D_turb = PSY.get_D_turb(tg)
+    gate_openings = PSY.get_gate_openings(tg)
+    power_gate_openings = PSY.get_power_gate_openings(tg)
+    A_tw = PSY.get_A_tw(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)
+    Pe = 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)
+        Tz = A_tw * Tw
+        ll_out, _ = lead_lag(power_at_gate, x7, 1.0, -Tz, Tz / 2.0)
+        Pe[ix] = ll_out - D_turb * (ω[ix] - ω_ref)
+    end
+    τm = Pe ./ ω
+    return ts, τm
+end