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
julia> using PowerSimulationsDynamics;julia> PSID = PowerSimulationsDynamicsPowerSimulationsDynamicsjulia> using PowerSystemCaseBuilderjulia> using PowerSystemsjulia> const PSY = PowerSystems;
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.
julia> sys = build_system(PSIDSystems, "2 Bus Load Tutorial Droop")┌ Info: Building new system 2 Bus Load Tutorial Droop from raw data └ sys_descriptor.raw_data = "/home/runner/.julia/artifacts/e72e0ad2722dc7be918ec2879d80adedb4d65576/PowerSystemsTestData-2.1/psid_tests/data_examples/Omib_Load.raw" [ Info: The PSS(R)E parser currently supports buses, loads, shunts, generators, branches, transformers, and dc lines [ Info: The PSS(R)E parser currently supports buses, loads, shunts, generators, branches, transformers, and dc lines [ Info: Parsing PSS(R)E Bus data into a PowerModels Dict... [ Info: Parsing PSS(R)E Load data into a PowerModels Dict... [ Info: Parsing PSS(R)E Shunt data into a PowerModels Dict... [ Info: Parsing PSS(R)E Generator data into a PowerModels Dict... [ Info: Parsing PSS(R)E Branch data into a PowerModels Dict... [ Info: Parsing PSS(R)E Transformer data into a PowerModels Dict... [ Info: Parsing PSS(R)E Two-Terminal and VSC DC line data into a PowerModels Dict... ┌ Warning: This PSS(R)E parser currently doesn't support Storage data parsing... └ @ PowerSystems ~/.julia/packages/PowerSystems/xh3fM/src/parsers/pm_io/psse.jl:985 ┌ Warning: This PSS(R)E parser currently doesn't support Switch data parsing... └ @ PowerSystems ~/.julia/packages/PowerSystems/xh3fM/src/parsers/pm_io/psse.jl:991 [ Info: angmin and angmax values are 0, widening these values on branch 1 to +/- 60.0 deg. ┌ Info: Constructing System from Power Models │ data["name"] = "omib_load" └ data["source_type"] = "pti" [ Info: Reading bus data [ Info: Reading Load data in PowerModels dict to populate System ... [ Info: Reading LoadZones data in PowerModels dict to populate System ... [ Info: Reading generator data [ Info: Reading branch data [ Info: Reading branch data [ Info: Reading DC Line data [ Info: Reading storage data [ Info: Serialized System to /home/runner/.julia/packages/PowerSystemCaseBuilder/vlTjU/data/serialized_system/NoArgs/2 Bus Load Tutorial Droop.json [ Info: Serialized System metadata to /home/runner/.julia/packages/PowerSystemCaseBuilder/vlTjU/data/serialized_system/NoArgs/2 Bus Load Tutorial Droop_metadata.json System ┌───────────────────┬─────────────┐ │ Property │ Value │ ├───────────────────┼─────────────┤ │ Name │ │ │ Description │ │ │ System Units Base │ SYSTEM_BASE │ │ Base Power │ 100.0 │ │ Base Frequency │ 60.0 │ │ Num Components │ 9 │ └───────────────────┴─────────────┘ Static Components ┌─────────────────┬───────┬────────────────────────┬───────────────┐ │ Type │ Count │ Has Static Time Series │ Has Forecasts │ ├─────────────────┼───────┼────────────────────────┼───────────────┤ │ ACBus │ 2 │ false │ false │ │ Arc │ 1 │ false │ false │ │ Area │ 1 │ false │ false │ │ ExponentialLoad │ 1 │ false │ false │ │ Line │ 1 │ false │ false │ │ LoadZone │ 1 │ false │ false │ │ ThermalStandard │ 1 │ false │ false │ └─────────────────┴───────┴────────────────────────┴───────────────┘ Dynamic Components ┌─────────────────┬───────┐ │ Type │ Count │ ├─────────────────┼───────┤ │ DynamicInverter │ 1 │ └─────────────────┴───────┘
julia> first(get_components(DynamicInverter, sys))generator-101-1 (DynamicInverter{AverageConverter, OuterControl{ActivePowerDroop, ReactivePowerDroop}, VoltageModeControl, FixedDCSource, FixedFrequency, LCLFilter, Nothing}): 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
The load is an exponential load modeled as a constant power load since the coefficients are set to zero.
julia> first(get_components(PSY.ExponentialLoad, sys))load1021 (PowerSystems.ExponentialLoad): name: load1021 available: true bus: BUS 2 (ACBus) active_power: 0.1 reactive_power: 0.032799999999999996 active_power_coefficient: 0.0 reactive_power_coefficient: 0.0 base_power: 100.0 max_active_power: 0.1 max_reactive_power: 0.032799999999999996 services: 0-element Vector{Service} dynamic_injector: nothing ext: Dict{String, Any}() time_series_container: InfrastructureSystems.SystemUnitsSettings: base_value: 100.0 unit_system: UnitSystem.SYSTEM_BASE = 0
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.
julia> sim = Simulation(ResidualModel, sys, mktempdir(), (0.0, 1.0), frequency_reference = ConstantFrequency())[ Info: Unit System changed to UnitSystem.DEVICE_BASE = 1 Simulation Summary ┌─────────────────────────┬────────────────┐ │ Property │ Value │ ├─────────────────────────┼────────────────┤ │ 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:
julia> sm = small_signal_analysis(sim);julia> df = summary_eigenvalues(sm);julia> 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.431047 -2971.58 -6377.16 42.2368 1119.74 3 │ generator-101-1 ii_cnv 0.431047 -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.551981 -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:
julia> # Update load coefficients to 2.0 load = first(get_components(PSY.ExponentialLoad, sys));julia> PSY.set_active_power_coefficient!(load, 2.0);julia> PSY.set_reactive_power_coefficient!(load, 2.0);
We then re-run the small-signal analysis:
julia> sim = Simulation(ResidualModel, sys, mktempdir(), (0.0, 1.0), frequency_reference = ConstantFrequency())[ Info: Unit System changed to UnitSystem.DEVICE_BASE = 1 Simulation Summary ┌─────────────────────────┬────────────────┐ │ Property │ Value │ ├─────────────────────────┼────────────────┤ │ 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 │ └─────────────────────────┴────────────────┘
julia> sm = small_signal_analysis(sim);julia> df = summary_eigenvalues(sm);julia> 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.08 96.513 2867.34 2 │ generator-101-1 ir_filter 0.47969 -17387.8 4716.08 96.513 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:
julia> # 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;julia> Vb = 380;julia> Pb = 10000;julia> Ib = Pb / Vb;julia> Zb = Vb / Ib;julia> Lb = Zb / Ωb;julia> Cb = 1 / (Zb * Ωb);julia> Vb_dc = sqrt(8)/sqrt(3) * Vb;julia> Ib_dc = Pb / Vb_dc;julia> Zb_dc = Vb_dc / Ib_dc;julia> Cb_dc = 1/(Zb_dc * Ωb);julia> 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, ) endactive_cpl (generic function with 1 method)
We then attach the model to the system:
julia> load = first(get_components(PSY.ExponentialLoad, sys));julia> dyn_load = active_cpl(load)load1021 (ActiveConstantPowerLoad): name: load1021 r_load: 1.8178670360110796 c_dc: 29.61400952076375 rf: 0.006925207756232687 lf: 0.0600470617999157 cf: 0.04790501540123547 rg: 0.0020775623268698058 lg: 0.024279898901705045 kp_pll: 0.4 ki_pll: 4.69 kpv: 19.25333333333334 kiv: 5776.000000000003 kpc: 1.038781163434903 kic: 2077.562326869806 base_power: 100.0 ext: Dict{String, Any}() P_ref: 1.0 Q_ref: 1.0 V_ref: 1.0 ω_ref: 1.0 is_filter_differential: 1 states: [:θ_pll, :ϵ_pll, :η, :v_dc, :γd, :γq, :ir_cnv, :ii_cnv, :vr_filter, :vi_filter, :ir_filter, :ii_filter] n_states: 12 internal: InfrastructureSystems.InfrastructureSystemsInternaljulia> add_component!(sys, dyn_load, load)
Finally, we set up the simulation:
julia> sim = Simulation(ResidualModel, sys, mktempdir(), (0.0, 1.0), frequency_reference = ConstantFrequency())[ Info: Unit System changed to UnitSystem.DEVICE_BASE = 1 Simulation Summary ┌─────────────────────────┬────────────────┐ │ Property │ Value │ ├─────────────────────────┼────────────────┤ │ 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 │ └─────────────────────────┴────────────────┘
julia> sm = small_signal_analysis(sim);julia> df = summary_eigenvalues(sm);julia> 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.99 6 │ load1021 vi_filter 0.189073 -2497.54 7644.75 31.0548 1279.99 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.6514 99.625 135.487 12 │ generator-101-1 ξq_ic 0.370911 -848.095 73.6514 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.7406e6 0.00330534 2.77026e5 21 │ load1021 ir_filter 0.410934 -57.5328 1.7406e6 0.00330534 2.77026e5 22 │ generator-101-1 ϕd_ic 0.527692 -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.