Inverter Modeling simulation

To follow along, you can download this tutorial as a Julia script (.jl) or Jupyter notebook (.ipynb).

Originally Contributed by: José Daniel Lara

Introduction

This tutorial will introduce the modeling of an inverter with Virtual Inertia in a multi-machine model of the system. We will load the data directly from PSS/e dynamic files.

The tutorial uses a modified 14-bus system on which all the synchronous machines have been substituted by generators with ESAC1A AVR's and no Turbine Governors.

In the first portion of the tutorial we will simulate the system with the original data and cause a line trip between Buses 2 and 4. In the second part of the simulation, we will switch generator 6 with a battery using an inverter and perform the same fault.

Load the packages

using PowerSimulationsDynamics
using PowerSystemCaseBuilder
using PowerSystems
const PSY = PowerSystems
using PowerFlows
using Logging
using Sundials
using Plots

Create the system using PowerSystemCaseBuilder.jl:

sys = build_system(PSIDSystems, "14 Bus Base Case")

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.

Define Simulation Problem with a 20 second simulation period and the branch trip at t = 1.0:

sim = Simulation(
    ResidualModel, #Type of model used
    sys,         #system
    mktempdir(),       #path for the simulation output
    (0.0, 20.0), #time span
    BranchTrip(1.0, Line, "BUS 02-BUS 04-i_1");
    console_level = Logging.Info,
)

Now that the system is initialized, we can verify the system states for potential issues.

show_states_initial_value(sim)

Voltage Variables
====================
BUS 01
====================
Vm 1.06
θ 0.0
====================
BUS 02
====================
Vm 1.04
θ -0.0711
====================
BUS 03
====================
Vm 1.01
θ -0.1787
====================
BUS 04
====================
Vm 1.0129
θ -0.1458
====================
BUS 05
====================
Vm 1.0165
θ -0.1235
====================
BUS 06
====================
Vm 1.06
θ -0.1949
====================
BUS 07
====================
Vm 1.0438
θ -0.1812
====================
BUS 08
====================
Vm 1.08
θ -0.1656
====================
BUS 09
====================
Vm 1.0263
θ -0.2102
====================
BUS 10
====================
Vm 1.0245
θ -0.2125
====================
BUS 11
====================
Vm 1.0384
θ -0.2059
====================
BUS 12
====================
Vm 1.0436
θ -0.2105
====================
BUS 13
====================
Vm 1.0372
θ -0.2119
====================
BUS 14
====================
Vm 1.0126
θ -0.2291
====================
====================
Differential States
generator-3-1
====================
eq_p 1.0649
ed_p 0.1243
ψ_kd 0.9872
ψ_kq 0.2132
δ 0.034
ω 1.0
Vm 1.01
Vr1 0.006
Vr2 2.419
Ve 1.791
Vr3 -0.0726
====================
Differential States
generator-8-1
====================
eq_p 1.2657
ed_p 0.0462
ψ_kd 1.1584
ψ_kq 0.1748
δ 0.019
ω 1.0
Vm 1.08
Vr1 0.0097
Vr2 3.9162
Ve 2.8839
Vr3 -0.1175
====================
Differential States
generator-1-1
====================
eq_p 1.0604
ed_p -0.0111
ψ_kd 1.0563
ψ_kq 0.1134
δ 0.1684
ω 1.0
Vm 1.06
Vr1 0.0049
Vr2 1.951
Ve 1.4049
Vr3 -0.0585
x_g1 0.3144
x_g2 0.3144
x_g3 0.3144
====================
Differential States
generator-2-1
====================
eq_p 1.1038
ed_p 0.1491
ψ_kd 1.003
ψ_kq 0.2748
δ 0.1963
ω 1.0
Vm 1.04
Vr1 0.0071
Vr2 2.8613
Ve 2.1338
Vr3 -0.0858
====================
Differential States
generator-6-1
====================
eq_p 1.167
ed_p 0.0955
ψ_kd 1.08
ψ_kq 0.3084
δ 0.1387
ω 1.0
Vm 1.06
Vr1 0.0082
Vr2 3.2875
Ve 2.4472
Vr3 -0.0986
====================

We execute the simulation with an additional tolerance for the solver set at 1e-8:

execute!(sim, IDA(); abstol = 1e-8)

SIMULATION_FINALIZED::BUILD_STATUS = 6

Using PowerSimulationsDynamics tools for exploring the results, we can plot all the voltage results for the buses:

result = read_results(sim)
p = plot();
for b in get_components(ACBus, sys)
    voltage_series = get_voltage_magnitude_series(result, get_number(b))
    plot!(
        p,
        voltage_series;
        xlabel = "Time",
        ylabel = "Voltage Magnitude [pu]",
        label = "Bus - $(get_name(b))",
    )
end
p

We can also explore the frequency of the different generators

p2 = plot();
for g in get_components(ThermalStandard, sys)
    state_series = get_state_series(result, (get_name(g), :ω))
    plot!(
        p2,
        state_series;
        xlabel = "Time",
        ylabel = "Speed [pu]",
        label = "$(get_name(g)) - ω",
    )
end
p2

It is also possible to explore the small signal stability of this system we created.

res = small_signal_analysis(sim)

The system is small signal stable

The eigenvalues can be explored

res.eigenvalues

58-element Vector{ComplexF64}:
  -1000.0000000000017 + 0.0im
  -1000.0000000000011 + 0.0im
  -1000.0000000000003 + 0.0im
   -999.9999999999995 + 0.0im
   -999.9999999999983 + 0.0im
   -51.83226146210246 + 0.0im
  -51.702828654896265 + 0.0im
   -51.44178338585979 - 0.018280166539365712im
   -51.44178338585979 + 0.018280166539365712im
  -51.408038558577836 + 0.0im
                      ⋮
  -0.8294145253846161 - 0.04425340347913098im
  -0.8294145253846161 + 0.04425340347913098im
  -0.6364256141508065 + 0.0im
                 -0.5 + 0.0im
 -0.46850110521490756 + 0.0im
  -0.2825569450439846 + 0.0im
 -0.22460226197094324 - 7.678563001298193im
 -0.22460226197094324 + 7.678563001298193im
                  0.0 + 0.0im

Modifying the system and adding storage

Reload the system for this example:

sys = build_system(PSIDSystems, "14 Bus Base Case")

# We want to remove the generator 6 and the dynamic component attached to it.
thermal_gen = get_component(ThermalStandard, sys, "generator-6-1")
remove_component!(sys, get_dynamic_injector(thermal_gen))
remove_component!(sys, thermal_gen)
# We also change

# We can now define our storage device and add it to the system
storage = EnergyReservoirStorage(;
    name = "Battery",
    available = true,
    bus = get_component(Bus, sys, "BUS 06"),
    prime_mover_type = PrimeMovers.BA,
    storage_technology_type = StorageTech.LIB,
    storage_capacity = 100.0,
    storage_level_limits = (min = 0.05, max = 1.0),
    initial_storage_capacity_level = 0.5,
    rating = 1.0,
    active_power = 0.6,
    input_active_power_limits = (min = 0.0, max = 1.0),
    output_active_power_limits = (min = 0.0, max = 1.0),
    efficiency = (in = 0.80, out = 0.90),
    reactive_power = 0.16,
    reactive_power_limits = (min = -1.0, max = 1.0),
    base_power = 25.0,
)

add_component!(sys, storage)

A good sanity check it running a power flow on the system to make sure all the components are properly scaled and that the system is properly balanced. We can use PowerSystems to perform this check. We can get the results back and perform a sanity check.

res = solve_power_flow(ACPowerFlow(), sys)
res["bus_results"]

After verifying that the system works, we can define our inverter dynamics and add it to the battery that has already been stored in the system.

inverter = DynamicInverter(;
    name = get_name(storage),
    ω_ref = 1.0, # ω_ref,
    converter = AverageConverter(; rated_voltage = 138.0, rated_current = 100.0),
    outer_control = OuterControl(
        VirtualInertia(; Ta = 2.0, kd = 400.0, kω = 20.0),
        ReactivePowerDroop(; kq = 0.2, ωf = 1000.0),
    ),
    inner_control = VoltageModeControl(;
        kpv = 0.59,     #Voltage controller proportional gain
        kiv = 736.0,    #Voltage controller integral gain
        kffv = 0.0,     #Binary variable enabling the voltage feed-forward in output of current controllers
        rv = 0.0,       #Virtual resistance in pu
        lv = 0.2,       #Virtual inductance in pu
        kpc = 1.27,     #Current controller proportional gain
        kic = 14.3,     #Current controller integral gain
        kffi = 0.0,     #Binary variable enabling the current feed-forward in output of current controllers
        ωad = 50.0,     #Active damping low pass filter cut-off frequency
        kad = 0.2,
    ),
    dc_source = FixedDCSource(; voltage = 600.0),
    freq_estimator = KauraPLL(;
        ω_lp = 500.0, #Cut-off frequency for LowPass filter of PLL filter.
        kp_pll = 0.084,  #PLL proportional gain
        ki_pll = 4.69   #PLL integral gain
    ),
    filter = LCLFilter(; lf = 0.08, rf = 0.003, cf = 0.074, lg = 0.2, rg = 0.01),
)
add_component!(sys, inverter, storage)

These are the current system components:

sys

Define Simulation problem using the same parameters:

sim = Simulation(
    ResidualModel, #Type of model used
    sys,         #system
    mktempdir(),       #path for the simulation output
    (0.0, 20.0), #time span
    BranchTrip(1.0, Line, "BUS 02-BUS 04-i_1");
    console_level = Logging.Info,
)

We can verify the small signal stability of the system before running the simulation:

res = small_signal_analysis(sim)

The system is small signal stable

Exploring the eigenvalues:

res.eigenvalues

66-element Vector{ComplexF64}:
  -2268.0867822853024 - 6832.691928678376im
  -2268.0867822853024 + 6832.691928678376im
   -2094.489756467523 - 6527.127108333034im
   -2094.489756467523 + 6527.127108333034im
  -1615.3517316545517 - 290.9591007369999im
  -1615.3517316545517 + 290.9591007369999im
  -1000.0000000000019 + 0.0im
  -1000.0000000000015 + 0.0im
  -1000.0000000000007 + 0.0im
  -1000.0000000000002 + 0.0im
                      ⋮
  -1.0079673594546996 + 0.0im
  -0.9363180292979066 + 0.0im
  -0.7915474210230488 + 0.0im
  -0.6203962737229913 + 0.0im
                 -0.5 + 0.0im
 -0.25257135557914234 + 0.0im
  -0.2362820394352636 - 7.659994104964647im
  -0.2362820394352636 + 7.659994104964647im
                  0.0 + 0.0im

We execute the simulation

execute!(sim, IDA(); abstol = 1e-6)

SIMULATION_FINALIZED::BUILD_STATUS = 6

Using PowerSimulationsDynamics tools for exploring the results, we can plot all the voltage results for the buses

result = read_results(sim)
p = plot();
for b in get_components(ACBus, sys)
    voltage_series = get_voltage_magnitude_series(result, get_number(b))
    plot!(
        p,
        voltage_series;
        xlabel = "Time",
        ylabel = "Voltage Magnitude [pu]",
        label = "Bus - $(get_name(b))",
    )
end
p

We can also explore the frequency of the different static generators and storage

p2 = plot();
for g in get_components(ThermalStandard, sys)
    state_series = get_state_series(result, (get_name(g), :ω))
    plot!(
        p2,
        state_series;
        xlabel = "Time",
        ylabel = "Speed [pu]",
        label = "$(get_name(g)) - ω",
    )
end
state_series = get_state_series(result, ("Battery", :ω_oc))
plot!(p2, state_series; xlabel = "Time", ylabel = "Speed [pu]", label = "Battery - ω_vsm");
p2