diff --git a/.gitignore b/.gitignore
index 46c911c..7fa4a21 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,7 @@
 #Files generated by invoking Julia with --code-coverage
 *.jl.cov
 *.jl.*.cov
+*.log
 
 # Files generated by invoking Julia with --track-allocation
 *.jl.mem
@@ -113,6 +114,5 @@ $RECYCLE.BIN/
 # Windows shortcuts
 *.lnk
 
-
 ## Acknowledgements
 # Many thanks to `https://gitignore.io/`, written and maintained by Joe Blau, which contributed much material to this gitignore file.
diff --git a/Project.toml b/Project.toml
index 1d120b9..9a995c5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -14,7 +15,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 [compat]
 DiffEqBase = "6"
 ForwardDiff = "~v0.10"
-InfrastructureSystems = "~0.6"
+InfrastructureSystems = "~0.6.3"
 NLsolve = "4"
-PowerSystems = "~0.16"
+PowerSystems = "~0.17"
 julia = "^1.2"
diff --git a/src/LITS.jl b/src/LITS.jl
index 9248502..f7d3167 100644
--- a/src/LITS.jl
+++ b/src/LITS.jl
@@ -19,6 +19,7 @@ export get_voltagemag_series
 export print_init_states
 
 ####################################### Package Imports ####################################
+import Logging
 import InfrastructureSystems
 import DiffEqBase
 import ForwardDiff
diff --git a/src/base/definitions.jl b/src/base/definitions.jl
index 76261e5..beea8a3 100644
--- a/src/base/definitions.jl
+++ b/src/base/definitions.jl
@@ -40,36 +40,48 @@ Inverter Inner Vars:
 md_var :: Modulation signal on the d-component
 mq_var :: Modulation signal on the q-component
 Vdc_var :: DC voltage supplied by the DC source
-Vd_filter_var :: Voltage seen in the capacitor of the filter in the d-component
-Vq_filter_var :: Voltage seen in the capacitor of the filter in the q-component
+Vr_filter_var :: Voltage seen in the capacitor of the filter in the R-component
+Vi_filter_var :: Voltage seen in the capacitor of the filter in the I-component
 ω_freq_estimator_var :: Frequency estimated by the frequency estimator.
 V_oc_var :: Control voltage supplied from the outer loop control to the inner loop
 ω_oc_var :: Control frequency supplied from the outer loop control the inner loop
 θ_oc_var :: Variation of the angle (PLL or VSM) of the inverter
 VR_inv_var :: Real terminal voltage on the inverter
 VI_inv_var :: Imaginary terminal voltage on the inverter
-Vd_cnv_var :: Voltage supplied from the converter in the d-component
-Vq_cnv_var :: Voltage supplied from the converter in the q-component
+Vr_cnv_var :: Voltage supplied from the converter in the R-component
+Vi_cnv_var :: Voltage supplied from the converter in the I-component
 """
 
 @enum inverter_inner_vars begin
     md_var = 1
     mq_var = 2
     Vdc_var = 3
-    Vd_filter_var = 4
-    Vq_filter_var = 5
+    Vr_filter_var = 4
+    Vi_filter_var = 5
     ω_freq_estimator_var = 6
     V_oc_var = 7
     ω_oc_var = 8
     θ_oc_var = 9
     VR_inv_var = 10
     VI_inv_var = 11
-    Vd_cnv_var = 12
-    Vq_cnv_var = 13
+    Vr_cnv_var = 12
+    Vi_cnv_var = 13
 end
 
 Base.to_index(ix::inverter_inner_vars) = Int(ix)
 
+@enum dq_ref begin
+    d = 1
+    q = 2
+end
+@enum RI_ref begin
+    R = 1
+    I = 2
+end
+
+Base.to_index(ix::dq_ref) = Int(ix)
+Base.to_index(ix::RI_ref) = Int(ix)
+
 const V_ref_index = 1
 const ω_ref_index = 2
 const P_ref_index = 3
diff --git a/src/base/ports.jl b/src/base/ports.jl
index a10f180..e317ee4 100644
--- a/src/base/ports.jl
+++ b/src/base/ports.jl
@@ -45,7 +45,7 @@ end
 #### Converter Ports ####
 function Ports(::PSY.Converter)
     state_input = Vector{Symbol}()
-    inner_input = [md_var, mq_var, Vdc_var, Vd_cnv_var, Vq_cnv_var]
+    inner_input = [md_var, mq_var, Vdc_var, Vr_cnv_var, Vi_cnv_var]
     return Ports(state_input, inner_input)
 end
 
@@ -63,11 +63,11 @@ function Ports(::PSY.Filter)
     inner_input = [
         VR_inv_var,
         VI_inv_var,
-        Vd_cnv_var,
-        Vq_cnv_var,
+        Vr_cnv_var,
+        Vi_cnv_var,
         θ_oc_var,
-        Vd_filter_var,
-        Vq_filter_var,
+        Vr_filter_var,
+        Vi_filter_var,
     ]
     return Ports(state_input, inner_input)
 end
@@ -75,22 +75,22 @@ end
 #### Freq. Estimator Ports ####
 
 function Ports(::PSY.FrequencyEstimator)
-    state_input = [:vd_filter, :vq_filter, :θ_oc]
+    state_input = [:vr_filter, :vi_filter, :θ_oc]
     #TODO: Move PLL to PCC, i.e. move v_cap (D'Arco v_o), to inner inputs
-    inner_input = [Vd_filter_var, Vq_filter_var, θ_oc_var, ω_freq_estimator_var]
+    inner_input = [Vr_filter_var, Vi_filter_var, θ_oc_var, ω_freq_estimator_var]
     return Ports(state_input, inner_input)
 end
 
 #### Outer Control Ports ####
 function Ports(::PSY.OuterControl)
-    state_input = [:vd_pll, :vq_pll, :ε_pll, :vd_filter, :vq_filter, :id_filter, :iq_filter]
-    inner_input = [Vd_filter_var, Vq_filter_var, ω_freq_estimator_var]
+    state_input = [:vd_pll, :vq_pll, :ε_pll, :vr_filter, :vi_filter, :ir_filter, :ii_filter]
+    inner_input = [Vr_filter_var, Vi_filter_var, ω_freq_estimator_var]
     return Ports(state_input, inner_input)
 end
 
 #### Inner Control Ports ####
 function Ports(::PSY.InnerControl)
-    state_input = [:id_filter, :iq_filter, :id_cnv, :iq_cnv, :vd_filter, :vq_filter]
-    inner_input = [Vd_filter_var, Vq_filter_var, V_oc_var, ω_oc_var]
+    state_input = [:ir_filter, :ii_filter, :ir_cnv, :ii_cnv, :vr_filter, :vi_filter]
+    inner_input = [Vr_filter_var, Vi_filter_var, V_oc_var, ω_oc_var]
     return Ports(state_input, inner_input)
 end
diff --git a/src/models/inverter_models/converter_models.jl b/src/models/inverter_models/converter_models.jl
index 4aebba6..87d0923 100644
--- a/src/models/inverter_models/converter_models.jl
+++ b/src/models/inverter_models/converter_models.jl
@@ -13,13 +13,13 @@ function mdl_converter_ode!(
     #Obtain inner variables for component
     md = get_inner_vars(device)[md_var]
     mq = get_inner_vars(device)[mq_var]
-    VDC = get_inner_vars(device)[Vdc_var]
+    Vdc = get_inner_vars(device)[Vdc_var]
+    θ_oc = get_inner_vars(device)[θ_oc_var]
+
+    #Transform reference frame to grid reference frame
+    m_ri = dq_ri(θ_oc + pi / 2) * [md; mq]
 
     #Update inner_vars
-    get_inner_vars(device)[Vd_cnv_var] = md * VDC
-    get_inner_vars(device)[Vq_cnv_var] = mq * VDC
+    get_inner_vars(device)[Vr_cnv_var] = m_ri[R] * Vdc
+    get_inner_vars(device)[Vi_cnv_var] = m_ri[I] * Vdc
 end
-
-#TODO: Same as above, but:
-## VDC is state of DC source
-## IDC is algebraic inner_var output from Vcnv real output power
diff --git a/src/models/inverter_models/filter_models.jl b/src/models/inverter_models/filter_models.jl
index 9ee8234..5b11944 100644
--- a/src/models/inverter_models/filter_models.jl
+++ b/src/models/inverter_models/filter_models.jl
@@ -15,19 +15,14 @@ function mdl_filter_ode!(
     P <: PSY.FrequencyEstimator,
 }
 
-    #Obtain external states inputs for component
-    #TODO: If converter has dynamics, need to reference states:
-    #external_ix = device.input_port_mapping[device.converter]
-    #Vd_cnv = device_states[external_ix[1]]
-    #Vq_cnv = device_states[external_ix[2]]
-    external_ix = get_input_port_ix(device, PSY.LCLFilter)
-    δ = device_states[external_ix[1]]
+    #external_ix = get_input_port_ix(device, PSY.LCLFilter)
+    #θ_oc = device_states[external_ix[1]]
 
     #Obtain inner variables for component
     V_tR = get_inner_vars(device)[VR_inv_var]
     V_tI = get_inner_vars(device)[VI_inv_var]
-    Vd_cnv = get_inner_vars(device)[Vd_cnv_var]
-    Vq_cnv = get_inner_vars(device)[Vq_cnv_var]
+    Vr_cnv = get_inner_vars(device)[Vr_cnv_var]
+    Vi_cnv = get_inner_vars(device)[Vi_cnv_var]
 
     #Get parameters
     filter = PSY.get_filter(device)
@@ -39,21 +34,17 @@ function mdl_filter_ode!(
     rg = PSY.get_rg(filter)
     MVABase = PSY.get_inverter_Sbase(device)
 
-    #RI to dq transformation
-    V_dq = ri_dq(δ + pi / 2) * [V_tR; V_tI]
-    V_g = sqrt(V_tR^2 + V_tI^2)
-
     #Obtain indices for component w/r to device
     local_ix = get_local_state_ix(device, PSY.LCLFilter)
 
     #Define internal states for filter
     internal_states = @view device_states[local_ix]
-    Id_cnv = internal_states[1]
-    Iq_cnv = internal_states[2]
-    Vd_filter = internal_states[3]
-    Vq_filter = internal_states[4]
-    Id_filter = internal_states[5]
-    Iq_filter = internal_states[6]
+    Ir_cnv = internal_states[1]
+    Ii_cnv = internal_states[2]
+    Vr_filter = internal_states[3]
+    Vi_filter = internal_states[4]
+    Ir_filter = internal_states[5]
+    Ii_filter = internal_states[6]
 
     #Inputs (control signals) - N/A
 
@@ -61,42 +52,40 @@ function mdl_filter_ode!(
     #Inverter Output Inductor (internal state)
     #𝜕id_c/𝜕t
     output_ode[local_ix[1]] = (
-        ωb / lf * Vd_cnv - ωb / lf * Vd_filter - ωb * rf / lf * Id_cnv +
-        ωb * ω_sys * Iq_cnv
+        ωb / lf * Vr_cnv - ωb / lf * Vr_filter - ωb * rf / lf * Ir_cnv +
+        ωb * ω_sys * Ii_cnv
     )
     #𝜕iq_c/𝜕t
     output_ode[local_ix[2]] = (
-        ωb / lf * Vq_cnv - ωb / lf * Vq_filter - ωb * rf / lf * Iq_cnv -
-        ωb * ω_sys * Id_cnv
+        ωb / lf * Vi_cnv - ωb / lf * Vi_filter - ωb * rf / lf * Ii_cnv -
+        ωb * ω_sys * Ir_cnv
     )
     #LCL Capacitor (internal state)
     #𝜕vd_o/𝜕t
     output_ode[local_ix[3]] =
-        (ωb / cf * Id_cnv - ωb / cf * Id_filter + ωb * ω_sys * Vq_filter)
+        (ωb / cf * Ir_cnv - ωb / cf * Ir_filter + ωb * ω_sys * Vi_filter)
     #𝜕vq_o/𝜕t
     output_ode[local_ix[4]] =
-        (ωb / cf * Iq_cnv - ωb / cf * Iq_filter - ωb * ω_sys * Vd_filter)
+        (ωb / cf * Ii_cnv - ωb / cf * Ii_filter - ωb * ω_sys * Vr_filter)
     #Grid Inductance (internal state)
     #𝜕id_o/𝜕t
     output_ode[local_ix[5]] = (
-        ωb / lg * Vd_filter - ωb / lg * V_dq[1] - ωb * rg / lg * Id_filter +
-        ωb * ω_sys * Iq_filter
+        ωb / lg * Vr_filter - ωb / lg * V_tR - ωb * rg / lg * Ir_filter +
+        ωb * ω_sys * Ii_filter
     )
-    #𝜕iq_o/𝜕t (Multiply Vq by -1 to lag instead of lead)
+    #𝜕iq_o/𝜕t
     output_ode[local_ix[6]] = (
-        ωb / lg * Vq_filter + ωb / lg * (-V_dq[2]) - ωb * rg / lg * Iq_filter -
-        ωb * ω_sys * Id_filter
+        ωb / lg * Vi_filter - ωb / lg * V_tI - ωb * rg / lg * Ii_filter -
+        ωb * ω_sys * Ir_filter
     )
 
     #Update inner_vars
-    get_inner_vars(device)[Vd_filter_var] = Vd_filter
-    get_inner_vars(device)[Vq_filter_var] = Vq_filter
-    #TODO: If PLL models at PCC, need to update inner vars:
-    #get_inner_vars(device)[Vd_filter_var] = V_dq[q::dq_ref]
-    #get_inner_vars(device)[Vq_filter_var] = V_dq[q::dq_ref]
+    get_inner_vars(device)[Vr_filter_var] = Vr_filter
+    get_inner_vars(device)[Vi_filter_var] = Vi_filter
 
     #Compute current from the inverter to the grid
-    I_RI = (MVABase / sys_Sbase) * dq_ri(δ + pi / 2) * [Id_filter; Iq_filter]
+    I_RI = (MVABase / sys_Sbase) * [Ir_filter; Ii_filter]
+
     #Update current
     current_r[1] += I_RI[1]
     current_i[1] += I_RI[2]
diff --git a/src/models/inverter_models/frequency_estimator_models.jl b/src/models/inverter_models/frequency_estimator_models.jl
index e943424..c3bd5a2 100644
--- a/src/models/inverter_models/frequency_estimator_models.jl
+++ b/src/models/inverter_models/frequency_estimator_models.jl
@@ -14,14 +14,11 @@ function mdl_freq_estimator_ode!(
 
     #Obtain external states inputs for component
     external_ix = get_input_port_ix(device, PSY.KauraPLL)
-    Vd_filter = device_states[external_ix[1]]
-    Vq_filter = device_states[external_ix[2]]
-    θ_oc = device_states[external_ix[3]]
+    Vr_filter = device_states[external_ix[1]]
+    Vi_filter = device_states[external_ix[2]]
 
-    #Obtain inner variables for component
-    #Vd_filter = device.inner_vars[Vd_filter_var]
-    #Vq_filter = device.inner_vars[Vq_filter_var]
-    #θ_oc = device.inner_vars[θ_oc_var]
+    #V_tR = get_inner_vars(device)[VR_inv_var]
+    #V_tI = get_inner_vars(device)[VI_inv_var]
 
     #Get parameters
     pll_control = PSY.get_freq_estimator(device)
@@ -40,20 +37,17 @@ function mdl_freq_estimator_ode!(
     ϵ_pll = internal_states[3]
     θ_pll = internal_states[4]
 
+    #Transform to internal dq-PLL reference frame
+    V_dq_pll = ri_dq(θ_pll + pi / 2) * [Vr_filter; Vi_filter]
+
     #Inputs (control signals)
 
     #Compute 6 states ODEs (D'Arco EPSR122 Model)
     #Output Voltage LPF (internal state)
     #𝜕vpll_d/𝜕t, D'Arco ESPR122 eqn. 12
-    output_ode[local_ix[1]] = (
-        ω_lp * Vd_filter * cos(θ_pll - θ_oc) + ω_lp * Vq_filter * sin(θ_pll - θ_oc) -
-        ω_lp * vpll_d
-    )
+    output_ode[local_ix[1]] = ω_lp * (V_dq_pll[d] - vpll_d)
     #𝜕vpll_q/𝜕t, D'Arco ESPR122 eqn. 12
-    output_ode[local_ix[2]] = (
-        -ω_lp * Vd_filter * sin(θ_pll - θ_oc) + ω_lp * Vq_filter * cos(θ_pll - θ_oc) -
-        ω_lp * vpll_q
-    )
+    output_ode[local_ix[2]] = ω_lp * (V_dq_pll[q] - vpll_q)
     #PI Integrator (internal state)
     #𝜕dϵ_pll/𝜕t, D'Arco ESPR122 eqn. 13
     output_ode[local_ix[3]] = atan(vpll_q / vpll_d)
diff --git a/src/models/inverter_models/inner_control_models.jl b/src/models/inverter_models/inner_control_models.jl
index e7a1384..db5b78c 100644
--- a/src/models/inverter_models/inner_control_models.jl
+++ b/src/models/inverter_models/inner_control_models.jl
@@ -12,19 +12,20 @@ function mdl_inner_ode!(
 
     #Obtain external states inputs for component
     external_ix = get_input_port_ix(device, PSY.CurrentControl)
-    Id_filter = device_states[external_ix[1]]
-    Iq_filter = device_states[external_ix[2]]
-    Id_cnv = device_states[external_ix[3]]
-    Iq_cnv = device_states[external_ix[4]]
-    Vd_filter = device_states[external_ix[5]] #TODO: Should be inner reference after initialization
-    Vq_filter = device_states[external_ix[6]] #TODO: Should be inner reference after initialization
+    Ir_filter = device_states[external_ix[1]]
+    Ii_filter = device_states[external_ix[2]]
+    Ir_cnv = device_states[external_ix[3]]
+    Ii_cnv = device_states[external_ix[4]]
+    Vr_filter = device_states[external_ix[5]] #TODO: Should be inner reference after initialization
+    Vi_filter = device_states[external_ix[6]] #TODO: Should be inner reference after initialization
 
     #Obtain inner variables for component
-    #Vd_filter = get_inner_vars(device)[Vd_filter_var]
-    #Vq_filter = get_inner_vars(device)[Vq_filter_var]
+    #Vd_filter = get_inner_vars(device)[Vr_filter_var]
+    #Vq_filter = get_inner_vars(device)[Vi_filter_var]
     ω_oc = get_inner_vars(device)[ω_oc_var]
+    θ_oc = get_inner_vars(device)[θ_oc_var]
     v_refr = get_inner_vars(device)[V_oc_var]
-    vdc = get_inner_vars(device)[Vdc_var]
+    Vdc = get_inner_vars(device)[Vdc_var]
 
     #Get Voltage Controller parameters
     inner_control = PSY.get_inner_control(device)
@@ -56,55 +57,59 @@ function mdl_inner_ode!(
     ϕ_d = internal_states[5]
     ϕ_q = internal_states[6]
 
-    #Inputs (control signals)
+    #Transformations to dq frame
+    I_dq_filter = ri_dq(θ_oc + pi / 2) * [Ir_filter; Ii_filter]
+    I_dq_cnv = ri_dq(θ_oc + pi / 2) * [Ir_cnv; Ii_cnv]
+    V_dq_filter = ri_dq(θ_oc + pi / 2) * [Vr_filter; Vi_filter]
 
     ### Compute 6 states ODEs (D'Arco EPSR122 Model) ###
     ## SRF Voltage Control w/ Virtual Impedance ##
-    #Virtual Impedance - Computation but not DAE
-    #v_refr = V_ref + kq*(q_ref - qm)
-    Vd_filter_ref = (v_refr - rv * Id_filter + ω_oc * lv * Iq_filter)
-    Vq_filter_ref = (-rv * Iq_filter - ω_oc * lv * Id_filter)
+    #Virtual Impedance
+    Vd_filter_ref = (v_refr - rv * I_dq_filter[d] + ω_oc * lv * I_dq_filter[q])
+    Vq_filter_ref = (-rv * I_dq_filter[q] - ω_oc * lv * I_dq_filter[d])
+
+    #Voltage Control ODEs
+    #PI Integrator (internal state)
+    output_ode[local_ix[1]] = (Vd_filter_ref - V_dq_filter[d])
+    output_ode[local_ix[2]] = (Vq_filter_ref - V_dq_filter[q])
+
     #Output Control Signal - Links to SRF Current Controller
     Id_cnv_ref = (
-        kpv * (Vd_filter_ref - Vd_filter) + kiv * ξ_d - cf * ω_oc * Vq_filter +
-        kffi * Id_filter
+        kpv * (Vd_filter_ref - V_dq_filter[d]) + kiv * ξ_d - cf * ω_oc * V_dq_filter[q] + kffi * I_dq_filter[d]
     )
 
     Iq_cnv_ref = (
-        kpv * (Vq_filter_ref - Vq_filter) +
+        kpv * (Vq_filter_ref - V_dq_filter[q]) +
         kiv * ξ_q +
-        cf * ω_oc * Vd_filter +
-        kffi * Iq_filter
+        cf * ω_oc * V_dq_filter[d] +
+        kffi * I_dq_filter[q]
     )
-    #Voltage Control ODEs
-    #PI Integrator (internal state)
-    output_ode[local_ix[1]] = (Vd_filter_ref - Vd_filter)
-    output_ode[local_ix[2]] = (Vq_filter_ref - Vq_filter)
 
     ## SRF Current Control ##
-    #Active Damping
-    #vad_d = kad*(Vd_filter-ϕ_d)
-    #vad_q = kad*(Vq_filter-ϕ_q)
+    #Current Control ODEs
+    #PI Integrator (internal state)
+    output_ode[local_ix[3]] = Id_cnv_ref - I_dq_cnv[d]
+    output_ode[local_ix[4]] = Iq_cnv_ref - I_dq_cnv[q]
+
     #References for Converter Output Voltage
     Vd_cnv_ref = (
-        kpc * (Id_cnv_ref - Id_cnv) + kic * γ_d - ω_oc * lf * Iq_cnv + kffv * Vd_filter - kad * (Vd_filter - ϕ_d)
+        kpc * (Id_cnv_ref - I_dq_cnv[d]) + kic * γ_d - ω_oc * lf * I_dq_cnv[q] +
+        kffv * V_dq_filter[d] - kad * (V_dq_filter[d] - ϕ_d)
     )
     Vq_cnv_ref = (
-        kpc * (Iq_cnv_ref - Iq_cnv) + kic * γ_q + ω_oc * lf * Id_cnv + kffv * Vq_filter - kad * (Vq_filter - ϕ_q)
+        kpc * (Iq_cnv_ref - I_dq_cnv[q]) +
+        kic * γ_q +
+        ω_oc * lf * I_dq_cnv[d] +
+        kffv * V_dq_filter[q] - kad * (V_dq_filter[q] - ϕ_q)
     )
-    #Modulation Commands to Converter
-    #md = Vd_cnv_ref/vdc
-    #mq = Vq_cnv_ref/vdc
-    #Current Control ODEs
-    #PI Integrator (internal state)
-    output_ode[local_ix[3]] = Id_cnv_ref - Id_cnv
-    output_ode[local_ix[4]] = Iq_cnv_ref - Iq_cnv
+
     #Active Damping LPF (internal state)
-    output_ode[local_ix[5]] = ωad * Vd_filter - ωad * ϕ_d
-    output_ode[local_ix[6]] = ωad * Vq_filter - ωad * ϕ_q
+    output_ode[local_ix[5]] = ωad * V_dq_filter[d] - ωad * ϕ_d
+    output_ode[local_ix[6]] = ωad * V_dq_filter[q] - ωad * ϕ_q
 
     #Update inner_vars
-    get_inner_vars(device)[md_var] = Vd_cnv_ref / vdc
-    get_inner_vars(device)[mq_var] = Vq_cnv_ref / vdc
+    #Modulation Commands to Converter
+    get_inner_vars(device)[md_var] = Vd_cnv_ref / Vdc
+    get_inner_vars(device)[mq_var] = Vq_cnv_ref / Vdc
 
 end
diff --git a/src/models/inverter_models/outer_control_models.jl b/src/models/inverter_models/outer_control_models.jl
index ede1b74..79cf190 100644
--- a/src/models/inverter_models/outer_control_models.jl
+++ b/src/models/inverter_models/outer_control_models.jl
@@ -27,13 +27,15 @@ function mdl_outer_ode!(
     vpll_d = device_states[external_ix[1]]
     vpll_q = device_states[external_ix[2]]
     ϵ_pll = device_states[external_ix[3]]
-    Vd_filter = device_states[external_ix[4]] #TODO: Should be inner reference after initialization
-    Vq_filter = device_states[external_ix[5]] #TODO: Should be inner reference after initialization
-    Id_filter = device_states[external_ix[6]]
-    Iq_filter = device_states[external_ix[7]]
+    Vr_filter = device_states[external_ix[4]]
+    Vi_filter = device_states[external_ix[5]]
+    Ir_filter = device_states[external_ix[6]]
+    Ii_filter = device_states[external_ix[7]]
 
     #Obtain inner variables for component
     ω_pll = get_inner_vars(device)[ω_freq_estimator_var]
+    V_tR = get_inner_vars(device)[VR_inv_var]
+    V_tI = get_inner_vars(device)[VI_inv_var]
 
     #Get Active Power Controller parameters
     outer_control = PSY.get_outer_control(device)
@@ -70,24 +72,17 @@ function mdl_outer_ode!(
     qm = internal_states[3]
 
     #Obtain additional expressions
-    p_elec_out = Id_filter * Vd_filter + Iq_filter * Vq_filter
+    p_elec_out = Ir_filter * Vr_filter + Ii_filter * Vi_filter
+    q_elec_out = -Ii_filter * Vr_filter + Ir_filter * Vi_filter
 
     #Compute 3 states ODEs
     output_ode[local_ix[1]] =
         (p_ref / Ta - p_elec_out / Ta - kd * (ω_oc - ω_pll) / Ta - kω * (ω_oc - ω_ref) / Ta)
     output_ode[local_ix[2]] = ωb * (ω_oc - ω_sys)
-    output_ode[local_ix[3]] =
-        (-ωf * Iq_filter * Vd_filter + ωf * Id_filter * Vq_filter - ωf * qm)
+    output_ode[local_ix[3]] = (ωf * (q_elec_out - qm))
 
     #Update inner vars
     get_inner_vars(device)[θ_oc_var] = θ_oc
     get_inner_vars(device)[ω_oc_var] = ω_oc
     get_inner_vars(device)[V_oc_var] = V_ref + kq * (q_ref - qm)
 end
-#output_ode[local_ix[1]] = (
-#    -Id_filter * Vd_filter / Ta - Iq_filter * Vq_filter / Ta +
-#    kd * kp_pll * atan(vpll_q / vpll_d) / Ta +
-#    kd * ki_pll * ϵ_pll / Ta - (kd + kω) * (ω_oc - ω_sys) / Ta +
-#    p_ref / Ta +
-#    kω * ω_ref / Ta - kω * ω_sys / Ta
-#)
diff --git a/src/models/load_models.jl b/src/models/load_models.jl
index 9dc3296..c532cda 100644
--- a/src/models/load_models.jl
+++ b/src/models/load_models.jl
@@ -7,25 +7,25 @@ function mdl_Zload!(
     sys::PSY.System,
 )
 
-    #Get parameters
+    #Load squared voltage magnitude at steady state
+    bus = PSY.get_bus(device)
+    #Vmag_sq = PSY.get_voltage(bus)^2
+    Vmag_sq = 1.0
+
+    #Load device parameters
     P = PSY.get_activepower(device)
     Q = PSY.get_reactivepower(device)
 
-    #Compute impedance and load
-    #Z_load =  1.0/(P-Q*1im) #in pu
+    #Compute impedance and load: Z = |V|^2/conj(S)
+    #Z_load =  Vmag_sq/(P-Q*1im) #in pu
     #I = -(voltage_r[1] + voltage_i[1]*1im)/Z_load #in pu flowing out
     #current_r[1] += real(I)
     #current_i[1] += imag(I)
 
-    #Compute current given constant power
-    #I = (P-Q*1im)/(voltage_r[1] - voltage_i[1]*1im)
-    #current_r[1] += real(I)
-    #current_i[1] += imag(I)
-
     #Update current
-    #This model creates an equivalent RL/RC circuit assuming nominal voltage (1 pu)
-    current_r[1] += -(voltage_r[1] * P + voltage_i[1] * Q) #in system pu flowing out
-    current_i[1] += -(voltage_i[1] * P - voltage_r[1] * Q) #in system pu flowing out
+    #This model creates an equivalent RL/RC circuit based on steady state voltage
+    current_r[1] += -(1.0 / Vmag_sq) * (voltage_r[1] * P + voltage_i[1] * Q) #in system pu flowing out
+    current_i[1] += -(1.0 / Vmag_sq) * (voltage_i[1] * P - voltage_r[1] * Q) #in system pu flowing out
 
     return
 end
diff --git a/src/models/ref_transformations.jl b/src/models/ref_transformations.jl
index 537fbc8..a68706e 100644
--- a/src/models/ref_transformations.jl
+++ b/src/models/ref_transformations.jl
@@ -1,28 +1,4 @@
-@enum dq_ref begin
-    q = 1
-    d = 2
-end
-@enum RI_ref begin
-    R = 1
-    I = 2
-end
-function dq_ri(δ::Real)
-    ## Uses the referenceframe of the Kundur page 852 of dq to RI
-    return [
-        sin(δ) cos(δ)
-        -cos(δ) sin(δ)
-    ]
-end
-
-function ri_dq(δ::Real)
-    #Uses the reference frame of the Kundur page 852 of RI to dq
-    return [
-        sin(δ) -cos(δ)
-        cos(δ) sin(δ)
-    ]
-end
-
-function dq_ri(δ::Float64)
+function dq_ri(δ)
     ## Uses the referenceframe of the Kundur page 852 of dq to RI
     return [
         sin(δ) cos(δ)
@@ -30,7 +6,7 @@ function dq_ri(δ::Float64)
     ]
 end
 
-function ri_dq(δ::Float64)
+function ri_dq(δ)
     #Uses the reference frame of the Kundur page 852 of RI to dq
     return [
         sin(δ) -cos(δ)
diff --git a/src/models/source_models.jl b/src/models/source_models.jl
index 7deecb6..80c46f1 100644
--- a/src/models/source_models.jl
+++ b/src/models/source_models.jl
@@ -9,8 +9,8 @@ function mdl_source!(
 
     #Load device parameters
     bus = PSY.get_bus(device)
-    V_R = PSY.get_voltage(bus)*cos(PSY.get_angle(bus))
-    V_I = PSY.get_voltage(bus)*sin(PSY.get_angle(bus))
+    V_R = PSY.get_voltage(bus) * cos(PSY.get_angle(bus))
+    V_I = PSY.get_voltage(bus) * sin(PSY.get_angle(bus))
     X_th = PSY.get_X_th(device)
 
     #I = ( (V_R + V_I*1im) - (V_tR + V_tI*1im) )/(X_th*1im)
diff --git a/src/models/system.jl b/src/models/system.jl
index 88e379f..a77705a 100644
--- a/src/models/system.jl
+++ b/src/models/system.jl
@@ -51,8 +51,20 @@ function system!(out::Vector{<:Real}, dx, x, sys::PSY.System, t::Float64)
         out[ix_range] = injection_ode[ode_range] - dx[ix_range]
     end
 
-    # TODO: remove filtered get_component
-    for d in PSY.get_components(PSY.StaticInjection, sys, d -> !isa(d, PSY.Generator))
+    for d in PSY.get_components(PSY.ElectricLoad, 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
+
+    for d in PSY.get_components(PSY.Source, sys)
         bus_n = PSY.get_number(PSY.get_bus(d))
 
         device!(
diff --git a/src/utils/logging.jl b/src/utils/logging.jl
new file mode 100644
index 0000000..1661551
--- /dev/null
+++ b/src/utils/logging.jl
@@ -0,0 +1,41 @@
+"""
+    configure_logging(;
+        console_level = Logging.Error,
+        file_level = Logging.Info,
+        filename = "power-simulations.log",
+    )
+
+Creates console and file loggers.
+
+**Note:** Log messages may not be written to the file until flush() or close() is called on
+the returned logger.
+
+# Arguments
+- `console_level = Logging.Error`: level for console messages
+- `file_level = Logging.Info`: level for file messages
+- `filename::String = power-simulations.log`: log file
+
+# Example
+```julia
+logger = configure_logging(console_level = Logging.Info)
+@info "log message"
+close(logger)
+```
+"""
+function configure_logging(;
+    console_level = Logging.Error,
+    file_level = Logging.Info,
+    filename = "power-simulations.log",
+)
+    return IS.configure_logging(
+        console = true,
+        console_stream = stderr,
+        console_level = console_level,
+        file = true,
+        filename = filename,
+        file_level = file_level,
+        file_mode = "w+",
+        tracker = nothing,
+        set_global = true,
+    )
+end
diff --git a/test/Project.toml b/test/Project.toml
index bdf85ea..595921e 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -8,3 +8,4 @@ PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
diff --git a/test/data_tests/OMIB.raw b/test/data_tests/OMIB.raw
new file mode 100644
index 0000000..9a43cf8
--- /dev/null
+++ b/test/data_tests/OMIB.raw
@@ -0,0 +1,28 @@
+0, 100, 33, 0, 0, 60  / 24-Apr-2020 17:05:49 - MATPOWER 7.0.1-dev
+
+
+     1, 'BUS 1       ',       230, 3,    1,    1, 1,        1.02,           0, 1.06, 0.94, 1.06, 0.94
+     2, 'BUS 2       ',       230, 2,    1,    1, 1,        1.05,           0, 1.06, 0.94, 1.06, 0.94
+0 / END OF BUS DATA, BEGIN LOAD DATA
+     2,  1, 1,    1,    1,        30,         1, 0, 0, 0, 0, 1, 1, 0
+0 / END OF LOAD DATA, BEGIN FIXED SHUNT DATA
+0 / END OF FIXED SHUNT DATA, BEGIN GENERATOR DATA
+     2,  1,        50,         0,       100,      -100,     1.02, 0,     100, 0, 1, 0, 0, 1, 1, 100,       100,         0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1
+0 / END OF GENERATOR DATA, BEGIN BRANCH DATA
+     1,      2, 1,    0.01,     0.05,   0.001,     100,     100,     100, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1
+0 / END OF BRANCH DATA, BEGIN TRANSFORMER DATA
+0 / END OF TRANSFORMER DATA, BEGIN AREA DATA
+0 / END OF AREA DATA, BEGIN TWO-TERMINAL DC DATA
+0 / END OF TWO-TERMINAL DC DATA, BEGIN VOLTAGE SOURCE CONVERTER DATA
+0 / END OF VOLTAGE SOURCE CONVERTER DATA, BEGIN IMPEDANCE CORRECTION DATA
+0 / END OF IMPEDANCE CORRECTION DATA, BEGIN MULTI-TERMINAL DC DATA
+0 / END OF MULTI-TERMINAL DC 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 CONTROL DEVICE DATA
+0 / END OF FACTS CONTROL 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
diff --git a/test/data_tests/ThreeBusNetwork.raw b/test/data_tests/ThreeBusNetwork.raw
new file mode 100644
index 0000000..d88af1f
--- /dev/null
+++ b/test/data_tests/ThreeBusNetwork.raw
@@ -0,0 +1,34 @@
+0, 100, 33, 0, 0, 60  / 24-Apr-2020 19:28:39 - MATPOWER 7.0.1-dev
+
+
+     1, 'BUS 1       ',       138, 3,    1,    1, 1,           1.02,        0,  1.1,  0.9,  1.1,  0.9
+     2, 'BUS 2       ',       138, 2,    1,    1, 1,           1,           0,  1.1,  0.9,  1.1,  0.9
+     3, 'BUS 3       ',       138, 2,    1,    1, 1,           1,           0,  1.1,  0.9,  1.1,  0.9
+0 / END OF BUS DATA, BEGIN LOAD DATA
+     1,  1, 1,    1,    1,        50,       15, 0, 0, 0, 0, 1, 1, 0
+     2,  1, 1,    1,    1,       170,       15, 0, 0, 0, 0, 1, 1, 0
+     3,  1, 1,    1,    1,       200,       5, 0, 0, 0, 0, 1, 1, 0
+0 / END OF LOAD DATA, BEGIN FIXED SHUNT DATA
+0 / END OF FIXED SHUNT DATA, BEGIN GENERATOR DATA
+     2,  1,       200,         0,       100,      -100,   1.0142, 0,     100, 0, 1, 0, 0, 1, 1, 100,       318,         0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1
+     3,  1,       200,         0,       100,      -100,   1.0059, 0,     100, 0, 1, 0, 0, 1, 1, 100,       318,         0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1
+0 / END OF GENERATOR DATA, BEGIN BRANCH DATA
+     1,      2, 1,  0.01008,     0.12,      0.2,     250,     250,     250, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1
+     1,      3, 1,  0.01008,     0.12,      0.2,     250,     250,     250, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1
+     2,      3, 1,  0.01008,     0.12,        1,     250,     250,     250, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1
+0 / END OF BRANCH DATA, BEGIN TRANSFORMER DATA
+0 / END OF TRANSFORMER DATA, BEGIN AREA DATA
+0 / END OF AREA DATA, BEGIN TWO-TERMINAL DC DATA
+0 / END OF TWO-TERMINAL DC DATA, BEGIN VOLTAGE SOURCE CONVERTER DATA
+0 / END OF VOLTAGE SOURCE CONVERTER DATA, BEGIN IMPEDANCE CORRECTION DATA
+0 / END OF IMPEDANCE CORRECTION DATA, BEGIN MULTI-TERMINAL DC DATA
+0 / END OF MULTI-TERMINAL DC 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 CONTROL DEVICE DATA
+0 / END OF FACTS CONTROL 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
diff --git a/test/data_tests/data_utils.jl b/test/data_tests/data_utils.jl
new file mode 100644
index 0000000..c6d9d61
--- /dev/null
+++ b/test/data_tests/data_utils.jl
@@ -0,0 +1,19 @@
+function add_source_to_ref(sys::PSY.System)
+    for g in PSY.get_components(StaticInjection, sys)
+        isa(g, ElectricLoad) && continue
+        g.bus.bustype == BusTypes.REF &&
+            error("A device is already attached to the REF bus")
+    end
+
+    slack_bus = [b for b in PSY.get_components(Bus, sys) if b.bustype == BusTypes.REF][1]
+    inf_source = Source(
+        name = "InfBus", #name
+        available = true, #availability
+        activepower = 0.0,
+        reactivepower = 0.0,
+        bus = slack_bus, #bus
+        X_th = 0.000005,
+    ) #Xth
+    PSY.add_component!(sys, inf_source)
+    return
+end
diff --git a/test/data_tests/dynamic_test_data.jl b/test/data_tests/dynamic_test_data.jl
index 87bb73d..5cb1b17 100644
--- a/test/data_tests/dynamic_test_data.jl
+++ b/test/data_tests/dynamic_test_data.jl
@@ -572,6 +572,37 @@ end
 ######## System Constructors #########
 ######################################
 
+function system_OMIB(nodes, branches, loads, sources, 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 infinite source
+    for source in sources
+        PSY.add_component!(sys, source)
+    end
+
+    #Add generator
+    for gen in gens
+        PSY.add_component!(sys, gen)
+    end
+    return sys
+end
+
 function system_no_inv(nodes, branches, loads, sources, gens)
     #Create system with BasePower = 100 MVA and nominal frequency 60 Hz.
     sys = PSY.System(100.0, frequency = 60.0)
diff --git a/test/data_tests/network_test_data.jl b/test/data_tests/network_test_data.jl
index b91af9b..e66634f 100644
--- a/test/data_tests/network_test_data.jl
+++ b/test/data_tests/network_test_data.jl
@@ -61,7 +61,7 @@ branches_OMIB(nodes_OMIB) = [PSY.Line(
     Arc(from = nodes_OMIB[1], to = nodes_OMIB[2]), #Connection between buses
     0.01, #resistance in pu
     0.05, #reactance in pu
-    (from = 0.0, to = 0.0), #susceptance in pu
+    (from = 0.001, to = 0.001), #susceptance in pu
     18.046, #rate in MW
     1.04,
 )]  #angle limits (-min and max)
@@ -74,7 +74,7 @@ branches_OMIB_fault(nodes_OMIB) = [PSY.Line(
     Arc(from = nodes_OMIB[1], to = nodes_OMIB[2]), #Connection between buses
     0.02, #resistance in pu
     0.1, #reactance in pu
-    (from = 0.0, to = 0.0), #susceptance in pu
+    (from = 0.001, to = 0.001), #susceptance in pu
     18.046, #rate in MW
     1.04,
 )]  #angle limits (-min and max)
diff --git a/test/data_tests/test01.jl b/test/data_tests/test01.jl
new file mode 100644
index 0000000..d3285a0
--- /dev/null
+++ b/test/data_tests/test01.jl
@@ -0,0 +1,39 @@
+using PowerSystems
+using NLsolve
+const PSY = PowerSystems
+
+include(joinpath(dirname(@__FILE__), "dynamic_test_data.jl"))
+include(joinpath(dirname(@__FILE__), "data_utils.jl"))
+############### Data Network ########################
+omib_file_dir = joinpath(dirname(@__FILE__), "OMIB.raw")
+omib_sys = System(PowerModelsData(omib_file_dir), runchecks = false)
+add_source_to_ref(omib_sys)
+res = solve_powerflow!(omib_sys, nlsolve)
+############### Data Dynamic devices ########################
+function dyn_gen_first_order(generator)
+    return PSY.DynamicGenerator(
+        1, #Number
+        "Case1Gen",
+        get_bus(generator), #bus
+        1.0, # ω_ref,
+        get_voltage(get_bus(generator)), #V_ref
+        get_activepower(generator), #P_ref
+        get_reactivepower(generator), #Q_ref
+        machine_OMIB(), #machine
+        shaft_damping(), #shaft
+        avr_none(), #avr
+        tg_none(), #tg
+        pss_none(),
+    ) #pss
+end
+
+#Attach dynamic generator. Currently use PSS/e format based on bus #.
+gen = [g for g in get_components(Generator, omib_sys)][1]
+case_gen = dyn_gen_first_order(gen)
+add_component!(omib_sys, case_gen)
+
+#Compute Y_bus after fault
+fault_branch = deepcopy(collect(get_components(Branch, omib_sys))[1])
+fault_branch.r = 0.02;
+fault_branch.x = 0.1;
+Ybus_fault = PSY.Ybus([fault_branch], get_components(Bus, omib_sys))[:, :]
diff --git a/test/data_tests/test02.jl b/test/data_tests/test02.jl
new file mode 100644
index 0000000..7d89e89
--- /dev/null
+++ b/test/data_tests/test02.jl
@@ -0,0 +1,40 @@
+using PowerSystems
+using NLsolve
+const PSY = PowerSystems
+
+############### Data Network ########################
+include(joinpath(dirname(@__FILE__), "dynamic_test_data.jl"))
+include(joinpath(dirname(@__FILE__), "data_utils.jl"))
+############### Data Network ########################
+threebus_file_dir = joinpath(dirname(@__FILE__), "ThreeBusNetwork.raw")
+threebus_sys = System(PowerModelsData(threebus_file_dir), runchecks = false)
+add_source_to_ref(threebus_sys)
+res = solve_powerflow!(threebus_sys, nlsolve)
+
+### Case 2 Generators ###
+
+function dyn_gen_second_order(generator)
+    return PSY.DynamicGenerator(
+        1, #Number
+        "Case2_$(get_name(generator))",
+        get_bus(generator), #bus
+        1.0, # ω_ref,
+        1.0, #V_ref
+        get_activepower(generator), #P_ref
+        get_reactivepower(generator), #Q_ref
+        machine_4th(), #machine
+        shaft_no_damping(), #shaft
+        avr_type1(), #avr
+        tg_none(), #tg
+        pss_none(),
+    ) #pss
+end
+
+for g in get_components(Generator, threebus_sys)
+    case_gen = dyn_gen_second_order(g)
+    add_component!(threebus_sys, case_gen)
+end
+
+#Compute Y_bus after fault
+fault_branches = deepcopy(collect(get_components(Branch, threebus_sys))[2:end])
+Ybus_fault = PSY.Ybus(fault_branches, get_components(Bus, threebus_sys))[:, :]
diff --git a/test/data_tests/test03.jl b/test/data_tests/test03.jl
new file mode 100644
index 0000000..3ba5462
--- /dev/null
+++ b/test/data_tests/test03.jl
@@ -0,0 +1,40 @@
+using PowerSystems
+using NLsolve
+const PSY = PowerSystems
+
+############### Data Network ########################
+include(joinpath(dirname(@__FILE__), "dynamic_test_data.jl"))
+include(joinpath(dirname(@__FILE__), "data_utils.jl"))
+############### Data Network ########################
+threebus_file_dir = joinpath(dirname(@__FILE__), "ThreeBusNetwork.raw")
+threebus_sys = System(PowerModelsData(threebus_file_dir), runchecks = false)
+add_source_to_ref(threebus_sys)
+res = solve_powerflow!(threebus_sys, nlsolve)
+
+### Case 2 Generators ###
+
+function dyn_gen_sixth_order(generator)
+    return PSY.DynamicGenerator(
+        1, #Number
+        "Case3_$(get_name(generator))",
+        get_bus(generator), #bus
+        1.0, # ω_ref,
+        get_voltage(get_bus(generator)), #V_ref
+        get_activepower(generator), #P_ref
+        get_reactivepower(generator), #Q_ref
+        machine_6th(), #machine
+        shaft_no_damping(), #shaft
+        avr_type1(), #avr
+        tg_none(), #tg
+        pss_none(),
+    ) #pss
+end
+
+for g in get_components(Generator, threebus_sys)
+    case_gen = dyn_gen_sixth_order(g)
+    add_component!(threebus_sys, case_gen)
+end
+
+#Compute Y_bus after fault
+fault_branches = deepcopy(collect(get_components(Branch, threebus_sys))[2:end])
+Ybus_fault = PSY.Ybus(fault_branches, get_components(Bus, threebus_sys))[:, :]
diff --git a/test/data_tests/test04.jl b/test/data_tests/test04.jl
new file mode 100644
index 0000000..e7445c6
--- /dev/null
+++ b/test/data_tests/test04.jl
@@ -0,0 +1,40 @@
+using PowerSystems
+using NLsolve
+const PSY = PowerSystems
+
+############### Data Network ########################
+include(joinpath(dirname(@__FILE__), "dynamic_test_data.jl"))
+include(joinpath(dirname(@__FILE__), "data_utils.jl"))
+############### Data Network ########################
+threebus_file_dir = joinpath(dirname(@__FILE__), "ThreeBusNetwork.raw")
+threebus_sys = System(PowerModelsData(threebus_file_dir), runchecks = false)
+add_source_to_ref(threebus_sys)
+res = solve_powerflow!(threebus_sys, nlsolve)
+
+### Case 2 Generators ###
+
+function dyn_gen_eight_order(generator)
+    return PSY.DynamicGenerator(
+        1, #Number
+        "Case4_$(get_name(generator))",
+        get_bus(generator), #bus
+        1.0, # ω_ref,
+        get_voltage(get_bus(generator)), #V_ref
+        get_activepower(generator), #P_ref
+        get_reactivepower(generator), #Q_ref
+        machine_8th(), #machine
+        shaft_no_damping(), #shaft
+        avr_type1(), #avr
+        tg_none(), #tg
+        pss_none(),
+    ) #pss
+end
+
+for g in get_components(Generator, threebus_sys)
+    case_gen = dyn_gen_eight_order(g)
+    add_component!(threebus_sys, case_gen)
+end
+
+#Compute Y_bus after fault
+fault_branches = deepcopy(collect(get_components(Branch, threebus_sys))[2:end])
+Ybus_fault = PSY.Ybus(fault_branches, get_components(Bus, threebus_sys))[:, :]
diff --git a/test/runtests.jl b/test/runtests.jl
index 13657a1..21109df 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,24 +4,69 @@ using Test
 using NLsolve
 using DiffEqBase
 using Sundials
+using InfrastructureSystems
+using Logging
 
+const IS = InfrastructureSystems
 const PSY = PowerSystems
 
-include("./data_tests/network_test_data.jl")
-include("./data_tests/dynamic_test_data.jl")
-
 tests = readdir(dirname(@__FILE__))
 tests = filter(
     f -> startswith(f, "test_") && endswith(f, ".jl") && f != basename(@__FILE__),
     tests,
 )
 
-@testset "BasicTests" begin
+const LOG_FILE = "lits-test.log"
+
+LOG_LEVELS = Dict(
+    "Debug" => Logging.Debug,
+    "Info" => Logging.Info,
+    "Warn" => Logging.Warn,
+    "Error" => Logging.Error,
+)
+
+function get_logging_level(env_name::String, default)
+    level = get(ENV, env_name, default)
+    log_level = get(LOG_LEVELS, level, nothing)
+    if isnothing(log_level)
+        error("Invalid log level $level: Supported levels: $(values(LOG_LEVELS))")
+    end
+
+    return log_level
+end
+
+function run_tests()
+    console_level = get_logging_level("SYS_CONSOLE_LOG_LEVEL", "Error")
+    console_logger = ConsoleLogger(stderr, console_level)
+    file_level = get_logging_level("SYS_LOG_LEVEL", "Info")
+
+    include("./data_tests/network_test_data.jl")
+    include("./data_tests/dynamic_test_data.jl")
 
-    for test in tests
-        print(splitext(test)[1], ": ")
-        include(test)
-        println()
+    IS.open_file_logger(LOG_FILE, file_level) do file_logger
+        multi_logger = IS.MultiLogger(
+            [console_logger, file_logger],
+            IS.LogEventTracker((Logging.Info, Logging.Warn, Logging.Error)),
+        )
+        global_logger(multi_logger)
+
+        @testset "BasicTests" begin
+            for test in tests
+                print(splitext(test)[1], ": ")
+                include(test)
+                println()
+            end
+        end
+        @info IS.report_log_summary(multi_logger)
     end
+end
+
+logger = global_logger()
 
+try
+    run_tests()
+finally
+    # Guarantee that the global logger is reset.
+    global_logger(logger)
+    nothing
 end
diff --git a/test/test_case01_OMIB.jl b/test/test_case01_OMIB.jl
index 8bae616..4a37f6f 100644
--- a/test/test_case01_OMIB.jl
+++ b/test/test_case01_OMIB.jl
@@ -8,39 +8,11 @@ drop a circuit on the (double circuit) line connecting the two buses, doubling i
 ############### LOAD DATA ########################
 ##################################################
 
-############### Data Network ########################
-
-nodes_case1 = nodes_OMIB()
-
-branch_case1 = branches_OMIB(nodes_case1)
-
-#Trip of a single circuit of Line 1 -> Resistance and Reactance doubled.
-branch_case1_fault = branches_OMIB_fault(nodes_case1)
-
-loads_case1 = loads_OMIB(nodes_case1)
-
-############### Data devices ########################
-
-inf_gen_case1 = inf_gen_105_pu(nodes_case1)
-
-### Case 1 Generator ###
-
-case1_gen = dyn_gen_OMIB(nodes_case1)
-
-######################### Dynamical System ########################
-
-#Create system with BasePower = 100 MVA and nominal frequency 60 Hz.
-sys = system_no_inv(nodes_case1, branch_case1, loads_case1, [inf_gen_case1], [case1_gen])
+include(joinpath(dirname(@__FILE__), "data_tests/test01.jl"))
 
 ##################################################
 ############### SOLVE PROBLEM ####################
 ##################################################
-
-#Compute Y_bus after fault
-Ybus_fault = get_admittance_matrix(nodes_case1, branch_case1_fault)
-
-tspan = (0.0, 30.0);
-
 #Define Fault: Change of YBus
 Ybus_change = ThreePhaseFault(
     1.0, #change at t = 1.0
@@ -49,8 +21,8 @@ Ybus_change = ThreePhaseFault(
 
 #Define Simulation Problem
 sim = Simulation(
-    sys, #system
-    tspan, #time span
+    omib_sys, #system
+    (0.0, 30.0), #time span
     Ybus_change,
 ) #Type of Fault
 
diff --git a/test/test_case02_4th_order.jl b/test/test_case02_4th_order.jl
index 757dc3b..47414bc 100644
--- a/test/test_case02_4th_order.jl
+++ b/test/test_case02_4th_order.jl
@@ -9,48 +9,12 @@ and the generator located in bus 3.
 ############### LOAD DATA ########################
 ##################################################
 
-############### Data Network ########################
-
-nodes_case234 = nodes_3bus()
-
-branch_case234 = branches_3lines(nodes_case234)
-
-#Trip of Line 1.
-branch_case234_fault = branches_3lines_fault(nodes_case234)
-
-loads_case234 = loads_3bus(nodes_case234)
-
-############### Data devices ########################
-
-inf_gen_case234 = inf_gen_102_pu(nodes_case234)
-
-### Case 2 Generators ###
-
-case2_gen2 = dyn_gen2_case2(nodes_case234)
-
-case2_gen3 = dyn_gen3_case2(nodes_case234)
-
-######################### Dynamical System ########################
-
-#Create system with BasePower = 100 MVA and nominal frequency 60 Hz.
-sys = system_no_inv(
-    nodes_case234,
-    branch_case234,
-    loads_case234,
-    [inf_gen_case234],
-    [case2_gen2, case2_gen3],
-)
+include(joinpath(dirname(@__FILE__), "data_tests/test02.jl"))
 
 ##################################################
 ############### SOLVE PROBLEM ####################
 ##################################################
 
-#Compute Y_bus after fault
-Ybus_fault = get_admittance_matrix(nodes_case234, branch_case234_fault)
-
-#time span
-tspan = (0.0, 30.0);
-
 #Initial guess
 x0_guess = [
     1.02,
@@ -85,16 +49,16 @@ Ybus_change = ThreePhaseFault(
 
 #Define Simulation Problem
 sim = Simulation(
-    sys, #system
-    tspan, #time span
+    threebus_sys, #system
+    (0.0, 30.0), #time span
     Ybus_change, #Type of Fault
     initial_guess = x0_guess,
 ) #initial guess
 
 #Solve problem in equilibrium
-run_simulation!(sim, IDA());
+run_simulation!(sim, IDA())
 
 #Obtain data for angles
-series = get_state_series(sim, ("Case2Gen2", :δ));
+series = get_state_series(sim, ("Case2_generator-2-1", :δ))
 
 @test sim.solution.retcode == :Success
diff --git a/test/test_case03_6th_order.jl b/test/test_case03_6th_order.jl
index 839630b..9e46cd1 100644
--- a/test/test_case03_6th_order.jl
+++ b/test/test_case03_6th_order.jl
@@ -9,45 +9,12 @@ and the generator located in bus 3.
 ############### LOAD DATA ########################
 ##################################################
 
-############### Data Network ########################
-
-nodes_case234 = nodes_3bus()
-
-branch_case234 = branches_3lines(nodes_case234)
-
-#Trip of Line 1.
-branch_case234_fault = branches_3lines_fault(nodes_case234)
-
-loads_case234 = loads_3bus(nodes_case234)
-
-############### Data devices ########################
-
-inf_gen_case234 = inf_gen_102_pu(nodes_case234)
-
-### Case 3 Generators ###
-
-case3_gen2 = dyn_gen2_case3(nodes_case234)
-
-case3_gen3 = dyn_gen3_case3(nodes_case234)
-
-######################### Dynamical System ########################
-
-#Create system with BasePower = 100 MVA and nominal frequency 60 Hz.
-sys = system_no_inv(
-    nodes_case234,
-    branch_case234,
-    loads_case234,
-    [inf_gen_case234],
-    [case3_gen2, case3_gen3],
-)
+include(joinpath(dirname(@__FILE__), "data_tests/test03.jl"))
 
 ##################################################
 ############### SOLVE PROBLEM ####################
 ##################################################
 
-#Compute Y_bus after fault
-Ybus_fault = get_admittance_matrix(nodes_case234, branch_case234_fault)
-
 #time span
 tspan = (0.0, 20.0);
 
@@ -89,7 +56,7 @@ Ybus_change = ThreePhaseFault(
 
 #Define Simulation Problem
 sim = Simulation(
-    sys, #system
+    threebus_sys, #system
     tspan, #time span
     Ybus_change, #Type of Fault
     initial_guess = x0_guess,
@@ -99,6 +66,6 @@ sim = Simulation(
 run_simulation!(sim, IDA());
 
 #Obtain data for angles
-series = get_state_series(sim, ("Case3Gen2", :δ));
+series = get_state_series(sim, ("Case3_generator-2-1", :δ));
 
 @test sim.solution.retcode == :Success
diff --git a/test/test_case04_8th_order.jl b/test/test_case04_8th_order.jl
index 90b4c7d..e0868fb 100644
--- a/test/test_case04_8th_order.jl
+++ b/test/test_case04_8th_order.jl
@@ -9,48 +9,11 @@ and the generator located in bus 3.
 ############### LOAD DATA ########################
 ##################################################
 
-############### Data Network ########################
-
-nodes_case234 = nodes_3bus()
-
-branch_case234 = branches_3lines(nodes_case234)
-
-#Trip of Line 1.
-branch_case234_fault = branches_3lines_fault(nodes_case234)
-
-loads_case234 = loads_3bus(nodes_case234)
-
-############### Data devices ########################
-
-inf_gen_case234 = inf_gen_102_pu(nodes_case234)
-
-### Case 4 Generators ###
-
-case4_gen2 = dyn_gen2_case4(nodes_case234)
-
-case4_gen3 = dyn_gen3_case4(nodes_case234)
-
-######################### Dynamical System ########################
-
-#Create system with BasePower = 100 MVA and nominal frequency 60 Hz.
-sys = system_no_inv(
-    nodes_case234,
-    branch_case234,
-    loads_case234,
-    [inf_gen_case234],
-    [case4_gen2, case4_gen3],
-)
+include(joinpath(dirname(@__FILE__), "data_tests/test04.jl"))
 
 ##################################################
 ############### SOLVE PROBLEM ####################
 ##################################################
-
-#Compute Y_bus after fault
-Ybus_fault = get_admittance_matrix(nodes_case234, branch_case234_fault)
-
-#time span
-tspan = (0.0, 20.0);
-
 #Initial guess
 x0_guess = [
     1.02,
@@ -93,8 +56,8 @@ Ybus_change = ThreePhaseFault(
 
 #Define Simulation Problem
 sim = Simulation(
-    sys, #system
-    tspan, #time span
+    threebus_sys, #system
+    (0.0, 20.0), #time span
     Ybus_change, #Type of Fault
     initial_guess = x0_guess,
 ) #initial guess
@@ -103,6 +66,6 @@ sim = Simulation(
 run_simulation!(sim, IDA());
 
 #Obtain data for angles
-series = get_state_series(sim, ("Case4Gen2", :δ));
+series = get_state_series(sim, ("Case4_generator-2-1", :δ));
 
 @test sim.solution.retcode == :Success