Tutorial Active Constant Power Load model

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

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

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

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.

sys = build_system(PSIDSystems, "2 Bus Load Tutorial Droop")

Property	Value
System
Name
Description
System Units Base	SYSTEM_BASE
Base Power	100.0
Base Frequency	60.0
Num Components	9

Type	Count
Static Components
ACBus	2
Arc	1
Area	1
ExponentialLoad	1
Line	1
LoadZone	1
ThermalStandard	1

Type	Count
Dynamic Components
DynamicInverter{AverageConverter, OuterControl, VoltageModeControl, FixedDCSource, FixedFrequency, LCLFilter, Nothing}	1

first(get_components(DynamicInverter, sys))

DynamicInverter: generator-101-1:
   name: generator-101-1
   ω_ref: 1.0
   converter: AverageConverter
   outer_control: OuterControl{ActivePowerDroop, ReactivePowerDroop}
   inner_control: VoltageModeControl
   dc_source: FixedDCSource
   freq_estimator: FixedFrequency
   filter: LCLFilter
   limiter: nothing
   base_power: 100.0
   n_states: 15
   states: [:θ_oc, :p_oc, :q_oc, :ξd_ic, :ξq_ic, :γd_ic, :γq_ic, :ϕd_ic, :ϕq_ic, :ir_cnv, :ii_cnv, :vr_filter, :vi_filter, :ir_filter, :ii_filter]
   ext: Dict{String, Any}()
   InfrastructureSystems.SystemUnitsSettings:
      base_value: 100.0
      unit_system: UnitSystem.SYSTEM_BASE = 0
   has_supplemental_attributes: false
   has_time_series: false

The load is an exponential load modeled as a constant power load since the coefficients are set to zero.

first(get_components(PSY.ExponentialLoad, sys))

ExponentialLoad: load1021:
   name: load1021
   available: true
   bus: ACBus: BUS 2
   active_power: 0.1
   reactive_power: 0.032799999999999996
   α: 0.0
   β: 0.0
   base_power: 100.0
   max_active_power: 0.1
   max_reactive_power: 0.032799999999999996
   conformity: LoadConformity.UNDEFINED = 2
   services: 0-element Vector{Service}
   dynamic_injector: nothing
   ext: Dict{String, Any}()
   InfrastructureSystems.SystemUnitsSettings:
      base_value: 100.0
      unit_system: UnitSystem.SYSTEM_BASE = 0
   has_supplemental_attributes: false
   has_time_series: false

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.

sim = Simulation(ResidualModel,
    sys,
    mktempdir(),
    (0.0, 1.0);
    frequency_reference = ConstantFrequency())

Property	Value
Simulation Summary
Status	BUILT
Simulation Type	Residual Model
Initialized?	Yes
Multimachine system?	No
Time Span	(0.0, 1.0)
Number of States	19
Number of Perturbations	0

The following provides a summary of eigenvalues for this droop system with a constant power load:

sm = small_signal_analysis(sim);
df = summary_eigenvalues(sm);
show(df; allrows = true, allcols = true)

15×6 DataFrame
 Row │ Most Associated            Part. Factor  Real Part  Imag. Part  Damping [%]  Freq [Hz] 
     │ Any                        Any           Any        Any         Any          Any       
─────┼────────────────────────────────────────────────────────────────────────────────────────
   1 │ generator-101-1 ii_filter  0.935337      -17394.7   0.0         100.0        2768.46
   2 │ generator-101-1 ii_cnv     0.431046      -2971.58   -6377.16    42.2368      1119.74
   3 │ generator-101-1 ii_cnv     0.431046      -2971.58   6377.16     42.2368      1119.74
   4 │ generator-101-1 ir_cnv     0.437046      -2495.77   -5920.94    38.8419      1022.64
   5 │ generator-101-1 ir_cnv     0.437046      -2495.77   5920.94     38.8419      1022.64
   6 │ generator-101-1 ξd_ic      0.81633       -532.334   0.0         100.0        84.7237
   7 │ generator-101-1 ξq_ic      0.82722       -459.814   0.0         100.0        73.1817
   8 │ generator-101-1 p_oc       0.99997       -125.664   0.0         100.0        20.0
   9 │ generator-101-1 q_oc       0.999564      -125.449   0.0         100.0        19.9659
  10 │ generator-101-1 ϕd_ic      0.551982      -50.8211   0.0         100.0        8.08843
  11 │ generator-101-1 ϕq_ic      0.551633      -50.7805   0.0         100.0        8.08197
  12 │ generator-101-1 γd_ic      0.663765      -11.3952   0.0         100.0        1.81359
  13 │ generator-101-1 γq_ic      0.664033      -11.3893   0.0         100.0        1.81266
  14 │ generator-101-1 θ_oc       1.0           -0.0       0.0         NaN          0.0
  15 │ generator-101-1 ir_filter  0.94952       17272.3    0.0         -100.0       2748.98

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:

# Update load coefficients to 2.0
load = first(get_components(PSY.ExponentialLoad, sys));
PSY.set_α!(load, 2.0);
PSY.set_β!(load, 2.0);

We then re-run the small-signal analysis:

sim = Simulation(ResidualModel,
    sys,
    mktempdir(),
    (0.0, 1.0);
    frequency_reference = ConstantFrequency())

Property	Value
Simulation Summary
Status	BUILT
Simulation Type	Residual Model
Initialized?	Yes
Multimachine system?	No
Time Span	(0.0, 1.0)
Number of States	19
Number of Perturbations	0

sm = small_signal_analysis(sim);
df = summary_eigenvalues(sm);
show(df; allrows = true, allcols = true)

15×6 DataFrame
 Row │ Most Associated            Part. Factor  Real Part  Imag. Part   Damping [%]  Freq [Hz] 
     │ Any                        Any           Any        Any          Any          Any       
─────┼─────────────────────────────────────────────────────────────────────────────────────────
   1 │ generator-101-1 ir_filter  0.47969       -17387.8   -4716.04     96.5131      2867.34
   2 │ generator-101-1 ir_filter  0.47969       -17387.8   4716.04      96.5131      2867.34
   3 │ generator-101-1 ii_cnv     0.246681      -3144.63   -6413.81     44.0225      1136.88
   4 │ generator-101-1 ii_cnv     0.246681      -3144.63   6413.81      44.0225      1136.88
   5 │ generator-101-1 ir_cnv     0.249406      -2819.8    -6135.34     41.7606      1074.66
   6 │ generator-101-1 ir_cnv     0.249406      -2819.8    6135.34      41.7606      1074.66
   7 │ generator-101-1 ξq_ic      0.454519      -464.905   -16.1851     99.9395      74.0368
   8 │ generator-101-1 ξq_ic      0.454519      -464.905   16.1851      99.9395      74.0368
   9 │ generator-101-1 q_oc       0.991037      -126.084   0.0          100.0        20.067
  10 │ generator-101-1 p_oc       0.99185       -125.663   0.0          100.0        20.0
  11 │ generator-101-1 ϕq_ic      0.491934      -50.7962   -0.0190757   100.0        8.08447
  12 │ generator-101-1 ϕq_ic      0.491934      -50.7962   0.0190757    100.0        8.08447
  13 │ generator-101-1 γq_ic      0.488883      -11.391    -0.00268224  100.0        1.81293
  14 │ generator-101-1 γq_ic      0.488883      -11.391    0.00268224   100.0        1.81293
  15 │ generator-101-1 θ_oc       1.0           -0.0       0.0          NaN          0.0

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:

# 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

active_cpl (generic function with 1 method)

We then attach the model to the system:

load = first(get_components(PSY.ExponentialLoad, sys));
dyn_load = active_cpl(load)
add_component!(sys, dyn_load, load)

Finally, we set up the simulation:

sim = Simulation(ResidualModel,
    sys,
    mktempdir(),
    (0.0, 1.0);
    frequency_reference = ConstantFrequency())

Property	Value
Simulation Summary
Status	BUILT
Simulation Type	Residual Model
Initialized?	Yes
Multimachine system?	No
Time Span	(0.0, 1.0)
Number of States	31
Number of Perturbations	0

sm = small_signal_analysis(sim);
df = summary_eigenvalues(sm);
show(df; allrows = true, allcols = true)

27×6 DataFrame
 Row │ Most Associated            Part. Factor  Real Part  Imag. Part  Damping [%]  Freq [Hz] 
     │ Any                        Any           Any        Any         Any          Any       
─────┼────────────────────────────────────────────────────────────────────────────────────────
   1 │ generator-101-1 vi_filter  0.166385      -2778.16   -6703.69    38.2849      1154.92
   2 │ generator-101-1 vi_filter  0.166385      -2778.16   6703.69     38.2849      1154.92
   3 │ generator-101-1 vr_filter  0.209778      -2607.5    -6520.11    37.1323      1117.61
   4 │ generator-101-1 vr_filter  0.209778      -2607.5    6520.11     37.1323      1117.61
   5 │ load1021 vi_filter         0.189073      -2497.54   -7644.75    31.0548      1279.98
   6 │ load1021 vi_filter         0.189073      -2497.54   7644.75     31.0548      1279.98
   7 │ load1021 vr_filter         0.209271      -2302.82   -8223.34    26.9661      1359.13
   8 │ load1021 vr_filter         0.209271      -2302.82   8223.34     26.9661      1359.13
   9 │ generator-101-1 ii_filter  0.200819      -970.335   -1219.65    62.2586      248.052
  10 │ generator-101-1 ii_filter  0.200819      -970.335   1219.65     62.2586      248.052
  11 │ generator-101-1 ξq_ic      0.370911      -848.095   -73.6513    99.625       135.487
  12 │ generator-101-1 ξq_ic      0.370911      -848.095   73.6513     99.625       135.487
  13 │ load1021 η                 0.379645      -326.767   -222.502    82.6572      62.9183
  14 │ load1021 η                 0.379645      -326.767   222.502     82.6572      62.9183
  15 │ load1021 v_dc              0.161176      -238.623   -832.086    27.5666      137.769
  16 │ load1021 v_dc              0.161176      -238.623   832.086     27.5666      137.769
  17 │ load1021 θ_pll             0.465358      -132.89    0.0         100.0        21.15
  18 │ generator-101-1 p_oc       0.995436      -125.653   0.0         100.0        19.9983
  19 │ generator-101-1 q_oc       0.677379      -122.337   0.0         100.0        19.4706
  20 │ load1021 ir_filter         0.410934      -57.5328   -1.74056e6  0.00330541   2.77019e5
  21 │ load1021 ir_filter         0.410934      -57.5328   1.74056e6   0.00330541   2.77019e5
  22 │ generator-101-1 ϕd_ic      0.527693      -50.8398   0.0         100.0        8.09141
  23 │ generator-101-1 ϕq_ic      0.523655      -50.7749   0.0         100.0        8.08108
  24 │ load1021 ϵ_pll             0.907233      -12.8266   0.0         100.0        2.04141
  25 │ generator-101-1 γd_ic      0.944549      -11.3935   0.0         100.0        1.81333
  26 │ generator-101-1 γq_ic      0.934217      -11.391    0.0         100.0        1.81294
  27 │ generator-101-1 θ_oc       1.0           -0.0       0.0         NaN          0.0

Observe the new states of the active load model and that the system is small-signal stable.