Skip to content

Commit c5550a7

Browse files
rodrigomharodrigomha
authored
Merge branch 'main' into rh/new_st8c
2 parents ba67ac4 + c97be58 commit c5550a7

9 files changed

Lines changed: 533 additions & 2 deletions

File tree

‎src/initialization/generator_components/init_tg.jl‎

Lines changed: 97 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -406,3 +406,100 @@ function initialize_tg!(
406406
end
407407
return
408408
end
409+
410+
function initialize_tg!(
411+
device_states,
412+
static::PSY.StaticInjection,
413+
dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, A, PSY.WPIDHY, P}},
414+
inner_vars::AbstractVector,
415+
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, P <: PSY.PSS}
416+
417+
#Get mechanical torque to SyncMach
418+
τm0 = inner_vars[τm_var]
419+
τe0 = inner_vars[τe_var]
420+
P0 = PSY.get_active_power(static)
421+
422+
#Get parameters
423+
tg = PSY.get_prime_mover(dynamic_device)
424+
reg = PSY.get_reg(tg)
425+
Kp = PSY.get_Kp(tg)
426+
Ki = PSY.get_Ki(tg)
427+
Kd = PSY.get_Kd(tg)
428+
Ta = PSY.get_Ta(tg)
429+
gate_openings = PSY.get_gate_openings(tg)
430+
power_gate_openings = PSY.get_power_gate_openings(tg)
431+
G_min, G_max = PSY.get_G_lim(tg)
432+
P_min, P_max = PSY.get_P_lim(tg)
433+
Tw = PSY.get_Tw(tg)
434+
435+
#Compute controller parameters for equivalent TF
436+
Kp_prime = (-Ta * Ki) + Kp
437+
Kd_prime = (Ta^2 * Ki) - (Ta * Kp) + Kd
438+
Ki_prime = Ki
439+
440+
function f!(out, x)
441+
P_ref = x[1]
442+
x_g1 = x[2]
443+
x_g2 = x[3]
444+
x_g3 = x[4]
445+
x_g4 = x[5]
446+
x_g5 = x[6]
447+
x_g6 = x[7]
448+
x_g7 = x[8]
449+
450+
pid_input = x_g1
451+
pi_out, dxg2_dt = pi_block(pid_input, x_g2, Kp_prime, Ki_prime)
452+
pd_out, dxg4_dt = high_pass(pid_input, x_g4, Kd_prime, Ta)
453+
454+
power_at_gate =
455+
three_level_gate_to_power_map(x_g6, gate_openings, power_gate_openings)
456+
457+
y_LL_out, dxg7_dt = lead_lag(power_at_gate, x_g7, 1.0, -Tw, Tw / 2.0)
458+
459+
out[1] = y_LL_out - τm0
460+
out[2] = (P0 - P_ref) * reg - x_g1
461+
out[3] = dxg2_dt
462+
out[4] = pi_out + pd_out - x_g3
463+
out[5] = dxg4_dt
464+
out[6] = x_g3 - x_g5
465+
out[7] = x_g5
466+
out[8] = dxg7_dt
467+
end
468+
gate0 = three_level_power_to_gate_map(τm0, gate_openings, power_gate_openings)
469+
P_ref_guess = P0
470+
xg1_guess = τe0 - P_ref_guess
471+
472+
#x0 = [P_ref_guess, xg1_guess, 0.0, 0.0, 0.0, 0.0, gate0, 3.0 * τm0]
473+
x0 = [P_ref_guess, 0.0, gate0, gate0, 0.0, gate0, gate0, 3.0 * τm0]
474+
sol = NLsolve.nlsolve(f!, x0; ftol = STRICT_NLSOLVE_F_TOLERANCE)
475+
if !NLsolve.converged(sol)
476+
@warn("Initialization of Turbine Governor $(PSY.get_name(static)) failed")
477+
else
478+
sol_x0 = sol.zero
479+
#Error if x_g3 is outside PI limits
480+
if sol_x0[7] > G_max || sol_x0[7] < G_min
481+
@error(
482+
"WPIDHY Turbine Governor $(PSY.get_name(static)) $(sol_x0[7]) outside its gate limits $G_min, $G_max. Consider updating the operating point.",
483+
)
484+
elseif τm0 > P_max || τm0 < P_min
485+
@error(
486+
"WPIDHY Turbine Governor $(PSY.get_name(static)) $(sol_x0[8]) outside its power limits $P_min, $P_max. Consider updating the operating point.",
487+
)
488+
end
489+
490+
#Update Control Refs
491+
PSY.set_P_ref!(tg, sol_x0[1])
492+
set_P_ref(dynamic_device, sol_x0[1])
493+
#Update states
494+
tg_ix = get_local_state_ix(dynamic_device, typeof(tg))
495+
tg_states = @view device_states[tg_ix]
496+
tg_states[1] = sol_x0[2]
497+
tg_states[2] = sol_x0[3]
498+
tg_states[3] = sol_x0[4]
499+
tg_states[4] = sol_x0[5]
500+
tg_states[5] = sol_x0[6]
501+
tg_states[6] = sol_x0[7]
502+
tg_states[7] = sol_x0[8]
503+
end
504+
return
505+
end

‎src/models/common_controls.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1606,7 +1606,7 @@ function three_level_gate_to_power_map(
16061606
p0 = 0.0
16071607
g3 = 1.0
16081608
if input <= g0
1609-
return zero(T)
1609+
return zero(Z)
16101610
elseif g0 < input <= g1
16111611
return ((p1 - p0) / (g1 - g0)) * (input - g0) + p0
16121612
elseif g1 < input <= g2

‎src/models/generator_models/tg_models.jl‎

Lines changed: 105 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,15 @@ function mass_matrix_tg_entries!(
4343
return
4444
end
4545

46+
function mass_matrix_tg_entries!(
47+
mass_matrix,
48+
tg::PSY.WPIDHY,
49+
global_index::Base.ImmutableDict{Symbol, Int64},
50+
)
51+
mass_matrix[global_index[:x_g1], global_index[:x_g1]] = PSY.get_T_reg(tg)
52+
return
53+
end
54+
4655
##################################
4756
##### Differential Equations #####
4857
##################################
@@ -409,6 +418,102 @@ function mdl_tg_ode!(
409418
return
410419
end
411420

421+
function mdl_tg_ode!(
422+
device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
423+
output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
424+
inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
425+
ω_sys::ACCEPTED_REAL_TYPES,
426+
device::DynamicWrapper{PSY.DynamicGenerator{M, S, A, PSY.WPIDHY, P}},
427+
h,
428+
t,
429+
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, P <: PSY.PSS}
430+
431+
#Obtain references
432+
P_ref = get_P_ref(device)
433+
434+
#Obtain indices for component w/r to device
435+
local_ix = get_local_state_ix(device, PSY.WPIDHY)
436+
437+
#Define internal states for component
438+
internal_states = @view device_states[local_ix]
439+
x_g1 = internal_states[1] # Filter Input
440+
x_g2 = internal_states[2] # PI Block
441+
x_g3 = internal_states[3] # Regulator After PID Block
442+
x_g4 = internal_states[4] # Derivative (High Pass) Block
443+
x_g5 = internal_states[5] # Second Regulator Block
444+
x_g6 = internal_states[6] # Gate State
445+
x_g7 = internal_states[7] # Water Inertia State
446+
447+
#Obtain external states inputs for component
448+
external_ix = get_input_port_ix(device, PSY.WPIDHY)
449+
ω = @view device_states[external_ix]
450+
451+
# Read Inner Vars
452+
τ_e = inner_vars[τe_var]
453+
454+
#Get Parameters
455+
tg = PSY.get_prime_mover(device)
456+
T_reg = PSY.get_T_reg(tg)
457+
reg = PSY.get_reg(tg)
458+
Kp = PSY.get_Kp(tg)
459+
Ki = PSY.get_Ki(tg)
460+
Kd = PSY.get_Kd(tg)#Kd>0
461+
Ta = PSY.get_Ta(tg)
462+
Tb = PSY.get_Tb(tg)
463+
V_min, V_max = PSY.get_V_lim(tg)
464+
G_min, G_max = PSY.get_G_lim(tg)
465+
P_min, P_max = PSY.get_P_lim(tg)
466+
Tw = PSY.get_Tw(tg)
467+
D = PSY.get_D(tg)
468+
gate_openings = PSY.get_gate_openings(tg)
469+
power_gate_openings = PSY.get_power_gate_openings(tg)
470+
471+
#Compute controller parameters for equivalent TF
472+
Kp_prime = (-Ta * Ki) + Kp
473+
Kd_prime = (Ta^2 * Ki) - (Ta * Kp) + Kd
474+
Ki_prime = Ki
475+
476+
x_in = τ_e - P_ref
477+
#Compute block derivatives
478+
_, dxg1_dt = low_pass_mass_matrix(x_in, x_g1, reg, T_reg)
479+
pid_input = x_g1 - (ω[1] - ω_sys)
480+
pi_out, dxg2_dt = pi_block(pid_input, x_g2, Kp_prime, Ki_prime)
481+
pd_out, dxg4_dt = high_pass(pid_input, x_g4, Kd_prime, Ta)
482+
pid_out = pi_out + pd_out
483+
_, dxg3_dt = low_pass(pid_out, x_g3, 1.0, Ta)
484+
_, dxg5_dt = low_pass(x_g3, x_g5, 1.0, Tb)
485+
486+
#Set clamping for G_vel.
487+
G_vel_sat = clamp(x_g5, V_min, V_max)
488+
489+
# Compute integrator
490+
xg6_sat, dxg6_dt = integrator_windup(G_vel_sat, x_g6, 1.0, 1.0, G_min, G_max)
491+
492+
power_at_gate =
493+
three_level_gate_to_power_map(xg6_sat, gate_openings, power_gate_openings)
494+
495+
# Compute Lead-Lag Block
496+
ll_out, dxg7_dt = lead_lag(power_at_gate, x_g7, 1.0, -Tw, Tw / 2.0)
497+
Power_sat = clamp(ll_out, P_min, P_max)
498+
499+
#Compute output torque
500+
P_m = Power_sat - (D * (ω[1] - ω_sys))
501+
502+
#Compute 1 State TG ODE:
503+
output_ode[local_ix[1]] = dxg1_dt
504+
output_ode[local_ix[2]] = dxg2_dt
505+
output_ode[local_ix[3]] = dxg3_dt
506+
output_ode[local_ix[4]] = dxg4_dt
507+
output_ode[local_ix[5]] = dxg5_dt
508+
output_ode[local_ix[6]] = dxg6_dt
509+
output_ode[local_ix[7]] = dxg7_dt
510+
511+
#Update mechanical torque
512+
inner_vars[τm_var] = P_m / ω[1]
513+
514+
return
515+
end
516+
412517
function mdl_tg_ode!(
413518
device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
414519
output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},

‎src/post_processing/post_proc_generator.jl‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1215,3 +1215,34 @@ function _mechanical_torque(
12151215
τm = Pe ./ ω
12161216
return ts, τm
12171217
end
1218+
1219+
function _mechanical_torque(
1220+
tg::PSY.WPIDHY,
1221+
name::String,
1222+
res::SimulationResults,
1223+
dt::Union{Nothing, Float64, Vector{Float64}},
1224+
)
1225+
# Get params
1226+
D = PSY.get_D(tg)
1227+
gate_openings = PSY.get_gate_openings(tg)
1228+
power_gate_openings = PSY.get_power_gate_openings(tg)
1229+
Tw = PSY.get_Tw(tg)
1230+
setpoints = get_setpoints(res)
1231+
ω_ref = setpoints[name]["ω_ref"]
1232+
1233+
# Get state results
1234+
ts, x_g7 = post_proc_state_series(res, (name, :x_g7), dt)
1235+
_, x_g6 = post_proc_state_series(res, (name, :x_g6), dt)
1236+
_, ω = post_proc_state_series(res, (name, :ω), dt)
1237+
Pm = similar(x_g7)
1238+
1239+
for (ix, x7) in enumerate(x_g7)
1240+
x6 = x_g6[ix]
1241+
power_at_gate =
1242+
three_level_gate_to_power_map(x6, gate_openings, power_gate_openings)
1243+
ll_out, _ = lead_lag(power_at_gate, x7, 1.0, -Tw, Tw / 2.0)
1244+
Pm[ix] = ll_out - D * (ω[ix] - ω_ref)
1245+
end
1246+
τm = Pm ./ ω
1247+
return ts, τm
1248+
end

‎test/benchmarks/psse/WPIDHY/ThreeBusMulti.raw‎

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
0, 100.00, 33, 0, 0, 60.00 / PSS(R)E 33 RAW created by rawd33 TUE, JUL 21 2020 17:55
2+
3+
4+
101,'BUS 1', 138.0000,3, 1, 1, 1,1.05000, 0.0000,1.10000,0.90000,1.10000,0.90000
5+
102,'BUS 2', 138.0000,2, 1, 1, 1,1.02000, -0.9440,1.10000,0.90000,1.10000,0.90000
6+
103,'BUS 3', 138.0000,1, 1, 1, 1,0.99341, -8.7697,1.10000,0.90000,1.10000,0.90000
7+
0 /End of Bus data, Begin Load data
8+
103,'1 ',1, 1, 1, 200.000, 30.000, 0.000, 0.000, 0.000, 0.000, 1,1,0
9+
0 /End of Load data, Begin Fixed shunt data
10+
0 /End of Fixed shunt data, Begin Generator data
11+
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
12+
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
13+
0 /End of Generator data, Begin Branch data
14+
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
15+
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
16+
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
17+
0 /End of Branch data, Begin Transformer data
18+
0 /End of Transformer data, Begin Area interchange data
19+
0 /End of Area interchange data, Begin Two-terminal dc line data
20+
0 /End of Two-terminal dc line data, Begin VSC dc line data
21+
0 /End of VSC dc line data, Begin Impedance correction table data
22+
0 /End of Impedance correction table data, Begin Multi-terminal dc line data
23+
0 /End of Multi-terminal dc line data, Begin Multi-section line data
24+
0 /End of Multi-section line data, Begin Zone data
25+
0 /End of Zone data, Begin Inter-area transfer data
26+
0 /End of Inter-area transfer data, Begin Owner data
27+
0 /End of Owner data, Begin FACTS device data
28+
0 /End of FACTS device data, Begin Switched shunt data
29+
0 /End of Switched shunt data, Begin GNE device data
30+
0 /End of GNE device data, Begin Induction machine data
31+
0 /End of Induction machine data
32+
Q

‎test/benchmarks/psse/WPIDHY/ThreeBus_WPIDHY.dyr‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
101 'GENCLS' 1 0.0000 0.0000 /
2+
102 'GENROU' 1 8.0000 0.30000E-01 0.40000 0.50000E-01
3+
6.1750 0.50000E-01 1.8000 1.7000 0.30000
4+
0.55000 0.25000 0.20000 0.0000 1.0000 /
5+
102 'SEXS' 1 0.4 5.0 20.0 1.0 -50.0 50.0 /
6+
102 'WPIDHY' 1 0.10000E-01 -0.50000E-01 3.0000 0.50000
7+
3.0000 0.20000E-01 0.10000E-01 0.25000E-01 -0.25000E-01
8+
1.0000 0.20000 3.2500 52001. 0.0000
9+
0.0000 0.12000 0.54500 0.40000 0.82700
10+
0.83000 1.0000 /

‎test/results/results_eigenvalues.jl‎

Lines changed: 19 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -985,6 +985,24 @@ test59_eigvals = [
985985
]
986986

987987
test60_eigvals = [
988+
-106.13231871102334 + 0.0im
989+
-91.13024819676329 + 0.0im
990+
-63.8688081304463 + 0.0im
991+
-40.03565106938418 + 0.0im
992+
-38.75323547447537 - 8.125818822049656im
993+
-38.75323547447537 + 8.125818822049656im
994+
-6.104826779980879 + 0.0im
995+
-1.395127240875047 - 5.966607970662552im
996+
-1.395127240875047 + 5.966607970662552im
997+
-0.9449531063097173 + 0.0im
998+
-0.48179149983179287 + 0.0im
999+
-0.1941504953444676 - 0.18213498194497163im
1000+
-0.1941504953444676 + 0.18213498194497163im
1001+
-0.03057660237077197 - 0.30118043163774183im
1002+
-0.03057660237077197 + 0.30118043163774183im
1003+
]
1004+
1005+
test61_eigvals = [
9881006
-351.6684302703973 + 0.0im
9891007
-163.4858530295313 + 0.0im
9901008
-98.32270507099068 + 0.0im
@@ -995,5 +1013,5 @@ test60_eigvals = [
9951013
-2.752172823341148 + 8.057381046225618im
9961014
-1.5898057550991769 + 0.0im
9971015
-0.3644442941887857 + 0.0im
998-
-0.025119358661025975 + 0.0im
1016+
-0.025119358661025975 + 0.0im
9991017
]

‎test/results/results_initial_conditions.jl‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1807,6 +1807,36 @@ test59_x0_init = Dict(
18071807
],
18081808
)
18091809

1810+
test60_x0_init = Dict(
1811+
"V_R" => [
1812+
1.05
1813+
1.0197944718502572
1814+
0.9923907751848658
1815+
],
1816+
"V_I" => [
1817+
0.0
1818+
-0.020475233421259967
1819+
-0.1243212484458825
1820+
],
1821+
"generator-102-1" => [
1822+
0.7538967192769663
1823+
0.5555487379562241
1824+
0.7042452148052648
1825+
0.7246287886385532
1826+
0.9158416020737494
1827+
1.0
1828+
1.4986692863524897
1829+
0.044960078577949814
1830+
0.0
1831+
3.871458709170383e-12
1832+
3.871569731472846e-12
1833+
0.0
1834+
3.871569731472846e-12
1835+
0.7417441860465115
1836+
2.099999999999999
1837+
],
1838+
)
1839+
18101840
test60_x0_init = Dict(
18111841
"V_R" => [
18121842
1.05

0 commit comments

Comments
 (0)