Skip to content

Commit

Permalink
Add wpidhy (#391)
Browse files Browse the repository at this point in the history 
* 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 <[email protected]>
  • Loading branch information
@SebastianManriqueM @rodrigomha
SebastianManriqueM and rodrigomha authored Sep 19, 2024
1 parent 0c375d2 commit c97be58
Show file tree
Hide file tree
Showing 9 changed files with 532 additions and 1 deletion.
97 changes: 97 additions & 0 deletions src/initialization/generator_components/init_tg.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
2 changes: 1 addition & 1 deletion src/models/common_controls.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
105 changes: 105 additions & 0 deletions src/models/generator_models/tg_models.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 #####
##################################
Expand Down Expand Up @@ -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},
Expand Down
31 changes: 31 additions & 0 deletions src/post_processing/post_proc_generator.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
32 changes: 32 additions & 0 deletions test/benchmarks/psse/WPIDHY/ThreeBusMulti.raw
View file
Original file line number Diff line number Diff line change
@@ -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
10 changes: 10 additions & 0 deletions test/benchmarks/psse/WPIDHY/ThreeBus_WPIDHY.dyr
View file
Original file line number Diff line number Diff line change
@@ -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 /
18 changes: 18 additions & 0 deletions test/results/results_eigenvalues.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
]
30 changes: 30 additions & 0 deletions test/results/results_initial_conditions.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
],
)
Loading

0 comments on commit c97be58

Please sign in to comment.