Merge pull request #40 from Energy-MAC/rh/update_tg

Add multimachine test case

README.md

@@ -13,4 +13,6 @@ To install and develop run the following:
 julia> ] add LITS
 ```
 
-The use examples can be found in the [LITS-Examples](https://github.com/Energy-MAC/LITS-Examples) repository
+The use examples can be found in the [LITS-Examples](https://github.com/Energy-MAC/LITS-Examples) repository.
+
+The pre-print paper can be found in our [arXiv preprint](https://arxiv.org/abs/2003.02957).

src/base/simulation.jl

@@ -23,15 +23,14 @@ function Simulation(
     DAE_vector = _index_dynamic_system!(system)
     var_count = get_variable_count(system)
 
+    flat_start = zeros(var_count)
+    bus_count = length(PSY.get_components(PSY.Bus, system))
+    flat_start[1:bus_count] .= 1.0
+    x0_init = get(kwargs, :initial_guess, flat_start)
+
     if initialize_simulation
         @info("Initializing Simulation States")
-        flat_start = zeros(var_count)
-        bus_count = length(PSY.get_components(PSY.Bus, system))
-        flat_start[1:bus_count] .= 1.0
-        x0_guess = get(kwargs, :initial_guess, flat_start)
-        x0_init, initialized = _calculate_initial_conditions(system, x0_guess)
-    else
-        x0_init = zeros(var_count)
+        initialized = _calculate_initial_conditions!(system, x0_init)
     end
 
     dx0 = zeros(var_count)
@@ -92,7 +91,7 @@ function _build_perturbations(perturbations::Vector{<:Perturbation})
     return callback_set, tstops
 end
 
-function _calculate_initial_conditions(sys::PSY.System, initial_guess::Vector{Float64})
+function _calculate_initial_conditions!(sys::PSY.System, initial_guess::Vector{Float64})
     # TODO: Code to refine initial_guess
     var_count = get_variable_count(sys)
     dx0 = zeros(var_count) #Define a vector of zeros for the derivative
@@ -113,8 +112,8 @@ function _calculate_initial_conditions(sys::PSY.System, initial_guess::Vector{Fl
     if !NLsolve.converged(sys_solve)
         @warn("Initialization failed, initial conditions do not meet conditions for an stable equilibrium")
     end
-
-    return sys_solve.zero, NLsolve.converged(sys_solve)
+    initial_guess .= sys_solve.zero
+    return NLsolve.converged(sys_solve)
 end
 
 function _index_local_states!(
@@ -251,7 +250,7 @@ function _index_dynamic_system!(sys::PSY.System)
         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][:ω] #To define 0 if infinite source, bus_number otherwise,
+        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

src/models/generator_models/tg_models.jl

@@ -97,7 +97,7 @@ function mdl_tg_ode!(
     τ_max = PSY.get_τ_max(tg)
 
     #Compute auxiliary parameters
-    τ_m = xg + (1.0 / R) * (ω_ref - ω[1]) + P_ref / 1.0
+    τ_m = (1.0 / R) * (T1 / T2) * (ω_ref - ω[1]) + P_ref / 1.0 + xg
 
     #Set anti-windup for τ_m. NOT WORKING
     #if τ_m > τ_max

test/data_tests/dynamic_test_data.jl

@@ -105,6 +105,20 @@ machine_full_kundur() = PSY.FullMachine(
     555.0,
 ) #MVABase
 
+machine_multi_ref() = BaseMachine(
+    0.0, #R
+    0.2995, #Xd_p
+    1.0901, #eq_p
+    100.0,
+)  #MVABase
+
+machine_multi() = BaseMachine(
+    0.0, #R
+    0.2995, #Xd_p
+    0.9516, #eq_p
+    100.0,
+)  #MVABase
+
 ######## Shaft Data #########
 
 shaft_damping() = PSY.SingleMass(
@@ -159,8 +173,8 @@ tg_type1() = PSY.TGTypeI(
 
 tg_type2() = PSY.TGTypeII(
     0.05, #R
-    1.0, #T1
-    2.0, #T2
+    2.0, #T1
+    1.0, #T2
     1.5, #τ_max
     0.1,
 ) #τ_min
@@ -390,6 +404,40 @@ function dyn_gen_case9(nodes)
     ) #pss
 end
 
+function dyn_gen_case10_ref(nodes)
+    return PSY.DynamicGenerator(
+        1, #Number
+        "Case10Gen1",
+        nodes[1], #bus
+        1.0, # ω_ref,
+        1.0, #V_ref
+        0.41, #P_ref
+        0.0, #Q_ref
+        machine_multi_ref(), #machine
+        shaft_damping(), #shaft
+        avr_none(), #avr
+        tg_none(), #tg
+        pss_none(),
+    ) #pss
+end
+
+function dyn_gen_case10_2(nodes)
+    return PSY.DynamicGenerator(
+        1, #Number
+        "Case10Gen2",
+        nodes[3], #bus
+        1.0, # ω_ref,
+        1.0, #V_ref
+        0.4, #P_ref
+        0.0, #Q_ref
+        machine_multi(), #machine
+        shaft_damping(), #shaft
+        avr_none(), #avr
+        tg_type2(), #tg
+        pss_none(),
+    ) #pss
+end
+
 ######################################
 ############# Inverters ##############
 ######################################
@@ -537,6 +585,32 @@ function system_no_inv(nodes, branches, loads, sources, gens)
     return sys
 end
 
+function system_no_inv_no_sources(nodes, branches, loads, gens)
+    #Create system with BasePower = 100 MVA and nominal frequency 60 Hz.
+    sys = PSY.System(100.0, frequency = 60.0)
+
+    #Add buses
+    for bus in nodes
+        PSY.add_component!(sys, bus)
+    end
+
+    #Add lines
+    for lines in branches
+        PSY.add_component!(sys, lines)
+    end
+
+    #Add loads
+    for load in loads
+        PSY.add_component!(sys, load)
+    end
+
+    #Add generator
+    for gen in gens
+        PSY.add_component!(sys, gen)
+    end
+    return sys
+end
+
 function system_DAIB(nodes, branches, sources, invs)
     #Create system with BasePower = 100 MVA and nominal frequency 50 Hz.
     sys = PSY.System(100.0, frequency = 50.0)

test/data_tests/network_test_data.jl

@@ -45,6 +45,12 @@ nodes_DArco_IB() = [
     PSY.Bus(2, "Bus 2", "PV", 0, 1.0, (min = 0.94, max = 1.06), 0.69, nothing, nothing),
 ]
 
+nodes_multimachine() = [
+    PSY.Bus(1, "Bus 1", "REF", 0, 1.05, (min = 0.94, max = 1.06), 138, nothing, nothing),
+    PSY.Bus(2, "Bus 2", "PQ", 0, 1.00, (min = 0.94, max = 1.06), 138, nothing, nothing),
+    PSY.Bus(3, "Bus 3", "PV", 0, 1.00, (min = 0.94, max = 1.06), 138, nothing, nothing),
+]
+
 ############### Branches Data ########################
 
 branches_OMIB(nodes_OMIB) = [PSY.Line(
@@ -386,6 +392,60 @@ branch_3bus_case9_fault(nodes) = [
     ),
 ]
 
+branch_multimachine(nodes) = [
+    Line(
+        "Line1", #name
+        true, #available
+        0.0, #active power flow initial condition (from-to)
+        0.0, #reactive power flow initial condition (from-to)
+        Arc(from = nodes[1], to = nodes[2]), #Connection between buses
+        0.01, #resistance in pu
+        0.05, #reactance in pu
+        (from = 0.0, to = 0.0), #susceptance in pu
+        18.046, #rate in MW
+        1.04,
+    ),
+    Line(
+        "Line2", #name
+        true, #available
+        0.0, #active power flow initial condition (from-to)
+        0.0, #reactive power flow initial condition (from-to)
+        Arc(from = nodes[2], to = nodes[3]), #Connection between buses
+        0.01, #resistance in pu
+        0.05, #reactance in pu
+        (from = 0.0, to = 0.0), #susceptance in pu
+        18.046, #rate in MW
+        1.04,
+    ),  #angle limits (-min and max)
+]
+
+branch_multimachine_fault(nodes) = [
+    Line(
+        "Line1", #name
+        true, #available
+        0.0, #active power flow initial condition (from-to)
+        0.0, #reactive power flow initial condition (from-to)
+        Arc(from = nodes[1], to = nodes[2]), #Connection between buses
+        0.01, #resistance in pu
+        0.05, #reactance in pu
+        (from = 0.0, to = 0.0), #susceptance in pu
+        18.046, #rate in MW
+        1.04,
+    ),
+    Line(
+        "Line2", #name
+        true, #available
+        0.0, #active power flow initial condition (from-to)
+        0.0, #reactive power flow initial condition (from-to)
+        Arc(from = nodes[2], to = nodes[3]), #Connection between buses
+        0.02, #resistance in pu
+        0.10, #reactance in pu
+        (from = 0.0, to = 0.0), #susceptance in pu
+        18.046, #rate in MW
+        1.04,
+    ),  #angle limits (-min and max)
+]
+
 ############### Load Data ########################
 
 loads_OMIB(nodes_OMIB) = [PSY.PowerLoad(
@@ -527,6 +587,17 @@ loads_3bus_case9(nodes) = [
     PowerLoad("Bus3", true, nodes[3], PSY.LoadModels.ConstantPower, 0.3, 0.1, 0.5, 0.3),
 ]
 
+loads_multimachine(nodes) = [PowerLoad(
+    "LBus1", #name
+    true, #availability
+    nodes[2], #bus
+    LoadModels.ConstantPower, #type
+    0.8, #P
+    0.01, #Q
+    0.8, #P_max
+    0.01, #Q_max
+)]
+
 ############### Infinite Sources Data ########################
 
 inf_gen_1_pu(nodes) = PSY.Source(

test/test_case10_multimachine.jl

@@ -0,0 +1,82 @@
+"""
+Case 10:
+This case study a three bus system with 2 machine (Fixed EMF - Single Shaft: 2 State model), and a load in between.
+The perturbation trips one of the two circuits of line between buses 2 and 3, duplicating its impedance.
+"""
+
+##################################################
+############### LOAD DATA ########################
+##################################################
+
+############### Data Network ########################
+
+nodes_case10 = nodes_multimachine()
+
+branch_case10 = branch_multimachine(nodes_case10)
+
+#Trip of Line 1.
+branch_case10_fault = branch_multimachine_fault(nodes_case10)
+
+loads_case10 = loads_multimachine(nodes_case10)
+
+######## Machine Data #########
+
+### Case 9 Generators ###
+case10_gen1 = dyn_gen_case10_ref(nodes_case10)
+case10_gen2 = dyn_gen_case10_2(nodes_case10)
+
+######################### Dynamical System ########################
+
+sys = system_no_inv_no_sources(
+    nodes_case10,
+    branch_case10,
+    loads_case10,
+    [case10_gen1, case10_gen2],
+);
+
+##################################################
+############### SOLVE PROBLEM ####################
+##################################################
+
+#Compute Y_bus after fault
+Ybus_fault = Ybus(branch_case10_fault, nodes_case10)[:, :]
+
+#Time span
+tspan = (0.0, 30.0)
+
+Ybus_change = ThreePhaseFault(
+    1.0, #change at t = 1.0
+    Ybus_fault,
+); #New YBus
+
+x0_guess = [
+    1.02, #V1_R
+    1.005, #V2_R
+    1.0, #V3_R
+    0.0, #V1_I
+    0.01, #V2_I
+    0.01, #V3_I
+    0.1, #δ_1
+    1.0, #ω_1
+    0.1, #xg
+    1.0, #δ_2
+    0.0, #ω_2
+]
+
+sim = Simulation(
+    sys, #system
+    tspan, #time span
+    Ybus_change, #Type of Fault
+    initial_guess = x0_guess,
+)
+
+#Run simulation
+run_simulation!(
+    sim, #simulation structure
+    IDA(),#Sundials DAE Solver
+    dtmax = 0.02, #keywords arguments
+)
+
+series = get_state_series(sim, ("Case10Gen1", :ω));
+
+@test sim.solution.retcode == :Success