diff --git a/Project.toml b/Project.toml
index e1d4560..09af653 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "LITS"
 uuid = "86b0dc02-7903-11e9-325f-f195ca7e6c1a"
-version = "0.3.1"
+version = "0.4.0"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
@@ -11,11 +11,10 @@ NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
-
 [compat]
 DiffEqBase = "6"
 ForwardDiff = "~v0.10"
-InfrastructureSystems = "~0.5"
+InfrastructureSystems = "~0.6"
 NLsolve = "4"
-PowerSystems = "~0.11"
+PowerSystems = "~0.14"
 julia = "^1.2"
diff --git a/dev_files/init_cond/Twogen2bus.m b/dev_files/init_cond/Twogen2bus.m
index 03a0908..f940ecc 100644
--- a/dev_files/init_cond/Twogen2bus.m
+++ b/dev_files/init_cond/Twogen2bus.m
@@ -15,7 +15,7 @@
 %% generator data
 %	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin	Pc1	Pc2	Qc1min	Qc1max	Qc2min	Qc2max	ramp_agc	ramp_10	ramp_30	ramp_q	apf
 mpc.gen = [
-	1	700	70	100	-100	1.02	1000	1000	1000	0	0	0	0	0	0	0	0	0	0	0;
+	1	700	70	100	-100	1.02	1000	1	1000	0	0	0	0	0	0	0	0	0	0	0;
 	%2	50 -30	100	-100	1.0	100	1	150	37.5	0	0	0	0	0	0	0	0	0	0	0;
 ];
 
diff --git a/dev_files/init_cond/initial_conditions.ipynb b/dev_files/init_cond/initial_conditions.ipynb
index 40dc4e3..854480f 100644
--- a/dev_files/init_cond/initial_conditions.ipynb
+++ b/dev_files/init_cond/initial_conditions.ipynb
@@ -2,37 +2,156 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 8,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "PowerSystems"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
    "source": [
+    "using Revise\n",
     "using PowerModels\n",
     "using Ipopt\n",
     "using ForwardDiff\n",
-    "using OrdinaryDiffEq\n",
-    "using NLsolve"
+    "#using OrdinaryDiffEq\n",
+    "using NLsolve\n",
+    "using PowerSystems\n",
+    "const PSY = PowerSystems"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 9,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<style>.container { width:100% !important; }</style>"
+      ],
+      "text/plain": [
+       "HTML{String}(\"<style>.container { width:100% !important; }</style>\")"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
    "source": [
     "display(HTML(\"<style>.container { width:100% !important; }</style>\"))"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": false
-   },
-   "outputs": [],
+   "execution_count": 21,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "┌ Info: Constructing System from Power Models\n",
+      "│   data[\"name\"] = Twogen2bus\n",
+      "│   data[\"source_type\"] = matpower\n",
+      "└ @ PowerSystems /Users/jdlara/.julia/packages/PowerSystems/wHSk9/src/parsers/power_models_data.jl:33\n",
+      "┌ Info: Reading bus data\n",
+      "└ @ PowerSystems /Users/jdlara/.julia/packages/PowerSystems/wHSk9/src/parsers/power_models_data.jl:122\n",
+      "┌ Info: Reading generator data\n",
+      "└ @ PowerSystems /Users/jdlara/.julia/packages/PowerSystems/wHSk9/src/parsers/power_models_data.jl:368\n",
+      "┌ Warning: Generator cost data not included for Generator: gen-1-1\n",
+      "└ @ PowerSystems /Users/jdlara/.julia/packages/PowerSystems/wHSk9/src/parsers/power_models_data.jl:330\n",
+      "┌ Info: Reading branch data\n",
+      "└ @ PowerSystems /Users/jdlara/.julia/packages/PowerSystems/wHSk9/src/parsers/power_models_data.jl:520\n",
+      "┌ Warning: Rate provided for 1 is 32.0, larger the SIL 258.94053053874836 in the range of (min = 325.0, max = 425.0).\n",
+      "└ @ PowerSystems /Users/jdlara/.julia/packages/PowerSystems/wHSk9/src/utils/IO/branchdata_checks.jl:150\n",
+      "┌ Info: Reading branch data\n",
+      "└ @ PowerSystems /Users/jdlara/.julia/packages/PowerSystems/wHSk9/src/parsers/power_models_data.jl:583\n",
+      "┌ Info: Reading DC Line data\n",
+      "└ @ PowerSystems /Users/jdlara/.julia/packages/PowerSystems/wHSk9/src/parsers/power_models_data.jl:554\n",
+      "┌ Info: PowerFlow solve converged, the results have been stored in the system\n",
+      "└ @ PowerSystems /Users/jdlara/.julia/packages/PowerSystems/wHSk9/src/utils/power_flow/power_flow.jl:147\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "true"
+      ]
+     },
+     "execution_count": 21,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "file_dir = joinpath(pwd(), \"Twogen2bus.m\")\n",
+    "sys = System(PSY.PowerModelsData(file_dir))\n",
+    "solve_powerflow!(sys, nlsolve, method = :newton)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 26,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Vm = 1.02 Va= 0.0\n",
+      "Vm = 1.0 Va= -0.1279338657728345\n"
+     ]
+    }
+   ],
+   "source": [
+    "for b in get_components(Bus, sys)\n",
+    "    println(\"Vm = $(get_voltage(b))\", \" \",\"Va= $(get_angle(b))\")\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 19,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\u001b[1mMetadata\u001b[0m\n",
+      "  baseMVA: 1000\n",
+      "  per_unit: true\n",
+      "\n",
+      "\u001b[1mTable Counts\u001b[0m\n",
+      "  bus: 2\n",
+      "  gen: 1\n",
+      "\n",
+      "\n",
+      "\u001b[1mTable: bus\u001b[0m\n",
+      "            vm,     va\n",
+      "      1: 1.020,  0.000\n",
+      "      2: 1.000, -0.128\n",
+      "\n",
+      "\n",
+      "\u001b[1mTable: gen\u001b[0m\n",
+      "            pg,    qg\n",
+      "      1: 0.703, 0.051\n"
+     ]
+    }
+   ],
    "source": [
-    "data = PowerModels.parse_file(\"C:/Users/Ewa/Downloads/Twogen2bus.m\")\n",
-    "pf_result = run_ac_pf(data ,Ipopt.Optimizer)\n",
-    "#PowerModels.print_summary(pf_result[\"solution\"])"
+    "data = PowerModels.parse_file(file_dir)\n",
+    "opt = optimizer_with_attributes(Ipopt.Optimizer, \"print_level\" => 0)\n",
+    "pf_result = run_ac_pf(data, opt)\n",
+    "PowerModels.print_summary(pf_result[\"solution\"])"
    ]
   },
   {
@@ -301,17 +420,17 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Julia 1.3.1",
+   "display_name": "Julia 1.4.0",
    "language": "julia",
-   "name": "julia-1.3"
+   "name": "julia-1.4"
   },
   "language_info": {
    "file_extension": ".jl",
    "mimetype": "application/julia",
    "name": "julia",
-   "version": "1.3.1"
+   "version": "1.4.0"
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }
diff --git a/docs/src/Examples/example_lines.md b/docs/src/Examples/example_lines.md
index 492a86b..14241e5 100644
--- a/docs/src/Examples/example_lines.md
+++ b/docs/src/Examples/example_lines.md
@@ -368,8 +368,8 @@ x0_guess = [
     0.0, #V1_I
     -0.01, #V2_I
     -0.01, #V3_I
-    0.0, #δω_vsm
-    0.2, #δθ_vsm
+    0.0, #ω_oc
+    0.2, #θ_oc
     0.025, #qm
     0.0015, #ξ_d
     -0.07, #ξ_q
@@ -380,13 +380,13 @@ x0_guess = [
     1.004, #vpll_d
     0.0, #vpll_q
     0.0, #ε_pll
-    0.1, #δθ_pll
+    0.1, #θ_pll
     0.5, #id_cv
     0.0, #iq_cv
-    0.95, #vod
-    -0.1, #voq
-    0.49, #iod
-    -0.1, #ioq
+    0.95, #Vd_filter
+    -0.1, #Vq_filter
+    0.49, #Id_filter
+    -0.1, #Iq_filter
     1.0, #eq_p
     0.47, #ed_p
     0.6, #δ
diff --git a/src/LITS.jl b/src/LITS.jl
index a281f03..1e604d5 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")
 
@@ -57,7 +58,7 @@ include("models/inverter_models/DCside_models.jl")
 include("models/inverter_models/filter_models.jl")
 include("models/inverter_models/frequency_estimator_models.jl")
 include("models/inverter_models/outer_control_models.jl")
-include("models/inverter_models/voltage_source_control_models.jl")
+include("models/inverter_models/inner_control_models.jl")
 include("models/inverter_models/converter_models.jl")
 
 #Injection Models
@@ -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..469c857 100644
--- a/src/base/definitions.jl
+++ b/src/base/definitions.jl
@@ -40,32 +40,32 @@ 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
-Vdo_var :: Voltage seen in the capacitor of the filter in the d-component
-Vqo_var :: Voltage seen in the capacitor of the filter in the q-component
+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
 ω_freq_estimator_var :: Frequency estimated by the frequency estimator.
-v_control_var :: Control voltage supplied from the outer loop control to the inner loop
-ω_control_var :: Control frequency supplied from the outer loop control the inner loop
-δdqRI_var :: Variation of the angle (PLL or VSM) of the inverter
+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
-Vdcnv_var :: Voltage supplied from the converter in the d-component
-Vqcnv_var :: Voltage supplied from the converter in the q-component
+Vd_cnv_var :: Voltage supplied from the converter in the d-component
+Vq_cnv_var :: Voltage supplied from the converter in the q-component
 """
 
 @enum inverter_inner_vars begin
     md_var = 1
     mq_var = 2
     Vdc_var = 3
-    Vdo_var = 4
-    Vqo_var = 5
+    Vd_filter_var = 4
+    Vq_filter_var = 5
     ω_freq_estimator_var = 6
-    v_control_var = 7
-    ω_control_var = 8
-    δdqRI_var = 9
+    V_oc_var = 7
+    ω_oc_var = 8
+    θ_oc_var = 9
     VR_inv_var = 10
     VI_inv_var = 11
-    Vdcnv_var = 12
-    Vqcnv_var = 13
+    Vd_cnv_var = 12
+    Vq_cnv_var = 13
 end
 
 Base.to_index(ix::inverter_inner_vars) = Int(ix)
@@ -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/ports.jl b/src/base/ports.jl
index 7b0a5e1..a10f180 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, Vdcnv_var, Vqcnv_var]
+    inner_input = [md_var, mq_var, Vdc_var, Vd_cnv_var, Vq_cnv_var]
     return Ports(state_input, inner_input)
 end
 
@@ -59,31 +59,38 @@ end
 #### Filter Ports ####
 function Ports(::PSY.Filter)
     #TODO: If converter has dynamics, need to connect state_input
-    state_input = [:δθ_vsm] #[:Vd_c, :Vq_c] #, :Id_c, :Iq_c]
-    inner_input =
-        [VR_inv_var, VI_inv_var, Vdcnv_var, Vqcnv_var, δdqRI_var, Vdo_var, Vqo_var]
+    state_input = [:θ_oc] #[:Vd_c, :Vq_c] #, :Id_c, :Iq_c]
+    inner_input = [
+        VR_inv_var,
+        VI_inv_var,
+        Vd_cnv_var,
+        Vq_cnv_var,
+        θ_oc_var,
+        Vd_filter_var,
+        Vq_filter_var,
+    ]
     return Ports(state_input, inner_input)
 end
 
 #### Freq. Estimator Ports ####
 
 function Ports(::PSY.FrequencyEstimator)
-    state_input = [:vd_cap, :vq_cap, :δθ_vsm]
+    state_input = [:vd_filter, :vq_filter, :θ_oc]
     #TODO: Move PLL to PCC, i.e. move v_cap (D'Arco v_o), to inner inputs
-    inner_input = [Vdo_var, Vqo_var, δdqRI_var, ω_freq_estimator_var]
+    inner_input = [Vd_filter_var, Vq_filter_var, θ_oc_var, ω_freq_estimator_var]
     return Ports(state_input, inner_input)
 end
 
 #### Outer Control Ports ####
 function Ports(::PSY.OuterControl)
-    state_input = [:vpll_d, :vpll_q, :ε_pll, :vd_cap, :vq_cap, :id_o, :iq_o]
-    inner_input = [Vdo_var, Vqo_var, ω_freq_estimator_var]
+    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]
     return Ports(state_input, inner_input)
 end
 
 #### Inner Control Ports ####
-function Ports(::PSY.VSControl)
-    state_input = [:id_o, :iq_o, :id_c, :iq_c, :vd_cap, :vq_cap]
-    inner_input = [Vdo_var, Vqo_var, v_control_var, ω_control_var]
+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]
     return Ports(state_input, inner_input)
 end
diff --git a/src/base/simulation.jl b/src/base/simulation.jl
index 56f8cda..76ece36 100644
--- a/src/base/simulation.jl
+++ b/src/base/simulation.jl
@@ -134,7 +134,7 @@ end
 
 function _attach_inner_vars!(
     device::PSY.DynamicGenerator,
-    ::Type{T} = Float64,
+    ::Type{T} = Real,
 ) where {T <: Real}
     device.ext[INNER_VARS] = zeros(T, 8)
     return
@@ -142,7 +142,7 @@ end
 
 function _attach_inner_vars!(
     device::PSY.DynamicInverter,
-    ::Type{T} = Float64,
+    ::Type{T} = Real,
 ) where {T <: Real}
     device.ext[INNER_VARS] = zeros(T, 13)
     return
@@ -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/device.jl b/src/models/device.jl
index 96783fe..eaa9a48 100644
--- a/src/models/device.jl
+++ b/src/models/device.jl
@@ -106,14 +106,11 @@ function device!(
 
     #Update V_ref
     V_ref = PSY.get_ext(device)[CONTROL_REFS][V_ref_index]
-    get_inner_vars(device)[v_control_var] = V_ref
+    get_inner_vars(device)[V_oc_var] = V_ref
 
     #Obtain ODES for DC side
     mdl_DCside_ode!(device)
 
-    #Obtain ODEs for OuterLoop
-    mdl_outer_ode!(device_states, view(output_ode, ode_range), sys_f0, sys_ω, device)
-
     #Obtain ODEs for PLL
     mdl_freq_estimator_ode!(
         device_states,
@@ -123,8 +120,11 @@ function device!(
         device,
     )
 
+    #Obtain ODEs for OuterLoop
+    mdl_outer_ode!(device_states, view(output_ode, ode_range), sys_f0, sys_ω, device)
+
     #Obtain inner controller ODEs and modulation commands
-    mdl_VScontrol_ode!(device_states, view(output_ode, ode_range), device)
+    mdl_inner_ode!(device_states, view(output_ode, ode_range), device)
 
     #Obtain converter relations
     mdl_converter_ode!(device)
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/inverter_models/DCside_models.jl b/src/models/inverter_models/DCside_models.jl
index d0e0aae..49cae94 100644
--- a/src/models/inverter_models/DCside_models.jl
+++ b/src/models/inverter_models/DCside_models.jl
@@ -1,9 +1,9 @@
 function mdl_DCside_ode!(
-    device::PSY.DynamicInverter{C, O, VC, PSY.FixedDCSource, P, F},
+    device::PSY.DynamicInverter{C, O, IC, PSY.FixedDCSource, P, F},
 ) where {
     C <: PSY.Converter,
     O <: PSY.OuterControl,
-    VC <: PSY.VSControl,
+    IC <: PSY.InnerControl,
     P <: PSY.FrequencyEstimator,
     F <: PSY.Filter,
 }
diff --git a/src/models/inverter_models/converter_models.jl b/src/models/inverter_models/converter_models.jl
index 31e8420..4aebba6 100644
--- a/src/models/inverter_models/converter_models.jl
+++ b/src/models/inverter_models/converter_models.jl
@@ -1,8 +1,8 @@
 function mdl_converter_ode!(
-    device::PSY.DynamicInverter{PSY.AvgCnvFixedDC, O, VC, DC, P, F},
+    device::PSY.DynamicInverter{PSY.AverageConverter, O, IC, DC, P, F},
 ) where {
     O <: PSY.OuterControl,
-    VC <: PSY.VSControl,
+    IC <: PSY.InnerControl,
     DC <: PSY.DCSource,
     P <: PSY.FrequencyEstimator,
     F <: PSY.Filter,
@@ -16,8 +16,8 @@ function mdl_converter_ode!(
     VDC = get_inner_vars(device)[Vdc_var]
 
     #Update inner_vars
-    get_inner_vars(device)[Vdcnv_var] = md * VDC
-    get_inner_vars(device)[Vqcnv_var] = mq * VDC
+    get_inner_vars(device)[Vd_cnv_var] = md * VDC
+    get_inner_vars(device)[Vq_cnv_var] = mq * VDC
 end
 
 #TODO: Same as above, but:
diff --git a/src/models/inverter_models/filter_models.jl b/src/models/inverter_models/filter_models.jl
index 3c92a4b..eb21d51 100644
--- a/src/models/inverter_models/filter_models.jl
+++ b/src/models/inverter_models/filter_models.jl
@@ -6,11 +6,11 @@ function mdl_filter_ode!(
     sys_Sbase,
     f0,
     ω_sys,
-    device::PSY.DynamicInverter{C, O, VC, DC, P, PSY.LCLFilter},
+    device::PSY.DynamicInverter{C, O, IC, DC, P, PSY.LCLFilter},
 ) where {
     C <: PSY.Converter,
     O <: PSY.OuterControl,
-    VC <: PSY.VSControl,
+    IC <: PSY.InnerControl,
     DC <: PSY.DCSource,
     P <: PSY.FrequencyEstimator,
 }
@@ -18,16 +18,16 @@ function mdl_filter_ode!(
     #Obtain external states inputs for component
     #TODO: If converter has dynamics, need to reference states:
     #external_ix = device.input_port_mapping[device.converter]
-    #vcvd = device_states[external_ix[1]]
-    #vcvq = device_states[external_ix[2]]
+    #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]]
 
     #Obtain inner variables for component
     V_tR = get_inner_vars(device)[VR_inv_var]
     V_tI = get_inner_vars(device)[VI_inv_var]
-    vcvd = get_inner_vars(device)[Vdcnv_var]
-    vcvq = get_inner_vars(device)[Vqcnv_var]
+    Vd_cnv = get_inner_vars(device)[Vd_cnv_var]
+    Vq_cnv = get_inner_vars(device)[Vq_cnv_var]
 
     #Get parameters
     filter = PSY.get_filter(device)
@@ -48,45 +48,55 @@ function mdl_filter_ode!(
 
     #Define internal states for filter
     internal_states = @view device_states[local_ix]
-    icvd = internal_states[1]
-    icvq = internal_states[2]
-    vod = internal_states[3]
-    voq = internal_states[4]
-    iod = internal_states[5]
-    ioq = internal_states[6]
+    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]
 
     #Inputs (control signals) - N/A
 
     #Compute 6 states ODEs (D'Arco EPSR122 Model)
     #Inverter Output Inductor (internal state)
     #𝜕id_c/𝜕t
-    output_ode[local_ix[1]] =
-        (ωb / lf * vcvd - ωb / lf * vod - ωb * rf / lf * icvd + ωb * ω_sys * icvq)
+    output_ode[local_ix[1]] = (
+        ωb / lf * Vd_cnv - ωb / lf * Vd_filter - ωb * rf / lf * Id_cnv +
+        ωb * ω_sys * Iq_cnv
+    )
     #𝜕iq_c/𝜕t
-    output_ode[local_ix[2]] =
-        (ωb / lf * vcvq - ωb / lf * voq - ωb * rf / lf * icvq - ωb * ω_sys * icvd)
+    output_ode[local_ix[2]] = (
+        ωb / lf * Vq_cnv - ωb / lf * Vq_filter - ωb * rf / lf * Iq_cnv -
+        ωb * ω_sys * Id_cnv
+    )
     #LCL Capacitor (internal state)
     #𝜕vd_o/𝜕t
-    output_ode[local_ix[3]] = (ωb / cf * icvd - ωb / cf * iod + ωb * ω_sys * voq)
+    output_ode[local_ix[3]] =
+        (ωb / cf * Id_cnv - ωb / cf * Id_filter + ωb * ω_sys * Vq_filter)
     #𝜕vq_o/𝜕t
-    output_ode[local_ix[4]] = (ωb / cf * icvq - ωb / cf * ioq - ωb * ω_sys * vod)
+    output_ode[local_ix[4]] =
+        (ωb / cf * Iq_cnv - ωb / cf * Iq_filter - ωb * ω_sys * Vd_filter)
     #Grid Inductance (internal state)
     #𝜕id_o/𝜕t
-    output_ode[local_ix[5]] =
-        (ωb / lg * vod - ωb / lg * V_dq[2] - ωb * rg / lg * iod + ωb * ω_sys * ioq)
+    output_ode[local_ix[5]] = (
+        ωb / lg * Vd_filter - ωb / lg * V_dq[2] - ωb * rg / lg * Id_filter +
+        ωb * ω_sys * Iq_filter
+    )
     #𝜕iq_o/𝜕t
-    output_ode[local_ix[6]] =
-        (ωb / lg * voq + ωb / lg * V_dq[1] - ωb * rg / lg * ioq - ωb * ω_sys * iod)
+    output_ode[local_ix[6]] = (
+        ωb / lg * Vq_filter + ωb / lg * V_dq[1] - ωb * rg / lg * Iq_filter -
+        ωb * ω_sys * Id_filter
+    )
 
     #Update inner_vars
-    get_inner_vars(device)[Vdo_var] = vod
-    get_inner_vars(device)[Vqo_var] = voq
+    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)[Vdo_var] = V_dq[q::dq_ref]
-    #get_inner_vars(device)[Vqo_var] = V_dq[q::dq_ref]
+    #get_inner_vars(device)[Vd_filter_var] = V_dq[q::dq_ref]
+    #get_inner_vars(device)[Vq_filter_var] = V_dq[q::dq_ref]
 
-    #Compute current from the generator to the grid
-    I_RI = (MVABase / sys_Sbase) * dq_ri(δ) * [iod; ioq]
+    #Compute current from the inverter to the grid
+    I_RI = (MVABase / sys_Sbase) * dq_ri(δ) * [Id_filter; Iq_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 1c12246..e943424 100644
--- a/src/models/inverter_models/frequency_estimator_models.jl
+++ b/src/models/inverter_models/frequency_estimator_models.jl
@@ -3,25 +3,25 @@ function mdl_freq_estimator_ode!(
     output_ode,
     f0,
     ω_sys,
-    device::PSY.DynamicInverter{C, O, VC, DC, PSY.PLL, F},
+    device::PSY.DynamicInverter{C, O, IC, DC, PSY.KauraPLL, F},
 ) where {
     C <: PSY.Converter,
     O <: PSY.OuterControl,
-    VC <: PSY.VSControl,
+    IC <: PSY.InnerControl,
     DC <: PSY.DCSource,
     F <: PSY.Filter,
 }
 
     #Obtain external states inputs for component
-    external_ix = get_input_port_ix(device, PSY.PLL)
-    vod = device_states[external_ix[1]]
-    voq = device_states[external_ix[2]]
-    δθ_vsm = device_states[external_ix[3]]
+    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]]
 
     #Obtain inner variables for component
-    #vod = device.inner_vars[Vdo_var]
-    #voq = device.inner_vars[Vqo_var]
-    #δθ_vsm = device.inner_vars[δdqRI_var]
+    #Vd_filter = device.inner_vars[Vd_filter_var]
+    #Vq_filter = device.inner_vars[Vq_filter_var]
+    #θ_oc = device.inner_vars[θ_oc_var]
 
     #Get parameters
     pll_control = PSY.get_freq_estimator(device)
@@ -31,14 +31,14 @@ function mdl_freq_estimator_ode!(
     ωb = 2.0 * pi * f0
 
     #Obtain indices for component w/r to device
-    local_ix = get_local_state_ix(device, PSY.PLL)
+    local_ix = get_local_state_ix(device, PSY.KauraPLL)
 
     #Define internal states for frequency estimator
     internal_states = @view device_states[local_ix]
     vpll_d = internal_states[1]
     vpll_q = internal_states[2]
     ϵ_pll = internal_states[3]
-    δθ_pll = internal_states[4]
+    θ_pll = internal_states[4]
 
     #Inputs (control signals)
 
@@ -46,19 +46,19 @@ function mdl_freq_estimator_ode!(
     #Output Voltage LPF (internal state)
     #𝜕vpll_d/𝜕t, D'Arco ESPR122 eqn. 12
     output_ode[local_ix[1]] = (
-        ω_lp * vod * cos(δθ_pll - δθ_vsm) + ω_lp * voq * sin(δθ_pll - δθ_vsm) -
+        ω_lp * Vd_filter * cos(θ_pll - θ_oc) + ω_lp * Vq_filter * sin(θ_pll - θ_oc) -
         ω_lp * vpll_d
     )
     #𝜕vpll_q/𝜕t, D'Arco ESPR122 eqn. 12
     output_ode[local_ix[2]] = (
-        -ω_lp * vod * sin(δθ_pll - δθ_vsm) + ω_lp * voq * cos(δθ_pll - δθ_vsm) -
+        -ω_lp * Vd_filter * sin(θ_pll - θ_oc) + ω_lp * Vq_filter * cos(θ_pll - θ_oc) -
         ω_lp * vpll_q
     )
     #PI Integrator (internal state)
     #𝜕dϵ_pll/𝜕t, D'Arco ESPR122 eqn. 13
     output_ode[local_ix[3]] = atan(vpll_q / vpll_d)
     #PLL Frequency Deviation (internal state)
-    #𝜕δθ_pll/𝜕t, D'Arco ESPR122 eqn. 15
+    #𝜕θ_pll/𝜕t, D'Arco ESPR122 eqn. 15
     output_ode[local_ix[4]] = (ωb * kp_pll * atan(vpll_q / vpll_d) + ωb * ki_pll * ϵ_pll)
 
     #Update inner_vars
diff --git a/src/models/inverter_models/inner_control_models.jl b/src/models/inverter_models/inner_control_models.jl
new file mode 100644
index 0000000..e7a1384
--- /dev/null
+++ b/src/models/inverter_models/inner_control_models.jl
@@ -0,0 +1,110 @@
+function mdl_inner_ode!(
+    device_states,
+    output_ode,
+    device::PSY.DynamicInverter{C, O, PSY.CurrentControl, DC, P, F},
+) where {
+    C <: PSY.Converter,
+    O <: PSY.OuterControl,
+    DC <: PSY.DCSource,
+    P <: PSY.FrequencyEstimator,
+    F <: PSY.Filter,
+}
+
+    #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
+
+    #Obtain inner variables for component
+    #Vd_filter = get_inner_vars(device)[Vd_filter_var]
+    #Vq_filter = get_inner_vars(device)[Vq_filter_var]
+    ω_oc = get_inner_vars(device)[ω_oc_var]
+    v_refr = get_inner_vars(device)[V_oc_var]
+    vdc = get_inner_vars(device)[Vdc_var]
+
+    #Get Voltage Controller parameters
+    inner_control = PSY.get_inner_control(device)
+    filter = PSY.get_filter(device)
+    kpv = PSY.get_kpv(inner_control)
+    kiv = PSY.get_kiv(inner_control)
+    kffi = PSY.get_kffi(inner_control)
+    cf = PSY.get_cf(filter)
+    rv = PSY.get_rv(inner_control)
+    lv = PSY.get_lv(inner_control)
+
+    #Get Current Controller parameters
+    kpc = PSY.get_kpc(inner_control)
+    kic = PSY.get_kic(inner_control)
+    kffv = PSY.get_kffv(inner_control)
+    lf = PSY.get_lf(filter)
+    ωad = PSY.get_ωad(inner_control)
+    kad = PSY.get_kad(inner_control)
+
+    #Obtain indices for component w/r to device
+    local_ix = get_local_state_ix(device, PSY.CurrentControl)
+
+    #Define internal states for frequency estimator
+    internal_states = @view device_states[local_ix]
+    ξ_d = internal_states[1]
+    ξ_q = internal_states[2]
+    γ_d = internal_states[3]
+    γ_q = internal_states[4]
+    ϕ_d = internal_states[5]
+    ϕ_q = internal_states[6]
+
+    #Inputs (control signals)
+
+    ### 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)
+    #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
+    )
+
+    Iq_cnv_ref = (
+        kpv * (Vq_filter_ref - Vq_filter) +
+        kiv * ξ_q +
+        cf * ω_oc * Vd_filter +
+        kffi * Iq_filter
+    )
+    #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)
+    #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)
+    )
+    Vq_cnv_ref = (
+        kpc * (Iq_cnv_ref - Iq_cnv) + kic * γ_q + ω_oc * lf * Id_cnv + kffv * Vq_filter - kad * (Vq_filter - ϕ_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
+
+    #Update inner_vars
+    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 19b6104..ede1b74 100644
--- a/src/models/inverter_models/outer_control_models.jl
+++ b/src/models/inverter_models/outer_control_models.jl
@@ -5,15 +5,15 @@ function mdl_outer_ode!(
     ω_sys,
     device::PSY.DynamicInverter{
         C,
-        PSY.VirtualInertiaQdroop{PSY.VirtualInertia, PSY.ReactivePowerDroop},
-        VC,
+        PSY.OuterControl{PSY.VirtualInertia, PSY.ReactivePowerDroop},
+        IC,
         DC,
         P,
         F,
     },
 ) where {
     C <: PSY.Converter,
-    VC <: PSY.VSControl,
+    IC <: PSY.InnerControl,
     DC <: PSY.DCSource,
     P <: PSY.FrequencyEstimator,
     F <: PSY.Filter,
@@ -22,21 +22,21 @@ function mdl_outer_ode!(
     #Obtain external states inputs for component
     external_ix = get_input_port_ix(
         device,
-        PSY.VirtualInertiaQdroop{PSY.VirtualInertia, PSY.ReactivePowerDroop},
+        PSY.OuterControl{PSY.VirtualInertia, PSY.ReactivePowerDroop},
     )
     vpll_d = device_states[external_ix[1]]
     vpll_q = device_states[external_ix[2]]
     ϵ_pll = device_states[external_ix[3]]
-    vod = device_states[external_ix[4]]
-    voq = device_states[external_ix[5]]
-    iod = device_states[external_ix[6]]
-    ioq = device_states[external_ix[7]]
+    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]]
 
     #Obtain inner variables for component
     ω_pll = get_inner_vars(device)[ω_freq_estimator_var]
 
     #Get Active Power Controller parameters
-    outer_control = PSY.get_outercontrol(device)
+    outer_control = PSY.get_outer_control(device)
     active_power_control = PSY.get_active_power(outer_control)
     Ta = PSY.get_Ta(active_power_control) #VSM Inertia constant
     kd = PSY.get_kd(active_power_control) #VSM damping constant
@@ -60,28 +60,34 @@ function mdl_outer_ode!(
     #Obtain indices for component w/r to device
     local_ix = get_local_state_ix(
         device,
-        PSY.VirtualInertiaQdroop{PSY.VirtualInertia, PSY.ReactivePowerDroop},
+        PSY.OuterControl{PSY.VirtualInertia, PSY.ReactivePowerDroop},
     )
 
     #Define internal states for frequency estimator
     internal_states = @view device_states[local_ix]
-    δω_vsm = internal_states[1]
-    δθ_vsm = internal_states[2]
+    ω_oc = internal_states[1]
+    θ_oc = internal_states[2]
     qm = internal_states[3]
 
+    #Obtain additional expressions
+    p_elec_out = Id_filter * Vd_filter + Iq_filter * Vq_filter
+
     #Compute 3 states ODEs
-    output_ode[local_ix[1]] = (
-        -iod * vod / Ta - ioq * voq / Ta +
-        kd * kp_pll * atan(vpll_q / vpll_d) / Ta +
-        kd * ki_pll * ϵ_pll / Ta - (kd + kω) * δω_vsm / Ta +
-        p_ref / Ta +
-        kω * ω_ref / Ta - kω * ω_sys / Ta
-    )
-    output_ode[local_ix[2]] = ωb * δω_vsm
-    output_ode[local_ix[3]] = (-ωf * ioq * vod + ωf * iod * voq - ωf * qm)
+    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)
 
     #Update inner vars
-    get_inner_vars(device)[δdqRI_var] = δθ_vsm
-    get_inner_vars(device)[ω_control_var] = δω_vsm + 1.0
-    get_inner_vars(device)[v_control_var] = V_ref + kq * (q_ref - qm)
+    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/inverter_models/voltage_source_control_models.jl b/src/models/inverter_models/voltage_source_control_models.jl
deleted file mode 100644
index 16f1a77..0000000
--- a/src/models/inverter_models/voltage_source_control_models.jl
+++ /dev/null
@@ -1,104 +0,0 @@
-function mdl_VScontrol_ode!(
-    device_states,
-    output_ode,
-    device::PSY.DynamicInverter{C, O, PSY.CombinedVIwithVZ, DC, P, F},
-) where {
-    C <: PSY.Converter,
-    O <: PSY.OuterControl,
-    DC <: PSY.DCSource,
-    P <: PSY.FrequencyEstimator,
-    F <: PSY.Filter,
-}
-
-    #Obtain external states inputs for component
-    external_ix = get_input_port_ix(device, PSY.CombinedVIwithVZ)
-    iod = device_states[external_ix[1]] #TODO: Should be var referemce?
-    ioq = device_states[external_ix[2]] #TODO: Should be state reference?
-    icvd = device_states[external_ix[3]]
-    icvq = device_states[external_ix[4]]
-    vod = device_states[external_ix[5]]
-    voq = device_states[external_ix[6]]
-
-    #Obtain inner variables for component
-    #vod = get_inner_vars(device)[Vdo_var] #TODO: Should be state reference?
-    #voq = get_inner_vars(device)[Vqo_var] #TODO: Should be state reference?
-    ω_vsm = get_inner_vars(device)[ω_control_var]
-    v_refr = get_inner_vars(device)[v_control_var]
-    vdc = get_inner_vars(device)[Vdc_var]
-
-    #Get Voltage Controller parameters
-    vscontrol = PSY.get_vscontrol(device)
-    filter = PSY.get_filter(device)
-    kpv = PSY.get_kpv(vscontrol)
-    kiv = PSY.get_kiv(vscontrol)
-    kffi = PSY.get_kffi(vscontrol)
-    cf = PSY.get_cf(filter)
-    rv = PSY.get_rv(vscontrol)
-    lv = PSY.get_lv(vscontrol)
-
-    #Get Current Controller parameters
-    kpc = PSY.get_kpc(vscontrol)
-    kic = PSY.get_kic(vscontrol)
-    kffv = PSY.get_kffv(vscontrol)
-    lf = PSY.get_lf(filter)
-    ωad = PSY.get_ωad(vscontrol)
-    kad = PSY.get_kad(vscontrol)
-
-    #Obtain indices for component w/r to device
-    local_ix = get_local_state_ix(device, PSY.CombinedVIwithVZ)
-
-    #Define internal states for frequency estimator
-    internal_states = @view device_states[local_ix]
-    ξ_d = internal_states[1]
-    ξ_q = internal_states[2]
-    γ_d = internal_states[3]
-    γ_q = internal_states[4]
-    ϕ_d = internal_states[5]
-    ϕ_q = internal_states[6]
-
-    #Inputs (control signals)
-
-    ### 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)
-    vod_ref = (v_refr - rv * iod + ω_vsm * lv * ioq)
-    voq_ref = (-rv * ioq - ω_vsm * lv * iod)
-    #Output Control Signal - Links to SRF Current Controller
-    icvd_ref = (kpv * (vod_ref - vod) + kiv * ξ_d - cf * ω_vsm * voq + kffi * iod)
-
-    icvq_ref = (kpv * (voq_ref - voq) + kiv * ξ_q + cf * ω_vsm * vod + kffi * ioq)
-    #Voltage Control ODEs
-    #PI Integrator (internal state)
-    output_ode[local_ix[1]] = (vod_ref - vod)
-    output_ode[local_ix[2]] = (voq_ref - voq)
-
-    ## SRF Current Control ##
-    #Active Damping
-    #vad_d = kad*(vod-ϕ_d)
-    #vad_q = kad*(voq-ϕ_q)
-    #References for Converter Output Voltage
-    vcvd_ref = (
-        kpc * (icvd_ref - icvd) + kic * γ_d - ω_vsm * lf * icvq + kffv * vod -
-        kad * (vod - ϕ_d)
-    )
-    vcvq_ref = (
-        kpc * (icvq_ref - icvq) + kic * γ_q + ω_vsm * lf * icvd + kffv * voq -
-        kad * (voq - ϕ_q)
-    )
-    #Modulation Commands to Converter
-    #md = vcvd_ref/vdc
-    #mq = vcvq_ref/vdc
-    #Current Control ODEs
-    #PI Integrator (internal state)
-    output_ode[local_ix[3]] = icvd_ref - icvd
-    output_ode[local_ix[4]] = icvq_ref - icvq
-    #Active Damping LPF (internal state)
-    output_ode[local_ix[5]] = ωad * vod - ωad * ϕ_d
-    output_ode[local_ix[6]] = ωad * voq - ωad * ϕ_q
-
-    #Update inner_vars
-    get_inner_vars(device)[md_var] = vcvd_ref / vdc
-    get_inner_vars(device)[mq_var] = vcvq_ref / vdc
-
-end
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..87bb73d 100644
--- a/test/data_tests/dynamic_test_data.jl
+++ b/test/data_tests/dynamic_test_data.jl
@@ -444,12 +444,12 @@ end
 
 ###### Converter Data ######
 
-converter_DAIB() = PSY.AvgCnvFixedDC(
+converter_DAIB() = PSY.AverageConverter(
     690.0, #Rated Voltage
     2.75,
 ) #Rated MVA
 
-converter_case78() = PSY.AvgCnvFixedDC(
+converter_case78() = PSY.AverageConverter(
     138.0, #Rated Voltage
     100.0,
 ) #Rated MVA
@@ -472,7 +472,7 @@ filter_test() = PSY.LCLFilter(
 
 ###### PLL Data ######
 
-pll_test() = PSY.PLL(
+pll_test() = PSY.KauraPLL(
     500.0, #ω_lp: Cut-off frequency for LowPass filter of PLL filter.
     0.084, #k_p: PLL proportional gain
     4.69,
@@ -492,12 +492,11 @@ reactive_droop_test() = PSY.ReactivePowerDroop(
     1000.0,
 ) #ωf:: Reactive power cut-off low pass filter frequency
 
-outer_control_test() =
-    PSY.VirtualInertiaQdroop(virtual_inertia_test(), reactive_droop_test())
+outer_control_test() = PSY.OuterControl(virtual_inertia_test(), reactive_droop_test())
 
 ######## Inner Control ######
 
-vsc_test() = PSY.CombinedVIwithVZ(
+vsc_test() = PSY.CurrentControl(
     0.59, #kpv:: Voltage controller proportional gain
     736.0, #kiv:: Voltage controller integral gain
     0.0, #kffv:: Binary variable enabling the voltage feed-forward in output of current controllers
@@ -550,6 +549,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_case06_DAIB.jl b/test/test_case06_DAIB.jl
index 8b6f2d1..fa4e763 100644
--- a/test/test_case06_DAIB.jl
+++ b/test/test_case06_DAIB.jl
@@ -40,8 +40,8 @@ x0_guess = [
     1.0648, #V2_R
     0.0, #V1_I
     0.001, #V2_I
-    0.0, #δω_vsm
-    0.2, #δθ_vsm
+    1.0, #ω_oc
+    0.2, #θ_oc
     0.025, #qm
     0.0015, #ξ_d
     -0.07, #ξ_q
@@ -52,14 +52,14 @@ x0_guess = [
     1.004, #vpll_d
     0.0, #vpll_q
     0.0, #ε_pll
-    0.1, #δθ_pll
+    0.1, #θ_pll
     0.5, #id_cv
     0.0, #iq_cv
-    0.95, #vod
-    -0.1, #voq
-    0.49, #iod
+    0.95, #Vd_filter
+    -0.1, #Vq_filter
+    0.49, #Id_filter
     -0.1,
-] #ioq
+] #Iq_filter
 
 #Define Fault using Callbacks
 Pref_change = LITS.ControlReferenceChange(1.0, Darco_Inverter, LITS.P_ref_index, 0.7)
@@ -71,6 +71,6 @@ sim = LITS.Simulation(sys, tspan, Pref_change, initial_guess = x0_guess)
 run_simulation!(sim, Sundials.IDA());
 
 #Obtain data for angles
-series = get_state_series(sim, ("DARCO", :δω_vsm))
+series = get_state_series(sim, ("DARCO", :ω_oc))
 
 @test sim.solution.retcode == :Success
diff --git a/test/test_case07_4th_order_Inverter.jl b/test/test_case07_4th_order_Inverter.jl
index f78d8eb..ac8c284 100644
--- a/test/test_case07_4th_order_Inverter.jl
+++ b/test/test_case07_4th_order_Inverter.jl
@@ -57,8 +57,8 @@ x0_guess = [
     0.0, #V1_I
     -0.01, #V2_I
     -0.01, #V3_I
-    0.0, #δω_vsm
-    0.2, #δθ_vsm
+    1.0, #ω_oc
+    0.2, #θ_oc
     0.025, #qm
     0.0015, #ξ_d
     -0.07, #ξ_q
@@ -69,13 +69,13 @@ x0_guess = [
     1.004, #vpll_d
     0.0, #vpll_q
     0.0, #ε_pll
-    0.1, #δθ_pll
+    0.1, #θ_pll
     0.5, #id_cv
     0.0, #iq_cv
-    0.95, #vod
-    -0.1, #voq
-    0.49, #iod
-    -0.1, #ioq
+    0.95, #Vd_filter
+    -0.1, #Vq_filter
+    0.49, #Id_filter
+    -0.1, #Iq_filter
     1.0, #eq_p
     0.47, #ed_p
     0.6, #δ
diff --git a/test/test_case08_staticbranches.jl b/test/test_case08_staticbranches.jl
index 35bf1d6..9198239 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 ########################
 
@@ -61,8 +61,8 @@ x0_guess = [
     0.0, #V1_I
     -0.01, #V2_I
     -0.01, #V3_I
-    0.0, #δω_vsm
-    0.2, #δθ_vsm
+    1.0, #ω_oc
+    0.2, #θ_oc
     0.025, #qm
     0.0015, #ξ_d
     -0.07, #ξ_q
@@ -73,13 +73,13 @@ x0_guess = [
     1.004, #vpll_d
     0.0, #vpll_q
     0.0, #ε_pll
-    0.1, #δθ_pll
+    0.1, #θ_pll
     0.5, #id_cv
     0.0, #iq_cv
-    0.95, #vod
-    -0.1, #voq
-    0.49, #iod
-    -0.1, #ioq
+    0.95, #Vd_filter
+    -0.1, #Vq_filter
+    0.49, #Id_filter
+    -0.1, #Iq_filter
     1.0, #eq_p
     0.47, #ed_p
     0.6, #δ
diff --git a/test/test_case09_dynbranches.jl b/test/test_case09_dynbranches.jl
index d52bfe6..cf5a06a 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 ########################
 
@@ -61,8 +61,8 @@ x0_guess = [
     0.0, #V1_I
     -0.01, #V2_I
     -0.01, #V3_I
-    0.0, #δω_vsm
-    0.2, #δθ_vsm
+    1.0, #ω_oc
+    0.2, #θ_oc
     0.025, #qm
     0.0015, #ξ_d
     -0.07, #ξ_q
@@ -73,13 +73,13 @@ x0_guess = [
     1.004, #vpll_d
     0.0, #vpll_q
     0.0, #ε_pll
-    0.1, #δθ_pll
+    0.1, #θ_pll
     0.5, #id_cv
     0.0, #iq_cv
-    0.95, #vod
-    -0.1, #voq
-    0.49, #iod
-    -0.1, #ioq
+    0.95, #Vd_filter
+    -0.1, #Vq_filter
+    0.49, #Id_filter
+    -0.1, #Iq_filter
     1.0, #eq_p
     0.47, #ed_p
     0.6, #δ