Merge branch 'main' into rh/new_st8c

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

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

@@ -985,6 +985,24 @@ test59_eigvals = [
 ]
 
 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
+]
+
+test61_eigvals = [
  -351.6684302703973 + 0.0im
  -163.4858530295313 + 0.0im
  -98.32270507099068 + 0.0im
@@ -995,5 +1013,5 @@ test60_eigvals = [
  -2.752172823341148 + 8.057381046225618im
  -1.5898057550991769 + 0.0im
  -0.3644442941887857 + 0.0im
- -0.025119358661025975 + 0.0im
+ -0.025119358661025975 + 0.0im 
 ]

test/results/results_initial_conditions.jl

@@ -1807,6 +1807,36 @@ test59_x0_init = Dict(
  ],
 )
 
+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
+ ],
+)
+
 test60_x0_init = Dict(
  "V_R" => [
  1.05