diff --git a/docs/make.jl b/docs/make.jl
index 1c60205a3..78824ffa4 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -14,6 +14,8 @@ makedocs(;
"Dynamic Lines Simulation" => "tutorials/tutorial_dynamic_lines.md",
"Inverter Modeling" => "tutorials/tutorial_inverter_modeling.md",
"240 WECC solver comparison" => "tutorials/tutorial_240bus.md",
+ "Small-Signal Analysis" => "tutorials/tutorial_continuation_pf.md",
+ "Active Load Model" => "tutorials/tutorial_activeload.md",
],
"Models" => "models.md",
"Initialization" => "initialization.md",
diff --git a/docs/src/tutorials/figs/pv_curve_cpl.svg b/docs/src/tutorials/figs/pv_curve_cpl.svg
new file mode 100644
index 000000000..ed990d89c
--- /dev/null
+++ b/docs/src/tutorials/figs/pv_curve_cpl.svg
@@ -0,0 +1,44 @@
+
+
diff --git a/docs/src/tutorials/figs/pv_curve_cpl_genrou.svg b/docs/src/tutorials/figs/pv_curve_cpl_genrou.svg
new file mode 100644
index 000000000..ce357d7b3
--- /dev/null
+++ b/docs/src/tutorials/figs/pv_curve_cpl_genrou.svg
@@ -0,0 +1,214 @@
+
+
diff --git a/docs/src/tutorials/tutorial_activeload.md b/docs/src/tutorials/tutorial_activeload.md
new file mode 100644
index 000000000..e2c13ea50
--- /dev/null
+++ b/docs/src/tutorials/tutorial_activeload.md
@@ -0,0 +1,181 @@
+# Tutorial Active Constant Power Load model
+
+**Originally Contributed by**: Rodrigo Henriquez-Auba
+
+## Introduction
+
+This tutorial will introduce you to the functionality of `PowerSimulationsDynamics` and `PowerSystems` to explore active load components and a small-signal analysis.
+
+This tutorial presents a simulation of a two-bus system with a GFM inverter at bus 1, and a load on bus 2. We will change the model from a constant power load model, to a constant impedance model and then to a 12-state active constant power load model.
+
+## Dependencies
+
+```@repl tutorial_load
+using PowerSimulationsDynamics;
+PSID = PowerSimulationsDynamics
+using PowerSystemCaseBuilder
+using PowerSystems
+const PSY = PowerSystems;
+```
+
+!!! note
+ `PowerSystemCaseBuilder.jl` is a helper library that makes it easier to reproduce examples in the documentation and tutorials. Normally you would pass your local files to create the system data instead of calling the function `build_system`.
+ For more details visit [PowerSystemCaseBuilder Documentation](https://nrel-sienna.github.io/PowerSystems.jl/stable/tutorials/powersystembuilder/)
+
+`PowerSystems` (abbreviated with `PSY`) is used to properly define the data structure and establish an equilibrium point initial condition with a power flow routine using `PowerFlows`.
+
+## Load the system
+
+
+We load the system using `PowerSystemCaseBuilder.jl`. This system has an inverter located at bus 1.
+
+```@repl tutorial_load
+sys = build_system(PSIDSystems, "2 Bus Load Tutorial Droop")
+```
+
+```@repl tutorial_load
+first(get_components(DynamicInverter, sys))
+```
+
+The load is an exponential load modeled as a constant power load since the coefficients are set to zero.
+
+```@repl tutorial_load
+first(get_components(PSY.ExponentialLoad, sys))
+```
+
+## Run a small-signal analysis
+
+We set up the Simulation. Since the droop model does not have a frequency state, we use a constant frequency reference frame for the network.
+
+```@repl tutorial_load
+sim = Simulation(ResidualModel,
+ sys,
+ mktempdir(),
+ (0.0, 1.0),
+ frequency_reference = ConstantFrequency())
+```
+
+The following provides a summary of eigenvalues for this droop system with a constant power load:
+
+```@repl tutorial_load
+sm = small_signal_analysis(sim);
+df = summary_eigenvalues(sm);
+show(df, allrows = true, allcols = true)
+```
+
+In this inverter model, the filter is modeled using differential equations, and as described in the literature, interfacing a RL filter against an algebraic constant power load usually results in unstable behavior as observed with the positive real part eigenvalue.
+
+## Change to a constant impedance load model
+
+Since the load is an exponential load model we can change the exponent coefficients to 2.0 to behave as a constant impedance model:
+
+```@repl tutorial_load
+# Update load coefficients to 2.0
+load = first(get_components(PSY.ExponentialLoad, sys));
+PSY.set_active_power_coefficient!(load, 2.0);
+PSY.set_reactive_power_coefficient!(load, 2.0);
+```
+
+We then re-run the small-signal analysis:
+
+```@repl tutorial_load
+sim = Simulation(ResidualModel,
+ sys,
+ mktempdir(),
+ (0.0, 1.0),
+ frequency_reference = ConstantFrequency())
+```
+
+```@repl tutorial_load
+sm = small_signal_analysis(sim);
+df = summary_eigenvalues(sm);
+show(df, allrows = true, allcols = true)
+```
+
+Observe that now the system is small-signal stable (since there is only one device the angle of the inverter is used as a reference, and hence is zero).
+
+## Adding a dynamic active load model
+
+To consider a dynamic model in the load it is only required to attach a dynamic component to the static load model. When a dynamic load model is attached, the active and reactive power of the static model are used to define reference parameters to ensure that the dynamic load model matches the static load output power.
+
+Note that when a dynamic model is attached to a static model, the static model does not participate in the dynamic system equations, i.e. the only model interfacing to the network equations is the dynamic model and not the static model (the exponential load).
+
+We define a function to create a active load model with the specific parameters:
+
+```@repl tutorial_load
+# Parameters taken from active load model from N. Bottrell Masters
+# Thesis "Small-Signal Analysis of Active Loads and Large-signal Analysis
+# of Faults in Inverter Interfaced Microgrid Applications", 2014.
+
+# The parameters are then per-unitized to be scalable to represent an aggregation
+# of multiple active loads
+
+# Base AC Voltage: Vb = 380 V
+# Base Power (AC and DC): Pb = 10000 VA
+# Base AC Current: Ib = 10000 / 380 = 26.32 A
+# Base AC Impedance: Zb = 380 / 26.32 = 14.44 Ω
+# Base AC Inductance: Lb = Zb / Ωb = 14.44 / 377 = 0.3831 H
+# Base AC Capacitance: Cb = 1 / (Zb * Ωb) = 0.000183697 F
+# Base DC Voltage: Vb_dc = (√8/√3) Vb = 620.54 V
+# Base DC Current: Ib_dc = Pb / V_dc = 10000/620.54 = 16.12 A
+# Base DC Impedance: Zb_dc = Vb_dc / Ib_dc = 38.50 Ω
+# Base DC Capacitance: Cb_dc = 1 / (Zb_dc * Ωb) = 6.8886315e-5 F
+
+Ωb = 2*pi*60;
+Vb = 380;
+Pb = 10000;
+Ib = Pb / Vb;
+Zb = Vb / Ib;
+Lb = Zb / Ωb;
+Cb = 1 / (Zb * Ωb);
+Vb_dc = sqrt(8)/sqrt(3) * Vb;
+Ib_dc = Pb / Vb_dc;
+Zb_dc = Vb_dc / Ib_dc;
+Cb_dc = 1/(Zb_dc * Ωb);
+
+function active_cpl(load)
+ return PSY.ActiveConstantPowerLoad(
+ name = get_name(load),
+ r_load = 70.0 / Zb_dc,
+ c_dc = 2040e-6 / Cb_dc,
+ rf = 0.1 / Zb,
+ lf = 2.3e-3 / Lb,
+ cf = 8.8e-6 / Cb,
+ rg = 0.03 / Zb,
+ lg = 0.93e-3 / Lb,
+ kp_pll = 0.4,
+ ki_pll = 4.69,
+ kpv = 0.5 * (Vb_dc / Ib_dc),
+ kiv = 150.0 * (Vb_dc / Ib_dc),
+ kpc = 15.0 * (Ib / Vb),
+ kic = 30000.0 * (Ib / Vb),
+ base_power = 100.0,
+ )
+end
+```
+
+We then attach the model to the system:
+
+```@repl tutorial_load
+load = first(get_components(PSY.ExponentialLoad, sys));
+dyn_load = active_cpl(load)
+add_component!(sys, dyn_load, load)
+```
+
+Finally, we set up the simulation:
+
+```@repl tutorial_load
+sim = Simulation(ResidualModel,
+ sys,
+ mktempdir(),
+ (0.0, 1.0),
+ frequency_reference = ConstantFrequency())
+```
+
+```@repl tutorial_load
+sm = small_signal_analysis(sim);
+df = summary_eigenvalues(sm);
+show(df, allrows = true, allcols = true)
+```
+
+Observe the new states of the active load model and that the system is small-signal stable.
\ No newline at end of file
diff --git a/docs/src/tutorials/tutorial_continuation_pf.md b/docs/src/tutorials/tutorial_continuation_pf.md
new file mode 100644
index 000000000..8f327e3f1
--- /dev/null
+++ b/docs/src/tutorials/tutorial_continuation_pf.md
@@ -0,0 +1,198 @@
+# Tutorial Small Signal Analysis with Continuation Power Flow
+
+**Originally Contributed by**: Rodrigo Henriquez-Auba
+
+## Introduction
+
+This tutorial will introduce you to the functionality of `PowerSimulationsDynamics` and `PowerFlows`
+for running small signal analysis in a continuation power flow.
+
+This tutorial presents a simulation of a two-bus system with a generator (represented with a GENROU + SEXS + TGOV1 model) at bus 1, and a load on bus 2. We will increase the load demand to observe the P-V curve and run a small-signal analysis to check if the system satisfies small-signal stability at different operating points.
+
+## Dependencies
+
+```@repl tutorial_pf
+using PowerSimulationsDynamics
+PSID = PowerSimulationsDynamics
+using PowerSystemCaseBuilder
+using PowerSystems
+using PowerFlows
+const PSY = PowerSystems
+using Plots
+gr()
+
+# Disable Logging to avoid excessive information
+using Logging
+Logging.disable_logging(Logging.Info);
+Logging.disable_logging(Logging.Warn);
+```
+
+!!! note
+ `PowerSystemCaseBuilder.jl` is a helper library that makes it easier to reproduce examples in the documentation and tutorials. Normally you would pass your local files to create the system data instead of calling the function `build_system`.
+ For more details visit [PowerSystemCaseBuilder Documentation](https://nrel-sienna.github.io/PowerSystems.jl/stable/tutorials/powersystembuilder/)
+
+`PowerSystems` (abbreviated with `PSY`) is used to properly define the data structure and establish an equilibrium point initial condition with a power flow routine using `PowerFlows`.
+
+## Load the system
+
+We load the system using `PowerSystemCaseBuilder.jl`. This system only have a generator without dynamic data on which we can use `PowerFlows` to generate a P-V (or nose) curve.
+
+```@repl tutorial_pf
+sys_static = build_system(PSIDSystems, "2 Bus Load Tutorial")
+```
+
+Note that this system contains an Exponential Load, but the parameters are set up to zero, so it behaves a Constant Power Load:
+
+```@repl tutorial_pf
+first(get_components(PSY.ExponentialLoad, sys_static))
+```
+
+## Create a P-V curve
+
+The next step is to run multiple power flows and store the voltage at the load and the active power. For this example we will set up the power factor to be unitary (i.e. no reactive power at the load).
+
+```@repl tutorial_pf
+# Create a Power Range to change the power load active power
+P_range = 0.01:0.01:4.6;
+# Choose the power factor
+load_pf = 1.0;
+```
+
+Then create vectors to store the results
+
+```@repl tutorial_pf
+# PV Curve Results
+P_load_p = Vector{Float64}();
+V_load_p = Vector{Float64}();
+```
+
+Then, we run multiple power flows in a for loop by changing the active power of the load:
+
+```@repl tutorial_pf
+for p in P_range
+ # Change the active power and reactive power of the load
+ power = p * 1.0
+ load = get_component(PSY.ExponentialLoad, sys_static, "load1021")
+ set_active_power!(load, power)
+ q_power = power * tan(acos(load_pf))
+ set_reactive_power!(load, q_power)
+ # Run Power Flow
+ status = solve_powerflow!(sys_static)
+ if !status
+ # Finish the loop if the power flow fails
+ print("Power Flow failed at p = $(power)")
+ break
+ end
+ # Obtain the bus voltage information
+ bus = get_component(Bus, sys_static, "BUS 2")
+ Vm = get_magnitude(bus)
+ # Store values in the vectors
+ push!(V_load_p, Vm)
+ push!(P_load_p, power)
+end
+```
+
+The plot can be visualized with:
+
+```@repl tutorial_pf
+plot(P_load_p,
+ V_load_p,
+ label = "PV Curve",
+ xlabel = "Load Power [pu]",
+ ylabel = "Load Bus Voltage [pu]",
+ color = :black
+)
+```
+
+
+
+## Run Small-Signal Analysis besides the Continuation Power Flow
+
+To run a small-signal analysis we require a dynamic model of the machine. We can use `PowerSystemCaseBuilder` to the load the same system, but with a dynamic model for the generator, including a GENROU + SEXS exciter + TGOV1 governor.
+
+```@repl tutorial_pf
+sys = build_system(PSIDSystems, "2 Bus Load Tutorial GENROU")
+```
+
+Here are the components of the generator:
+
+```@repl tutorial_pf
+first(get_components(DynamicGenerator, sys))
+```
+
+Besides the results of the P-V curve, we need to store if the system is small-signal stable or not by looking if there is a positive real part eigenvalue.
+
+```@repl tutorial_pf
+# Vectors to store stability using a boolean (true for stable).
+stable_vec = Vector{Bool}();
+status_vec = Vector{Bool}();
+
+# PV Curve Results
+P_load_p = Vector{Float64}();
+V_load_p = Vector{Float64}();
+```
+
+We then run the main for loop by updating the load active power, but in addition we create a PowerSimulationsDynamics simulation on which we can run a small-signal analysis to check stability.
+
+```@repl tutorial_pf
+for p in P_range
+ # Change the active power and reactive power of the load
+ power = p * 1.0
+ load = get_component(PSY.ExponentialLoad, sys_static, "load1021")
+ set_active_power!(load, power)
+ q_power = power * tan(acos(load_pf))
+ set_reactive_power!(load, q_power)
+ # Run Power Flow
+ status = solve_powerflow!(sys_static)
+ if !status
+ # Finish the loop if the power flow fails
+ print("Power Flow failed at p = $(power)")
+ break
+ end
+ # Obtain the bus voltage information
+ bus = get_component(Bus, sys_static, "BUS 2")
+ Vm = get_magnitude(bus)
+ # Store values in the vectors
+ push!(V_load_p, Vm)
+ push!(P_load_p, power)
+
+ # Update Load Power in the GENROU system
+ load = get_component(PSY.ExponentialLoad, sys, "load1021")
+ set_active_power!(load, power)
+ q_power = power * tan(acos(load_pf))
+ set_reactive_power!(load, q_power)
+ # Construct Simulation
+ sim = Simulation(ResidualModel, sys, mktempdir(), (0.0, 1.0))
+ if sim.status == PSID.BUILT
+ # Check small-signal stability
+ sm = small_signal_analysis(sim).stable
+ # Push results of small-signal stability
+ push!(stable_vec, sm)
+ # Push results if the simulation was able to be constructed
+ push!(status_vec, true)
+ else
+ # Push results if the simulation was not able to be constructed
+ push!(status_vec, false)
+ end
+end
+```
+
+The following plot showcases the P-V curve, while also showcasing (in red) the regions on which the system is small-signal stable.
+
+```@repl tutorial_pf
+
+# Find where is stable and unstable
+dict_true_ixs_p = Vector();
+dict_false_ixs_p = Vector();
+dict_true_ixs_p = findall(x->x, stable_vec);
+dict_false_ixs_p = findall(x->!x, stable_vec);
+
+# Create plot
+true_ixs = dict_true_ixs_p;
+plot(P_load_p, V_load_p, color = :blue, label = "PV Curve", xlabel = "Load Power [pu]", ylabel = "Load Bus Voltage [pu]")
+plot!(Plots.scatter!(P_load_p[true_ixs] , V_load_p[true_ixs], markerstrokewidth= 0, label = "GENROU SSA"))
+```
+
+
+
+This results is consistent with most of the literature for dynamic generator models supplying constant power loads, on which by increasing the active power of the load, produce critical eigenvalues which cross the ``j\omega`` axis at some point. This is called a Hopf Bifurcation, in this case a subcritical one since the limit cycles are unstable.
\ No newline at end of file
diff --git a/src/base/simulation_inputs.jl b/src/base/simulation_inputs.jl
index 0291c59c9..925041dfe 100644
--- a/src/base/simulation_inputs.jl
+++ b/src/base/simulation_inputs.jl
@@ -17,6 +17,7 @@ struct SimulationInputs
mass_matrix::LinearAlgebra.Diagonal{Float64}
global_vars_update_pointers::Dict{Int, Int}
global_state_map::MAPPING_DICT
+ global_inner_var_map::Dict{String, Dict}
function SimulationInputs(
sys::PSY.System,
@@ -89,6 +90,7 @@ struct SimulationInputs
mass_matrix,
global_vars,
MAPPING_DICT(),
+ Dict{String, Dict}(),
)
end
end
diff --git a/src/post_processing/post_proc_common.jl b/src/post_processing/post_proc_common.jl
index e1a3acd04..60aefaf49 100644
--- a/src/post_processing/post_proc_common.jl
+++ b/src/post_processing/post_proc_common.jl
@@ -23,6 +23,27 @@ function make_global_state_map(inputs::SimulationInputs)
return dic
end
+function _get_inner_vars_map(wrapper::DynamicWrapper{T}) where {T <: PSY.DynamicGenerator}
+ index = get_inner_vars_index(wrapper)
+ inner_vars_enums = instances(generator_inner_vars)
+ return Dict(inner_vars_enums .=> index)
+end
+
+function _get_inner_vars_map(wrapper::DynamicWrapper{T}) where {T <: PSY.DynamicInverter}
+ index = get_inner_vars_index(wrapper)
+ inner_vars_enums = instances(inverter_inner_vars)
+ return Dict(inner_vars_enums .=> index)
+end
+
+function make_inner_vars_map(inputs::SimulationInputs)
+ map = inputs.global_inner_var_map
+ device_wrappers = get_dynamic_injectors(inputs)
+ for d in device_wrappers
+ map[PSY.get_name(d.device)] = _get_inner_vars_map(d)
+ end
+ return map
+end
+
function get_state_from_ix(global_index::MAPPING_DICT, idx::Int)
for (name, device_ix) in global_index
if idx ∈ values(device_ix)