Skip to content

Commit

Permalink
Add AVR ST6B (#388)
Browse files Browse the repository at this point in the history 
* add mass matrix entries for ST6B

* Add Comment

* Add ST6B model

* Add initialization for ST6B

* add mass matrix entries for ST6B

* Add Test 59

* Fix Typo

---------

Co-authored-by: Machado <smachado@nrel.gov>
  • Loading branch information
@SebastianManriqueM
SebastianManriqueM and Machado authored Sep 5, 2024
1 parent 29ae384 commit 0c375d2
Show file tree
Hide file tree
Showing 8 changed files with 470 additions and 0 deletions.
81 changes: 81 additions & 0 deletions src/initialization/generator_components/init_avr.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -571,3 +571,84 @@ function initialize_avr!(
avr_states[5] = -Kf / Tf * Vf0
return
end

function initialize_avr!(
device_states,
static::PSY.StaticInjection,
dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, PSY.ST6B, TG, P}},
inner_vars::AbstractVector,
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}
#Obtain Vf0 solved from Machine
Vf0 = inner_vars[Vf_var]
#Obtain measured terminal voltage
Vt0 = sqrt(inner_vars[VR_gen_var]^2 + inner_vars[VI_gen_var]^2)
#Obtain field winding current
Ifd0 = inner_vars[Xad_Ifd_var]

#Get parameters
avr = PSY.get_avr(dynamic_device)
Tr = PSY.get_Tr(avr)
K_pa = PSY.get_K_pa(avr) #k_pa>0
K_ia = PSY.get_K_ia(avr)
K_da = PSY.get_K_da(avr)
T_da = PSY.get_T_da(avr)
Va_min, Va_max = PSY.get_Va_lim(avr)
K_ff = PSY.get_K_ff(avr)
K_m = PSY.get_K_m(avr)
K_ci = PSY.get_K_ci(avr) #K_cl in pss
K_lr = PSY.get_K_lr(avr)
I_lr = PSY.get_I_lr(avr)
Vr_min, Vr_max = PSY.get_Vr_lim(avr)
Kg = PSY.get_Kg(avr)
Tg = PSY.get_Tg(avr) #T_g>0

#Compute auxiliary parameters
Vr = Vf0 / Vt0
Vg = Kg * Vf0

Vr2 = max(Vr_min, ((I_lr * K_ci) - Ifd0) * K_lr)
if (Vr - Vr2 > eps()) # or using STRICT_NLSOLVE_F_TOLERANCE?

#Check saturation

Va = (Vr + K_m * Vg) / (K_ff + K_m)

#Check ss error according to control parameters
if K_ia < eps()
Vref0 = Va / K_pa + Vt0
error(
"AVR in $(PSY.get_name(dynamic_device)) has non positive integrator gain. Please update it.",
)
else
Vref0 = Vt0
end

if !((Vr < Vr_max) && (Vr > Vr_min))
@error(
"AVR controller in $(PSY.get_name(dynamic_device)) is saturated. Consider updating the operating point."
)
end
else
error(
"Current limiter of AVR in $(PSY.get_name(dynamic_device)) is activated. Consider updating the operating point.",
)
end

PSY.set_V_ref!(avr, Vref0)
set_V_ref(dynamic_device, Vref0)
@warn(
"I_LR parameter was updated from $(I_lr) to $(Ifd0) of $(PSY.get_name(dynamic_device)) to have zero current field error."
)
PSY.set_I_lr!(avr, Ifd0 * 10.0)

#States of ST6B [:Vm, :x_i, :x_d, :Vg]
#Update AVR states
avr_ix = get_local_state_ix(dynamic_device, PSY.ST6B)
avr_states = @view device_states[avr_ix]
avr_states[1] = Vt0
avr_states[2] = Va / K_ia
avr_states[3] = 0
avr_states[4] = Vg

return
end
90 changes: 90 additions & 0 deletions src/models/generator_models/avr_models.jl
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,7 @@
##################################
###### Mass Matrix Entries #######
##################################

function mass_matrix_avr_entries!(
mass_matrix,
avr::T,
Expand Down Expand Up @@ -61,6 +65,20 @@ function mass_matrix_avr_entries!(
return
end

function mass_matrix_avr_entries!(
mass_matrix,
avr::PSY.ST6B,
global_index::Base.ImmutableDict{Symbol, Int64},
)
mass_matrix[global_index[:Vm], global_index[:Vm]] = PSY.get_Tr(avr)
mass_matrix[global_index[:x_d], global_index[:x_d]] = PSY.get_T_da(avr)
return
end

##################################
##### Differential Equations #####
##################################

function mdl_avr_ode!(
::AbstractArray{<:ACCEPTED_REAL_TYPES},
::AbstractArray{<:ACCEPTED_REAL_TYPES},
Expand Down Expand Up @@ -634,3 +652,75 @@ function mdl_avr_ode!(
inner_vars[Vf_var] = Vf
return
end

######################################################################
function mdl_avr_ode!(
device_states::AbstractArray,
output_ode::AbstractArray,
inner_vars::AbstractArray,
dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, PSY.ST6B, TG, P}},
h,
t,
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}

#Obtain references
V_ref = get_V_ref(dynamic_device)

#Obtain avr
avr = PSY.get_avr(dynamic_device)

#Obtain indices for component w/r to device
local_ix = get_local_state_ix(dynamic_device, typeof(avr))

#Define inner states for component
internal_states = @view device_states[local_ix]
Vm = internal_states[1]
x_i = internal_states[2]
x_d = internal_states[3]
Vg = internal_states[4]

#Define external states for device
V_th = sqrt(inner_vars[VR_gen_var]^2 + inner_vars[VI_gen_var]^2) # machine's terminal voltage
Vs = inner_vars[V_pss_var] # PSS output
Xad_Ifd = inner_vars[Xad_Ifd_var] # machine's field current in exciter base

#Get parameters
Tr = PSY.get_Tr(avr)
K_pa = PSY.get_K_pa(avr) #k_pa>0
K_ia = PSY.get_K_ia(avr)
K_da = PSY.get_K_da(avr)
T_da = PSY.get_T_da(avr)
Va_min, Va_max = PSY.get_Va_lim(avr)
K_ff = PSY.get_K_ff(avr)
K_m = PSY.get_K_m(avr)
K_ci = PSY.get_K_ci(avr) #K_cl in pss
K_lr = PSY.get_K_lr(avr)
I_lr = PSY.get_I_lr(avr)
Vr_min, Vr_max = PSY.get_Vr_lim(avr)
Kg = PSY.get_Kg(avr)
Tg = PSY.get_Tg(avr) #T_g>0

#Compute block derivatives
_, dVm_dt = low_pass_mass_matrix(V_th, Vm, 1.0, Tr)
pid_input = V_ref + Vs - Vm
pi_out, dx_i = pi_block_nonwindup(pid_input, x_i, K_pa, K_ia, Va_min, Va_max)
pd_out, dx_d = high_pass_mass_matrix(pid_input, x_d, K_da, T_da)
Va = pi_out + pd_out

ff_out = ((Va - Vg) * K_m) + (K_ff * Va)
V_r1 = max(((I_lr * K_ci) - Xad_Ifd) * K_lr, Vr_min)
V_r2 = clamp(ff_out, Vr_min, Vr_max)
V_r = min(V_r1, V_r2)
E_fd = V_r * Vm
_, dVg = low_pass(E_fd, Vg, Kg, Tg)

#Compute 4 States AVR ODE:
output_ode[local_ix[1]] = dVm_dt
output_ode[local_ix[2]] = dx_i
output_ode[local_ix[3]] = dx_d
output_ode[local_ix[4]] = dVg

#Update inner_vars
inner_vars[Vf_var] = E_fd
return
end
63 changes: 63 additions & 0 deletions src/post_processing/post_proc_generator.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -831,6 +831,69 @@ function _field_voltage(
return ts, Vf
end

"""
Function to obtain the field voltage time series of a Dynamic Generator with avr ESST1A.
"""
function _field_voltage(
avr::PSY.ST6B,
name::String,
res::SimulationResults,
dt::Union{Nothing, Float64, Vector{Float64}},
)
#ASSUMPTION THAT PSS IS NOT ACTIVATE
# Obtain state Vm
ts, Vm = post_proc_state_series(res, (name, :Vm), dt)

# Obtain state x_i
ts, x_i = post_proc_state_series(res, (name, :x_i), dt)

# Obtain state x_d
ts, x_d = post_proc_state_series(res, (name, :x_d), dt)

# Obtain state x_d
ts, Vg = post_proc_state_series(res, (name, :Vg), dt)

# Obtain state Xad_Ifd
ts, Xad_Ifd = post_proc_field_current_series(res, name, dt)

# Obtain PSS output
_, Vs = post_proc_pss_output_series(res, name, dt)

V_ref = PSY.get_V_ref(avr)

#Get parameters
Tr = PSY.get_Tr(avr)
K_pa = PSY.get_K_pa(avr) #k_pa>0
K_ia = PSY.get_K_ia(avr)
K_da = PSY.get_K_da(avr)
T_da = PSY.get_T_da(avr)
Va_min, Va_max = PSY.get_Va_lim(avr)
K_ff = PSY.get_K_ff(avr)
K_m = PSY.get_K_m(avr)
K_ci = PSY.get_K_ci(avr) #K_cl in pss
K_lr = PSY.get_K_lr(avr)
I_lr = PSY.get_I_lr(avr)
Vr_min, Vr_max = PSY.get_Vr_lim(avr)
Kg = PSY.get_Kg(avr)
Tg = PSY.get_Tg(avr) #T_g>0

Efd = zeros(length(ts))
for (ix, t) in enumerate(ts)
pid_input = V_ref + Vs[ix] - Vm[ix]
pi_out, dx_i = pi_block_nonwindup(pid_input, x_i[ix], K_pa, K_ia, Va_min, Va_max)
pd_out, dx_d = high_pass_mass_matrix(pid_input, x_d[ix], K_da, T_da)
Va = pi_out + pd_out
ff_out = ((Va - Vg[ix]) * K_m) + (K_ff * Va)
V_r1 = max(((I_lr * K_ci) - Xad_Ifd[ix]) * K_lr, Vr_min)
V_r2 = clamp(ff_out, Vr_min, Vr_max)
V_r = min(V_r1, V_r2)
Efd[ix] = V_r * Vm[ix]
end

return ts, Efd
end

"""
Function to obtain the pss output time series of a Dynamic Generator with pss PSSFixed.
Expand Down
32 changes: 32 additions & 0 deletions test/benchmarks/psse/ST6B/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, 250.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 ', 153.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 ', 100.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/ST6B/ThreeBus_ST6B.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.10000 0.80000 /
102 'ST6B' 1 0
0.12000E-01 50.038 45.094
0.0000 0.0000 4.8100
-3.8500 1.0000 1.0000
1.0577 17.330 2.1500
4.8100 -3.8500 1.0000 0.20000E-01 /
12 changes: 12 additions & 0 deletions test/results/results_eigenvalues.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -971,3 +971,15 @@ test58_eigvals_GateFlag = [
-0.19849913626541832 + 0.17939599481963514im
-0.04238523317744444 + 0.0im
]

test59_eigvals = [
-99.69916689593907 + 0.0im
-84.05122388240402 + 0.0im
-39.79684261365684 + 0.0im
-38.3444969887499 + 0.0im
-8.550153616541023 + 0.0im
-0.6772738854081002 - 8.637736379730404im
-0.6772738854081002 + 8.637736379730404im
-0.36823401395037775 - 0.5649228138733879im
-0.36823401395037775 + 0.5649228138733879im
]
25 changes: 25 additions & 0 deletions test/results/results_initial_conditions.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1781,3 +1781,28 @@ test58_x0_init_GateFlag = Dict(
2.099999999999999
],
)

test59_x0_init = Dict(
"V_R" => [
1.05
1.0198615749401696
0.9817964461177487
],
"V_I" => [
0.0
-0.01680380791834505
-0.15145839832533242
],
"generator-102-1" => [
0.8084146853601003
0.5302614530689579
0.7288775579091268
0.7311890654084042
0.9615873122791483
1.0
1.0199999999999998
0.047279162537359226
0.0
2.153115531256307
],
)
Loading

0 comments on commit 0c375d2

Please sign in to comment.