diff --git a/Project.toml b/Project.toml index e1d4560..404eb3b 100644 --- a/Project.toml +++ b/Project.toml @@ -11,7 +11,6 @@ NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" - [compat] DiffEqBase = "6" ForwardDiff = "~v0.10" diff --git a/src/LITS.jl b/src/LITS.jl index a281f03..00fe8da 100644 --- a/src/LITS.jl +++ b/src/LITS.jl @@ -4,7 +4,7 @@ module LITS ####################################### Structs Exports #################################### # Base Exports -export DynBranch +export DynamicLine export Simulation export run_simulation! export ThreePhaseFault @@ -23,10 +23,11 @@ import InfrastructureSystems import DiffEqBase import ForwardDiff import SparseArrays: SparseMatrixCSC -import LinearAlgebra: BLAS +#import LinearAlgebra: BLAS import LinearAlgebra: eigen import Base.to_index import NLsolve +import Base.ImmutableDict import PowerSystems const PSY = PowerSystems const IS = InfrastructureSystems @@ -41,7 +42,7 @@ include("base/simulation.jl") #Common Models include("models/branch.jl") include("models/device.jl") -include("models/kcl.jl") +include("models/kirchoff_laws.jl") include("models/dynline_model.jl") include("models/ref_transformations.jl") @@ -69,6 +70,7 @@ include("models/system.jl") #Utils include("utils/plot_utils.jl") +include("utils/immutable_dicts.jl") include("utils/print.jl") include("utils/kwargs_check.jl") diff --git a/src/base/definitions.jl b/src/base/definitions.jl index cc82974..12c659f 100644 --- a/src/base/definitions.jl +++ b/src/base/definitions.jl @@ -84,5 +84,8 @@ const INNER_VARS = "inner_vars" const YBUS = "Ybus" const CONTROL_REFS = "control_refs" const GLOBAL_VARS = "global_vars" +const VOLTAGE_BUSES_NO = "voltage_buses_no" +const CURRENT_BUSES_NO = "current_buses_no" +const TOTAL_SHUNTS = "total_shunts" const SIMULATION_ACCEPTED_KWARGS = [:initial_guess, :initialize_simulation] diff --git a/src/base/simulation.jl b/src/base/simulation.jl index 56f8cda..eafb8bd 100644 --- a/src/base/simulation.jl +++ b/src/base/simulation.jl @@ -210,18 +210,58 @@ function _index_dynamic_system!(sys::PSY.System) n_buses = length(PSY.get_components(PSY.Bus, sys)) DAE_vector = collect(falses(n_buses * 2)) global_state_index = Dict{String, Dict{Symbol, Int64}}() - n_buses = length(PSY.get_components(PSY.Bus, sys)) state_space_ix = [n_buses * 2] + current_buses_no = collect(1:n_buses) + static_bus_var_count = 2 * length(current_buses_no) + voltage_buses_no = Vector{Int}() total_states = 0 first_dyn_branch_point = -1 branches_n_states = 0 - static_bus_vars = 2 * n_buses global_vars = Dict{Symbol, Number}( :ω_sys => 1.0, :ω_sys_index => -1, #To define 0 if infinite source, bus_number otherwise, ) + total_shunts = Dict{Int, Float64}() found_ref_bus = false + dyn_branches = PSY.get_components(DynamicLine, sys) + if !(isempty(dyn_branches)) + first_dyn_branch_point = state_space_ix[1] + 1 + for br in dyn_branches + arc = PSY.get_arc(br) + n_states = PSY.get_n_states(br) + from_bus_number = PSY.get_number(arc.from) + to_bus_number = PSY.get_number(arc.to) + merge!( + +, + total_shunts, + Dict( + from_bus_number => 1 / PSY.get_b(br).from, + to_bus_number => 1 / PSY.get_b(br).to, + ), + ) + push!(voltage_buses_no, from_bus_number, to_bus_number) + DAE_vector[from_bus_number] = DAE_vector[from_bus_number + n_buses] = true + DAE_vector[to_bus_number] = DAE_vector[to_bus_number + n_buses] = true + DAE_vector = push!(DAE_vector, collect(trues(n_states))...) + total_states += n_states + _add_states_to_global!(global_state_index, state_space_ix, br) + end + + for (ix, val) in enumerate(DAE_vector[1:n_buses]) + if val + global_state_index["V_$(ix)"] = Dict(:R => ix, :I => ix + n_buses) + total_states += 2 + static_bus_var_count -= 2 + push!(voltage_buses_no, ix) + @assert static_bus_var_count >= 0 + end + end + branches_n_states = state_space_ix[1] - n_buses * 2 + else + @debug("System doesn't contain Dynamic Branches") + end + unique!(voltage_buses_no) sources = PSY.get_components(PSY.Source, sys) for s in sources btype = PSY.get_bustype(PSY.get_bus(s)) @@ -232,7 +272,6 @@ function _index_dynamic_system!(sys::PSY.System) global_vars[:ω_sys_index] = 0 #To define 0 if infinite source, bus_number otherwise, found_ref_bus = true end - dynamic_injection = PSY.get_components(PSY.DynamicInjection, sys) isempty(dynamic_injection) && error("System doesn't contain any DynamicInjection devices") @@ -240,68 +279,45 @@ function _index_dynamic_system!(sys::PSY.System) if !(:states in fieldnames(typeof(d))) continue end - btype = PSY.get_bustype(PSY.get_bus(d)) + device_bus = PSY.get_bus(d) + btype = PSY.get_bustype(device_bus) if (btype == PSY.BusTypes.REF) && found_ref_bus throw(IS.ConflictingInputsError("The system can't have more than one source or generator in the REF Bus")) end _make_device_index!(d) device_n_states = PSY.get_n_states(d) - DAE_vector = vcat(DAE_vector, collect(trues(device_n_states))) + DAE_vector = push!(DAE_vector, collect(trues(device_n_states))...) total_states += device_n_states _add_states_to_global!(global_state_index, state_space_ix, d) btype != PSY.BusTypes.REF && continue global_vars[:ω_sys_index] = global_state_index[d.name][:ω] #To define 0 if infinite source, bus_number otherwise, found_ref_bus = true end - injection_n_states = state_space_ix[1] - n_buses * 2 - - dyn_branches = PSY.get_components(DynamicLine, sys) - if !(isempty(dyn_branches)) - first_dyn_branch_point = state_space_ix[1] + 1 - for br in dyn_branches - arc = PSY.get_arc(br) - n_states = PSY.get_n_states(br) - from_bus_number = PSY.get_number(arc.from) - to_bus_number = PSY.get_number(arc.to) - DAE_vector[from_bus_number] = DAE_vector[from_bus_number + n_buses] = true - DAE_vector[to_bus_number] = DAE_vector[to_bus_number + n_buses] = true - DAE_vector = vcat(DAE_vector, collect(trues(n_states))) - total_states += n_states - _add_states_to_global!(global_state_index, state_space_ix, br) - end - - for (ix, val) in enumerate(DAE_vector[1:n_buses]) - if val - global_state_index["V_$(ix)"] = Dict(:R => ix, :I => ix + n_buses) - total_states += 2 - state_space_ix[1] += 2 - static_bus_vars -= 2 - @assert static_bus_vars >= 0 - end - end - branches_n_states = state_space_ix[1] - injection_n_states - n_buses * 2 - else - @debug("System doesn't contain Dynamic Branches") - end - @assert total_states == state_space_ix[1] - n_buses * 2 - + injection_n_states = state_space_ix[1] - branches_n_states - n_buses * 2 + @assert total_states == state_space_ix[1] - static_bus_var_count + @debug total_states + setdiff!(current_buses_no, voltage_buses_no) if !isempty(PSY.get_components(PSY.ACBranch, sys)) Ybus = PSY.Ybus(sys)[:, :] else Ybus = SparseMatrixCSC{Complex{Float64}, Int64}(zeros(n_buses, n_buses)) end - sys_ext = Dict{String, Any}() #I change it to be Any - counts = Dict{Symbol, Int64}( + sys_ext = Dict{String, Any}() + counts = Base.ImmutableDict( :total_states => total_states, :injection_n_states => injection_n_states, :branches_n_states => branches_n_states, - :first_dyn_injection_pointer => 2 * n_buses + 1, + :first_dyn_injection_pointer => 2 * n_buses + branches_n_states + 1, :first_dyn_branch_point => first_dyn_branch_point, - :total_variables => total_states + static_bus_vars, + :total_variables => total_states + static_bus_var_count, ) + sys_ext[LITS_COUNTS] = counts sys_ext[GLOBAL_INDEX] = global_state_index + sys_ext[VOLTAGE_BUSES_NO] = voltage_buses_no + sys_ext[CURRENT_BUSES_NO] = current_buses_no sys_ext[YBUS] = Ybus + sys_ext[TOTAL_SHUNTS] = total_shunts sys_ext[GLOBAL_VARS] = global_vars @assert sys_ext[GLOBAL_VARS][:ω_sys_index] != -1 sys.internal.ext = sys_ext @@ -318,9 +334,11 @@ get_system_state_count(sys::PSY.System) = PSY.get_ext(sys)[LITS_COUNTS][:total_s get_variable_count(sys::PSY.System) = PSY.get_ext(sys)[LITS_COUNTS][:total_variables] get_device_index(sys::PSY.System, device::D) where {D <: PSY.DynamicInjection} = PSY.get_ext(sys)[GLOBAL_INDEX][device.name] - get_inner_vars(device::PSY.DynamicInjection) = device.ext[INNER_VARS] get_ω_sys(sys::PSY.System) = PSY.get_ext(sys)[GLOBAL_VARS][:ω_sys] +get_current_bus_no(sys::PSY.System) = PSY.get_ext(sys)[CURRENT_BUSES_NO] +get_voltage_bus_no(sys::PSY.System) = PSY.get_ext(sys)[VOLTAGE_BUSES_NO] +get_total_shunts(sys::PSY.System) = PSY.get_ext(sys)[TOTAL_SHUNTS] function _get_internal_mapping( device::PSY.DynamicInjection, @@ -380,11 +398,6 @@ function small_signal_analysis(sim::Simulation; kwargs...) @error("Reset the simulation") end - dyn_branches = PSY.get_components(DynamicLine, sim.system) - if !(isempty(dyn_branches)) - @error("Small Signal Analysis is not currently supported for models with DynamicLines") - end - _change_vector_type(sim.system) var_count = LITS.get_variable_count(sim.system) dx0 = zeros(var_count) #Define a vector of zeros for the derivative @@ -398,15 +411,19 @@ function small_signal_analysis(sim::Simulation; kwargs...) out = zeros(var_count) #Define a vector of zeros for the output x_eval = get(kwargs, :operating_point, sim.x0_init) jacobian = ForwardDiff.jacobian(sysf!, out, x_eval) - first_dyn_injection_pointer = - PSY.get_ext(sim.system)[LITS_COUNTS][:first_dyn_injection_pointer] - bus_size = length(PSY.get_components(PSY.Bus, sim.system)) - alg_states = 1:(2 * bus_size) - diff_states = first_dyn_injection_pointer:var_count - fx = jacobian[diff_states, diff_states] + n_buses = length(PSY.get_components(PSY.Bus, sim.system)) + diff_states = collect(trues(var_count)) + diff_states[1:(2 * n_buses)] .= false + for b_no in get_voltage_bus_no(sim.system) + diff_states[b_no] = true + diff_states[b_no + n_buses] = true + end + alg_states = .!diff_states + fx = @view jacobian[diff_states, diff_states] gy = jacobian[alg_states, alg_states] - fy = jacobian[diff_states, alg_states] - gx = jacobian[alg_states, diff_states] + fy = @view jacobian[diff_states, alg_states] + gx = @view jacobian[alg_states, diff_states] + # TODO: Make operation using BLAS! reduced_jacobian = fx - fy * inv(gy) * gx vals, vect = eigen(reduced_jacobian) return SmallSignalOutput( diff --git a/src/models/branch.jl b/src/models/branch.jl index 12ae184..297c640 100644 --- a/src/models/branch.jl +++ b/src/models/branch.jl @@ -54,7 +54,7 @@ end function branch!( x, dx, - output_ode::Vector{Float64}, + output_ode::Vector{T}, V_r_from, V_i_from, V_r_to, @@ -64,17 +64,14 @@ function branch!( current_r_to, current_i_to, ix_range::UnitRange{Int64}, - ix_dx::Vector{Int64}, ode_range::UnitRange{Int64}, branch::DynamicLine{PSY.Line}, sys::PSY.System, -) +) where {T <: Real} #Obtain local device states n_states = PSY.get_n_states(branch) device_states = @view x[ix_range] - dv_from = view(dx, ix_dx[1:2]) - dv_to = view(dx, ix_dx[3:4]) #Obtain references Sbase = PSY.get_basepower(sys) @@ -92,8 +89,6 @@ function branch!( current_i_from, current_r_to, current_i_to, - dv_from, - dv_to, sys_f, branch, ) diff --git a/src/models/dynline_model.jl b/src/models/dynline_model.jl index 9c2d167..84bdfd2 100644 --- a/src/models/dynline_model.jl +++ b/src/models/dynline_model.jl @@ -9,8 +9,6 @@ function mdl_line_ode!( current_i_from, current_r_to, current_i_to, - dv_from, - dv_to, sys_f, branch::DynamicLine{PSY.Line}, ) @@ -18,13 +16,6 @@ function mdl_line_ode!( L = PSY.get_x(branch) R = PSY.get_r(branch) ω_b = sys_f * 2 * π - c_from = PSY.get_b(branch).from - c_to = PSY.get_b(branch).to - - current_r_from[1] -= (1.0 / ω_b) * c_from * dv_from[1] - c_from * V_i_from[1] - current_i_from[1] -= (1.0 / ω_b) * c_from * dv_from[2] + c_from * V_r_from[1] - current_r_to[1] -= (1.0 / ω_b) * c_to * dv_to[1] - c_to * V_i_to[1] - current_i_to[1] -= (1.0 / ω_b) * c_to * dv_to[2] + c_to * V_r_to[1] Il_r = device_states[1] Il_i = device_states[2] diff --git a/src/models/kcl.jl b/src/models/kcl.jl deleted file mode 100644 index bc67f8d..0000000 --- a/src/models/kcl.jl +++ /dev/null @@ -1,13 +0,0 @@ -""" - Kirchoff current law for buses. - I_gen_r[j]: Real current injection from all generators connected at bus j. - It is zero if no generator is connected at bus j. - I_gen_i[j]: Real current injection from all generators connected at bus j. - It is zero if no generator is connected at bus j. - -""" - -function kcl(Ybus, V_r, V_i, I_injections_r, I_injections_i) - I_bus = Ybus * (V_r + V_i .* 1im) - return [I_injections_r - real(I_bus); I_injections_i - imag(I_bus)] -end diff --git a/src/models/kirchoff_laws.jl b/src/models/kirchoff_laws.jl new file mode 100644 index 0000000..22c9d3d --- /dev/null +++ b/src/models/kirchoff_laws.jl @@ -0,0 +1,32 @@ +""" + Kirchoff law for buses. + I_gen_r[j]: Real current injection from all generators connected at bus j. + It is zero if no generator is connected at bus j. + I_gen_i[j]: Real current injection from all generators connected at bus j. + It is zero if no generator is connected at bus j. + +""" +function kirchoff_laws(sys, V_r, V_i, I_injections_r, I_injections_i, dx) + Ybus = PSY.get_ext(sys)[YBUS] + # TODO: Make more efficent calculation here. Note: BLAS doesn't work because the the type of Vr and Vi is not Matrix of Complex + I_bus = Ybus * (V_r + V_i .* 1im) + I_balance = [real(I_bus) - I_injections_r; imag(I_bus) - I_injections_i] + + voltage_buses = get_voltage_bus_no(sys) + isempty(voltage_buses) && return I_balance + + sys_f = PSY.get_frequency(sys) + ω_b = 2.0 * π * sys_f + n_buses = length(PSY.get_components(PSY.Bus, sys)) + shunts = get_total_shunts(sys) + for bus_no in voltage_buses + shunt_multiplier = shunts[bus_no] + I_balance[bus_no] = + -ω_b * I_balance[bus_no] * shunt_multiplier + ω_b * V_i[bus_no] - dx[bus_no] + I_balance[bus_no + n_buses] = + -ω_b * I_balance[bus_no + n_buses] * shunt_multiplier - ω_b * V_r[bus_no] - + dx[bus_no + n_buses] + end + + return I_balance +end diff --git a/src/models/system.jl b/src/models/system.jl index cb27f88..4950147 100644 --- a/src/models/system.jl +++ b/src/models/system.jl @@ -21,6 +21,7 @@ function system!(out::Vector{T}, dx, x, sys, t) where {T <: Real} V_r = @view x[1:bus_size] V_i = @view x[(bus_size + 1):bus_vars_count] Sbase = PSY.get_basepower(sys) + # TODO: Don't create this matrices here every iteration I_injections_r = zeros(T, bus_size) I_injections_i = zeros(T, bus_size) injection_ode = zeros(T, get_n_injection_states(sys)) @@ -48,6 +49,19 @@ function system!(out::Vector{T}, dx, x, sys, t) where {T <: Real} out[ix_range] = injection_ode[ode_range] - dx[ix_range] end + for d in PSY.get_components(PSY.StaticInjection, sys) + bus_n = PSY.get_number(PSY.get_bus(d)) + + device!( + view(V_r, bus_n), + view(V_i, bus_n), + view(I_injections_r, bus_n), + view(I_injections_i, bus_n), + d, + sys, + ) + end + dyn_branches = PSY.get_components(DynamicLine, sys) if !(isempty(dyn_branches)) for br in dyn_branches @@ -55,12 +69,6 @@ function system!(out::Vector{T}, dx, x, sys, t) where {T <: Real} n_states = PSY.get_n_states(br) from_bus_number = PSY.get_number(arc.from) to_bus_number = PSY.get_number(arc.to) - ix_dx = [ - from_bus_number, - from_bus_number + bus_size, - to_bus_number, - to_bus_number + bus_size, - ] ix_range = range(branches_start, length = n_states) ode_range = range(branches_count, length = n_states) branches_count = branches_count + n_states @@ -79,7 +87,6 @@ function system!(out::Vector{T}, dx, x, sys, t) where {T <: Real} view(I_injections_r, to_bus_number), view(I_injections_i, to_bus_number), ix_range, - ix_dx, ode_range, br, sys, @@ -88,19 +95,6 @@ function system!(out::Vector{T}, dx, x, sys, t) where {T <: Real} end end - for d in PSY.get_components(PSY.StaticInjection, sys) - bus_n = PSY.get_number(PSY.get_bus(d)) - - device!( - view(V_r, bus_n), - view(V_i, bus_n), - view(I_injections_r, bus_n), - view(I_injections_i, bus_n), - d, - sys, - ) - end - - out[bus_range] = kcl(PSY.get_ext(sys)[YBUS], V_r, V_i, I_injections_r, I_injections_i) + out[bus_range] = kirchoff_laws(sys, V_r, V_i, I_injections_r, I_injections_i, dx) end diff --git a/src/utils/immutable_dicts.jl b/src/utils/immutable_dicts.jl new file mode 100644 index 0000000..fa3b08f --- /dev/null +++ b/src/utils/immutable_dicts.jl @@ -0,0 +1,7 @@ +function Base.ImmutableDict(KV::Pair{K, V}, KVs::Pair{K, V}...) where {K, V} + d = Base.ImmutableDict(KV) + for p in KVs + d = Base.ImmutableDict(d, p) + end + return d +end diff --git a/test/data_tests/dynamic_test_data.jl b/test/data_tests/dynamic_test_data.jl index 1abd815..fd45a5c 100644 --- a/test/data_tests/dynamic_test_data.jl +++ b/test/data_tests/dynamic_test_data.jl @@ -550,6 +550,25 @@ function inv_case78(nodes) ) #filter end +function inv_case9(nodes) + return PSY.DynamicInverter( + 1, #Number + "DARCO", #name + nodes[3], #bus + 1.0, # ω_ref, + 0.8, #V_ref + 0.5, #P_ref + -0.3, #Q_ref + 100.0, #MVABase + converter_case78(), #converter + outer_control_test(), #outer control + vsc_test(), #inner control voltage source + dc_source_case78(), #dc source + pll_test(), #pll + filter_test(), + ) #filter +end + ###################################### ######## System Constructors ######### ###################################### diff --git a/test/test_case08_staticbranches.jl b/test/test_case08_staticbranches.jl index 35bf1d6..56b0186 100644 --- a/test/test_case08_staticbranches.jl +++ b/test/test_case08_staticbranches.jl @@ -29,7 +29,7 @@ case8_gen = dyn_gen_case8(nodes_case8) ############### Inverter Data ######################## -case8_inv = inv_case78(nodes_case8) +case8_inv = inv_case9(nodes_case8) ######################### Dynamical System ######################## diff --git a/test/test_case09_dynbranches.jl b/test/test_case09_dynbranches.jl index d52bfe6..18b5406 100644 --- a/test/test_case09_dynbranches.jl +++ b/test/test_case09_dynbranches.jl @@ -30,7 +30,7 @@ case9_gen = dyn_gen_case9(nodes_case9) ############### Inverter Data ######################## -case9_inv = inv_case78(nodes_case9) +case9_inv = inv_case9(nodes_case9) ######################### Dynamical System ########################