Consolidate Sauer Pai model and add Static Load post processing currents

src/base/simulation_results.jl

@@ -81,7 +81,11 @@ function post_proc_voltage_current_series(
     bus_ix = get(bus_lookup, PSY.get_number(PSY.get_bus(device)), -1)
     ts, V_R, V_I = post_proc_voltage_series(solution, bus_ix, n_buses, dt)
     dyn_device = PSY.get_dynamic_injector(device)
-    _, I_R, I_I = compute_output_current(res, dyn_device, V_R, V_I, dt)
+    if isnothing(dyn_device)
+        _, I_R, I_I = compute_output_current(res, device, V_R, V_I, dt)
+    else
+        _, I_R, I_I = compute_output_current(res, dyn_device, V_R, V_I, dt)
+    end
     return ts, V_R, V_I, I_R, I_I
 end
 

src/initialization/generator_components/init_machine.jl

@@ -149,6 +149,114 @@ function initialize_mach_shaft!(
     return
 end
 
+"""
+Initialitation of model of 6-state (SauerPai) synchronous machine in Julia.
+Refer to Power System Modelling and Scripting by F. Milano for the equations
+"""
+function initialize_mach_shaft!(
+    device_states,
+    static::PSY.StaticInjection,
+    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{PSY.SauerPaiMachine, S, A, TG, P}},
+    inner_vars::AbstractVector,
+) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
+    #PowerFlow Data
+
+    P0 = PSY.get_active_power(static)
+    Q0 = PSY.get_reactive_power(static)
+    Vm = PSY.get_magnitude(PSY.get_bus(static))
+    θ = PSY.get_angle(PSY.get_bus(static))
+    S0 = P0 + Q0 * 1im
+    V_R = Vm * cos(θ)
+    V_I = Vm * sin(θ)
+    V = V_R + V_I * 1im
+    I = conj(S0 / V)
+
+    #Machine Data
+    machine = PSY.get_machine(dynamic_device)
+    R = PSY.get_R(machine)
+    Xd = PSY.get_Xd(machine)
+    Xq = PSY.get_Xq(machine)
+    Xd_p = PSY.get_Xd_p(machine)
+    Xq_p = PSY.get_Xq_p(machine)
+    Xd_pp = PSY.get_Xd_pp(machine)
+    Xq_pp = PSY.get_Xq_pp(machine)
+    Xl = PSY.get_Xl(machine)
+    Td0_p = PSY.get_Td0_p(machine)
+    γ_d1 = PSY.get_γ_d1(machine)
+    γ_q1 = PSY.get_γ_q1(machine)
+    γ_d2 = PSY.get_γ_d2(machine)
+    γ_q2 = PSY.get_γ_q2(machine)
+
+    #States of SauerPaiMachine are [1] ψq, [2] ψd, [3] eq_p, [4] ed_p, [5] ψd_pp and [6] ψq_pp
+    δ0 = angle(V + (R + Xq * 1im) * I)
+    ω0 = 1.0
+    τm0 = real(V * conj(I))
+    @assert isapprox(τm0, P0; atol = STRICT_NLSOLVE_F_TOLERANCE)
+    #To solve: δ, τm, Vf0, eq_p, ed_p
+    function f!(out, x)
+        δ = x[1]
+        τm = x[2]
+        Vf0 = x[3]
+        ψq = x[4]
+        ψd = x[5]
+        eq_p = x[6]
+        ed_p = x[7]
+        ψd_pp = x[8]
+        ψq_pp = x[9]
+
+        V_dq = ri_dq(δ) * [V_R; V_I]
+        i_d = (1.0 / Xd_pp) * (γ_d1 * eq_p - ψd + (1 - γ_d1) * ψd_pp)      #15.15
+        i_q = (1.0 / Xq_pp) * (-γ_q1 * ed_p - ψq + (1 - γ_q1) * ψq_pp)     #15.15
+        τ_e = ψd * i_q - ψq * i_d               #15.6
+        out[1] = τm - τ_e #Mechanical Torque
+        out[2] = P0 - (V_dq[1] * i_d + V_dq[2] * i_q) #Output Power
+        out[3] = Q0 - (V_dq[2] * i_d - V_dq[1] * i_q) #Output Reactive Power
+        out[4] = R * i_q - ω0 * ψd + V_dq[2]                                    #15.9 ψq
+        out[5] = R * i_d + ω0 * ψq + V_dq[1]                                    #15.9 ψd
+        out[6] =
+            -eq_p - (Xd - Xd_p) * (i_d - γ_d2 * ψd_pp - (1 - γ_d1) * i_d + γ_d2 * eq_p) +
+            Vf0    #15.13
+        out[7] = -ed_p + (Xq - Xq_p) * (i_q - γ_q2 * ψq_pp - (1 - γ_q1) * i_q - γ_d2 * ed_p)          #15.13
+        out[8] = -ψd_pp + eq_p - (Xd_p - Xl) * i_d      #15.13
+        out[9] = -ψq_pp - ed_p - (Xq_p - Xl) * i_q #15.13
+    end
+
+    V_dq0 = ri_dq(δ0) * [V_R; V_I]
+    x0 = [δ0, τm0, 1.0, V_dq0[1], V_dq0[2], V_dq0[2], V_dq0[1], V_dq0[2], V_dq0[1]]
+    sol = NLsolve.nlsolve(f!, x0, ftol = STRICT_NLSOLVE_F_TOLERANCE)
+    if !NLsolve.converged(sol)
+        @warn("Initialization in Machine $(PSY.get_name(static)) failed")
+    else
+        sol_x0 = sol.zero
+        #Update terminal voltages
+        inner_vars[VR_gen_var] = V_R
+        inner_vars[VI_gen_var] = V_I
+        #Update δ and ω of Shaft. Works for every Shaft.
+        shaft_ix = get_local_state_ix(dynamic_device, S)
+        shaft_states = @view device_states[shaft_ix]
+        shaft_states[1] = sol_x0[1] #δ
+        shaft_states[2] = ω0 #ω
+        #Update Mechanical and Electrical Torque on Generator
+        inner_vars[τe_var] = sol_x0[2]
+        inner_vars[τm_var] = sol_x0[2]
+        #Update Vf for AVR in OneDOneQ Machine.
+        inner_vars[Vf_var] = sol_x0[3]
+        #Update states for Machine
+        machine_ix = get_local_state_ix(dynamic_device, PSY.SauerPaiMachine)
+        machine_states = @view device_states[machine_ix]
+        machine_states[1] = sol_x0[4] #ψq
+        machine_states[2] = sol_x0[5] #ψd
+        machine_states[3] = sol_x0[6] #eq_p
+        machine_states[4] = sol_x0[7] #ed_p
+        machine_states[5] = sol_x0[8] #eq_pp
+        machine_states[6] = sol_x0[9] #ed_pp
+        #Update fluxes inner vars
+        inner_vars[ψq_var] = sol_x0[4]
+        inner_vars[ψd_var] = sol_x0[5]
+    end
+    return
+end
+
 """
 Initialitation of model of 6-state (Marconato) synchronous machine in Julia.
 Refer to Power System Modelling and Scripting by F. Milano for the equations

src/models/generator_models/machine_models.jl

@@ -122,6 +122,96 @@ function mdl_machine_ode!(
     return
 end
 
+"""
+Model of 6-state (SauerPaiMachine) synchronous machine in Julia.
+Refer to Power System Modelling and Scripting by F. Milano for the equations
+"""
+function mdl_machine_ode!(
+    device_states::AbstractArray{<:ACCEPTED_REAL_TYPES},
+    output_ode::AbstractArray{<:ACCEPTED_REAL_TYPES},
+    inner_vars::AbstractArray{<:ACCEPTED_REAL_TYPES},
+    current_r::AbstractArray{<:ACCEPTED_REAL_TYPES},
+    current_i::AbstractArray{<:ACCEPTED_REAL_TYPES},
+    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{PSY.SauerPaiMachine, S, A, TG, P}},
+) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
+    Sbase = get_system_base_power(dynamic_device)
+    f0 = get_system_base_frequency(dynamic_device)
+
+    #Obtain indices for component w/r to device
+    local_ix = get_local_state_ix(dynamic_device, PSY.SauerPaiMachine)
+
+    #Define internal states for component
+    internal_states = @view device_states[local_ix]
+    ψq = internal_states[1]
+    ψd = internal_states[2]
+    eq_p = internal_states[3]
+    ed_p = internal_states[4]
+    ψd_pp = internal_states[5]
+    ψq_pp = internal_states[6]
+
+    #Obtain external states inputs for component
+    external_ix = get_input_port_ix(dynamic_device, PSY.SauerPaiMachine)
+    δ = device_states[external_ix[1]]
+    ω = device_states[external_ix[2]]
+
+    #Obtain inner variables for component
+    V_tR = inner_vars[VR_gen_var]
+    V_tI = inner_vars[VI_gen_var]
+    Vf = inner_vars[Vf_var]
+
+    #Get parameters
+    machine = PSY.get_machine(dynamic_device)
+    R = PSY.get_R(machine)
+    Xd = PSY.get_Xd(machine)
+    Xq = PSY.get_Xq(machine)
+    Xd_p = PSY.get_Xd_p(machine)
+    Xq_p = PSY.get_Xq_p(machine)
+    Xd_pp = PSY.get_Xd_pp(machine)
+    Xq_pp = PSY.get_Xq_pp(machine)
+    Xl = PSY.get_Xl(machine)
+    Td0_p = PSY.get_Td0_p(machine)
+    Tq0_p = PSY.get_Tq0_p(machine)
+    Td0_pp = PSY.get_Td0_pp(machine)
+    Tq0_pp = PSY.get_Tq0_pp(machine)
+    γ_d1 = PSY.get_γ_d1(machine)
+    γ_q1 = PSY.get_γ_q1(machine)
+    γ_d2 = PSY.get_γ_d2(machine)
+    γ_q2 = PSY.get_γ_q2(machine)
+    basepower = PSY.get_base_power(dynamic_device)
+    ``
+    #RI to dq transformation
+    V_dq = ri_dq(δ) * [V_tR; V_tI]
+
+    #Obtain electric variables
+    i_d = (1.0 / Xd_pp) * (γ_d1 * eq_p - ψd + (1 - γ_d1) * ψd_pp)      #15.15
+    i_q = (1.0 / Xq_pp) * (-γ_q1 * ed_p - ψq + (1 - γ_q1) * ψq_pp)     #15.15
+    τ_e = ψd * i_q - ψq * i_d               #15.6
+
+    #Compute ODEs
+    output_ode[local_ix[1]] = 2 * π * f0 * (R * i_q - ω * ψd + V_dq[2])                        #15.9 ψq
+    output_ode[local_ix[2]] = 2 * π * f0 * (R * i_d + ω * ψq + V_dq[1])                        #15.9 ψd
+    output_ode[local_ix[3]] =
+        (1.0 / Td0_p) *
+        (-eq_p - (Xd - Xd_p) * (i_d - γ_d2 * ψd_pp - (1 - γ_d1) * i_d + γ_d2 * eq_p) + Vf)     #15.13 eq_p
+    output_ode[local_ix[4]] =
+        (1.0 / Tq0_p) *
+        (-ed_p + (Xq - Xq_p) * (i_q - γ_q2 * ψq_pp - (1 - γ_q1) * i_q - γ_d2 * ed_p))          #15.13 ed_p
+    output_ode[local_ix[5]] = (1.0 / Td0_pp) * (-ψd_pp + eq_p - (Xd_p - Xl) * i_d)        #15.13 ψd_pp
+    output_ode[local_ix[6]] = (1.0 / Tq0_pp) * (-ψq_pp - ed_p - (Xq_p - Xl) * i_q)        #15.13 ψq_pp
+
+    #Update inner_vars
+    inner_vars[τe_var] = τ_e
+
+    #Compute current from the generator to the grid
+    I_RI = (basepower / Sbase) * dq_ri(δ) * [i_d; i_q]
+
+    #Update current
+    current_r[1] += I_RI[1]
+    current_i[1] += I_RI[2]
+
+    return
+end
+
 """
 Model of 6-state (MarconatoMachine) synchronous machine in Julia.
 Refer to Power System Modelling and Scripting by F. Milano for the equations

src/post_processing/post_proc_generator.jl

@@ -235,6 +235,48 @@ function _machine_current(
     return ts, I_R, I_I
 end
 
+"""
+Function to obtain the output current time series of a SauerPaiMachine model out of the DAE Solution. It is dispatched via the machine type.
+"""
+function _machine_current(
+    machine::PSY.SauerPaiMachine,
+    name::String,
+    ::Vector{Float64},
+    ::Vector{Float64},
+    base_power_ratio::Float64,
+    res::SimulationResults,
+    dt::Union{Nothing, Float64},
+)
+    ts, δ = post_proc_state_series(res, (name, :δ), dt)
+    _, eq_p = post_proc_state_series(res, (name, :eq_p), dt)
+    _, ed_p = post_proc_state_series(res, (name, :ed_p), dt)
+    _, ψd = post_proc_state_series(res, (name, :ψd), dt)
+    _, ψq = post_proc_state_series(res, (name, :ψq), dt)
+    _, ψd_pp = post_proc_state_series(res, (name, :ψd_pp), dt)
+    _, ψq_pp = post_proc_state_series(res, (name, :ψq_pp), dt)
+
+    #Get parameters
+    Xd_pp = PSY.get_Xd_pp(machine)
+    Xq_pp = PSY.get_Xq_pp(machine)
+    γ_d1 = PSY.get_γ_d1(machine)
+    γ_q1 = PSY.get_γ_q1(machine)
+
+    i_dq = Vector{Float64}(undef, 2)
+    I_R = similar(δ, Float64)
+    I_I = similar(δ, Float64)
+
+    for ix in 1:length(δ)
+        v = δ[ix]
+
+        #Obtain electric current
+        i_dq[1] = (1.0 / Xd_pp) * (γ_d1 * eq_p[ix] - ψd[ix] + (1 - γ_d1) * ψd_pp[ix])      #15.15
+        i_dq[2] = (1.0 / Xq_pp) * (-γ_q1 * ed_p[ix] - ψq[ix] + (1 - γ_q1) * ψq_pp[ix])     #15.15
+
+        I_R[ix], I_I[ix] = base_power_ratio * dq_ri(v) * i_dq
+    end
+    return ts, I_R, I_I
+end
+
 """
 Function to obtain the output current time series of a GENROU/GENROE model out of the DAE Solution. It is dispatched via the machine type.
 """

src/post_processing/post_proc_loads.jl

@@ -102,3 +102,85 @@ function compute_output_current(
 
     return ts, I_R, I_I
 end
+
+"""
+Function to obtain the output current time series of a PowerLoad model. 
+
+"""
+function compute_output_current(
+    res::SimulationResults,
+    device::PSY.PowerLoad,
+    V_R::Vector{Float64},
+    V_I::Vector{Float64},
+    dt::Union{Nothing, Float64},
+)
+    #TODO: We should dispatch this using the ZipLoad model that we have, but that would
+    #      require to properly have access to it in the SimResults.
+    #TODO: Load is assumed to be connected. We need proper ways of keep tracking when
+    #      something is disconnected
+    solution = res.solution
+    if dt === nothing
+        ix_t = unique(i -> solution.t[i], eachindex(solution.t))
+        ts = solution.t[ix_t]
+    else
+        ts = range(0, stop = solution.t[end], step = dt)
+    end
+
+    V0 = sqrt(V_R[1]^2 + V_I[1]^2)
+    V_mag = sqrt.(V_R .^ 2 + V_I .^ 2)
+    P = PSY.get_active_power(device)
+    Q = PSY.get_reactive_power(device)
+    I_R = similar(V_mag)
+    I_I = similar(V_mag)
+
+    if PSY.get_model(device) == PSY.LoadModels.ConstantImpedance
+        I_R = (1.0 / V0)^2 .* (P .* V_R + Q .* V_I)
+        I_I = (1.0 / V0)^2 .* (P .* V_I - Q .* V_R)
+    elseif PSY.get_model(device) == PSY.LoadModels.ConstantCurrent
+        I_R = (1.0 / V0) .* (P .* V_R + Q .* V_I) ./ V_mag
+        I_I = (1.0 / V0) .* (P .* V_I - Q .* V_R) ./ V_mag
+    elseif PSY.get_model(device) == PSY.LoadModels.ConstantPower
+        I_R = (P .* V_R + Q .* V_I) ./ V_mag .^ 2
+        I_I = (P .* V_I - Q .* V_R) ./ V_mag .^ 2
+    else
+        @error("Load Model not supported. Returning zeros")
+    end
+    return ts, I_R, I_I
+end
+
+"""
+Function to obtain the output current time series of a ExponentialLoad model. 
+
+"""
+function compute_output_current(
+    res::SimulationResults,
+    device::PSY.ExponentialLoad,
+    V_R::Vector{Float64},
+    V_I::Vector{Float64},
+    dt::Union{Nothing, Float64},
+)
+    #TODO: We should dispatch this using the ZipLoad model that we have, but that would
+    #      require to properly have access to it in the SimResults.
+    #TODO: Load is assumed to be connected. We need proper ways of keep tracking when
+    #      something is disconnected
+    solution = res.solution
+    if dt === nothing
+        ix_t = unique(i -> solution.t[i], eachindex(solution.t))
+        ts = solution.t[ix_t]
+    else
+        ts = range(0, stop = solution.t[end], step = dt)
+    end
+
+    V0 = sqrt(V_R[1]^2 + V_I[1]^2)
+    V_mag = sqrt.(V_R .^ 2 + V_I .^ 2)
+    P = PSY.get_active_power(device)
+    Q = PSY.get_reactive_power(device)
+    α = PSY.get_active_power_coefficient(device)
+    β = PSY.get_reactive_power_coefficient(device)
+
+    I_R =
+        P .* V_R .* (V_mag .^ (α - 2.0) ./ V0^α) + Q .* V_I .* (V_mag .^ (β - 2.0) ./ V0^β)
+    I_I =
+        P .* V_I .* (V_mag .^ (α - 2.0) ./ V0^α) - Q .* V_R .* (V_mag .^ (β - 2.0) ./ V0^β)
+    return ts, I_R, I_I
+end