Typical Simulation Workflow

The goal of this tutorial is to get you started with PowerDynamics.jl. We'll walk you through the different "stages" of a typical simulation workflow while introducing key terminology along the way.

This tutorial can be downloaded as a normal Julia script here.

The system to model is a simple 3 bus system:

• Bus 1: ideal droop inverter
• Bus 2: a constant Y load
• Bus 3: a second constant Y load

All buses are connected with standard pi-model power lines.

    ╭───────╮
2 ┯━┿       ┿━┯ 3
  ↓ │   ╭───╯ ↓
    ┷━┯━┷ 1
      │
     (~)

using PowerDynamics
using PowerDynamics: Library
using OrdinaryDiffEqRosenbrock
using OrdinaryDiffEqNonlinearSolve
using CairoMakie

Stage I: Defining Dynamical Models

The first phase of the workflow is typically to define your dynamical models. Models may either come from a library (like the IdealDroopInverter and ConstantYLoad below) or they can be defined by the user.

Inverter Model

In this case, we want the first bus to be an ideal droop inverter.

Often, model definition will be a multi step process:

First: we define an "injector model", in this case our inverter:

(t) ┌────────────────┐
 o──┤ Droop Inverter │
    └────────────────┘

inverter_model = Library.IdealDroopInverter(; name=:droop, Vset=1)

This model is an equation-based/symbolic model representing the dynamics. It is based on the great ModelingToolkit.jl Library.

Second: we build a "Bus Model", which connects the injector to a Busbar. This model still lives in the equation-based/symbolic domain.

┌───────────────────────────────┐
│BusModel                       │
│┌────────┐   ┌────────────────┐│
││ BusBar ├─o─┤ Droop inverter ││
│└────────┘   └────────────────┘│
└───────────────────────────────┘

bus_model = MTKBus(inverter_model; name=:invbus)

Third: we compile the symbolic model into a julia function for numeric simulation. Doing so, we get a VertexModel, which is an object from our backend NetworkDynamics.jl

           ╔════════════════════════════════╗
 Network   ║ VertexModel (compiled)         ║
interface  ║  ┌───────────────────────────┐ ║
           ║  │BusModel                   │ ║
 current ────→│┌──────┐ ┌────────────────┐│ ║
           ║  ││BusBar├o┤ Droop inverter ││ ║
 voltage ←────│└──────┘ └────────────────┘│ ║
           ║  └───────────────────────────┘ ║
           ╚════════════════════════════════╝

bus1 = compile_bus(MTKBus(inverter_model); vidx=1)

VertexModel :bus NoFeedForward() @ Vertex 1
 ├─ 2 inputs:  [busbar₊i_r, busbar₊i_i]
 ├─ 3 states:  [droop₊δ≈0, droop₊Qfilt≈1, droop₊Pfilt≈1]
 ├─ 2 outputs: [busbar₊u_r=1, busbar₊u_i=0]
 └─ 8 params:  [droop₊Pset≈1, droop₊Qset≈0, droop₊Vset=1, droop₊ω₀=1, droop₊Kp=1, droop₊Kq=0.1, droop₊τ_p=1, droop₊τ_q=1]

Note that this model is no longer symbolic. The equations have been reduced and transformed into a nonlinear descriptor model. For more information on the different model types, see the Modeling Concepts docs. You can check out the NetworkDynamics.jl doc on the underlying mathematical model.

In the printout above, you can see that we consider different types of variables in our models:

• the input is always the current flowing from the attached power lines into the bus,
• the output is always the voltage at the busbar,
• the states are dynamical or algebraic states in the sense of a Differential-Algebraic-Equation (DAE) model,
• and parameters are static values that stay mostly constant during simulation and define the system behavior.

There is a 5th class of states not shown above: observables. Observables are time dependent values, which are not states in the sense of a DAE but can be reconstructed from the states, inputs, outputs and parameters. Thus, they don't need to be "solved" for numerically, but they can be reconstructed in post-processing. A simple example of an "observed" state would be a voltage angle or the active and reactive power at some bus.

Load Models

For the two loads, we use the predefined ConstantYLoad model from the Library and compile them:

load_model = Library.ConstantYLoad(; name=:load)
bus2 = compile_bus(MTKBus(load_model); name=:loadbus, vidx=2)
bus3 = compile_bus(MTKBus(load_model); name=:loadbus, vidx=3)

VertexModel :loadbus PureStateMap() @ Vertex 3
 ├─ 2 inputs:  [busbar₊i_r, busbar₊i_i]
 ├─ 2 states:  [busbar₊u_r=1, busbar₊u_i=0]
 |    with diagonal mass matrix [0, 0]
 ├─ 2 outputs: [busbar₊u_r=1, busbar₊u_i=0]
 └─ 2 params:  [load₊B≈1, load₊G≈1]

Power Line Models

Lastly, we need to define three power lines. The workflow is similar to the bus models:

l = MTKLine(Library.PiLine(; name=:piline))
line12 = compile_line(l; src=1, dst=2, piline₊R=0.01)
line13 = compile_line(l; src=1, dst=3, piline₊R=0.01)
line23 = compile_line(l; src=2, dst=3, piline₊R=0.01)

EdgeModel :line PureFeedForward() @ Edge 2=>3
 ├─ 2/2 inputs:  src=[src₊u_r, src₊u_i] dst=[dst₊u_r, dst₊u_i]
 ├─   0 states:  []  
 ├─ 2/2 outputs: src=[src₊i_r, src₊i_i] dst=[dst₊i_r, dst₊i_i]
 └─   9 params:  [piline₊R=0.01, piline₊X=0.1, piline₊G_src=0, piline₊B_src=0, piline₊G_dst=0, piline₊B_dst=0, piline₊r_src=1, piline₊r_dst=1, piline₊active=1]

Now we're all set for the next stage.

Stage II: Initialization

When simulating power systems (or any large dynamical system for that matter), it is quite typical to start from a steady state/equilibrium point. In general, it is not trivial to find such a point for a large nonlinear system. In power systems specifically, it is common to solve a simpler system first – the so-called "power flow" problem.

In the power flow problem, we neglect all the node dynamics and consider only 4 variables at each bus:

• the active power $P$,
• the reactive power $Q$,
• the voltage magnitude $V$ and
• the voltage angle $\theta$.

In the simplest power flow, each bus can then be classified into one of three types:

• Slack Bus: The voltage magnitude and angle are fixed (typically used for one bus in the system), $P$ and $Q$ is considered free.
• PV Bus $P$ and $V$ are fixed, $Q$ and $\theta$ are free. Often used for generator buses or any buses with active voltage control.
• PQ Bus $P$ and $Q$ are fixed, $V$ and $\theta$ are free. Typically used for load buses.

So each component essentially introduces two algebraic equations and two free variables – the system is then solved for the free variables such that all equations are satisfied.

Attaching Power Flow Models

In PowerDynamics.jl, we can attach the power flow models to the dynamic bus models using the set_pfmodel! function.

set_pfmodel!(bus1, pfSlack(V=1))
set_pfmodel!(bus2, pfPQ(P=-0.4, Q=-0.3))
set_pfmodel!(bus3, pfPQ(P=-0.6, Q=-0.2))

Building the Network

Now we can build the network using the Network constructor from NetworkDynamics.jl This constructor takes a list of VertexModels (the buses) and a list of EdgeModels (the powerlines) and connects them to a network. In general, we also need to define the topology of the undelying graph. In this case however, this is not necessary because we told each component at the compile step where it is placed in the network (see the vidx, src and dst arguments above)

nw = Network([bus1, bus2, bus3], [line12, line13, line23])

Network with 7 states and 39 parameters
 ├─ 3 vertices (2 unique types)
 └─ 3 edges (1 unique type)
Edge-Aggregation using SequentialAggregator(+)

The Network object tells us that we've just defined a system with 7 States and 39 parameters. We have 3 vertices of 2 unique types (the inverter bus and the load bus) and 3 edges of a single unique type (all power lines are the same pi-line type).

The "states" and "parameters" already hint at a very important property of PowerDynamics/NetworkDynamics: in the end, the whole network is just a big DAE system of the form

\[\mathbf{M}\,\dot{\mathbf{x}} = f(\mathbf{x}, \mathbf{p}, t)\]

where $\mathbf{x}$ are the states and $\mathbf{p}$ the parameters. This is very important to keep in mind, because it allows us to integrate seamlessly with the whole SciML ecosystem and, most importantly, DifferentialEquations.jl.

The 7 states are essentially just the states of our models stacked on top of eachother. Look at the representation of our Vertex and EdgeModels above to see their contribution:

• Bus 1: 3 states, 8 parameters
• Bus 2 & 3: 2 states, 2 parameters each
• Lines 1,2 and 3: 0 states, 9 parameters each

In sum, we get the 7 states and 39 parameters.

Initializing the System via Power Flow

With the network constructed, we can finally find our equilibrium point. We do so using initialize_from_pf:

s0 = initialize_from_pf(nw; verbose=true, subverbose=true)

┌ Warning: The `alias_u0` keyword argument is deprecated. Please use a NonlinearAliasSpecifier, e.g. `alias = NonlinearAliasSpecifier(alias_u0 = true)`.
└ @ NonlinearSolveBase ~/.julia/packages/NonlinearSolveBase/2E600/src/solve.jl:57
Initializing vertex 1...
[ Info: Initialization problem is fully constrained. Created NonlinearLeastSquaresProblem for [:droop₊δ, :droop₊Qfilt, :droop₊Pfilt, :droop₊Pset, :droop₊Qset]
[ Info: Initialization successful with residual 6.126462320543899e-17

Initializing vertex 2...
[ Info: Initialization problem is overconstrained (2 vars for 4 equations). Create NonlinearLeastSquaresProblem for [:load₊B, :load₊G].
[ Info: Initialization successful with residual 1.7399681585579255e-13

Initializing vertex 3...
[ Info: Initialization problem is overconstrained (2 vars for 4 equations). Create NonlinearLeastSquaresProblem for [:load₊B, :load₊G].
[ Info: Initialization successful with residual 4.85722573273506e-17

Initializing edge 1...
[ Info: No free variables! Residual 0.0

Initializing edge 2...
[ Info: No free variables! Residual 0.0

Initializing edge 3...
[ Info: No free variables! Residual 0.0

Initialized network with residual 1.7402763534895576e-13!

This function actually does quite a lot.

1. It extracts the powerflow model for each component constructing the powerflow problem.
2. It solves the powerflow problem.
3. From the solution, we know the voltages and powers at each bus. Consequently, we also know the currents at each bus. This means we can directly "map" the powerflow solution to the inputs and outputs of each dynamic model.
4. For each dynamic model, we "fix" the inputs and outputs and try to find values for the states (and potentially parameters) such that the component model is in equilibrium. This is done using a nonlinear solver.

In the log statements, we see which variables/parameters were considered free during the initialization of each component. This behavior can be fine-tuned in a lot of ways, which are beyond the scope of this tutorial. However, here we see that, for example, the complex parameter $Y = B + j\,G$ of the constant Y-load was initially left free but then initialized from the powerflow solution. This means that $Y$ is now set in a way that it draws the correct amount of power at the given voltage.

Similarly, we see that the inverter bus had the parameters $P_{set}$ and $Q_{set}$ free, which were also initialized from the powerflow solution. This is important, because we need to achieve power balance in the system, and due to the losses in the lines it's not possible to know the exact power injections a priori.

The return value of the initialize_from_pf function is a so-called NWState object, which wraps flat $x$ and $p$ vectors and provides a lot of helper functions to access and modify states (including observables) and parameters.

Lets inspect this object further:

s0

NWState{Vector{Float64}} of Network (3 vertices, 3 edges)
  ├─ VIndex(1, :droop₊δ)     => -6.126462320543899e-17
  ├─ VIndex(1, :droop₊Qfilt) => 0.567646497521268
  ├─ VIndex(1, :droop₊Pfilt) => 1.0067646497521268
  ├─ VIndex(2, :busbar₊u_r)  => 0.9654281128037676
  ├─ VIndex(2, :busbar₊u_i)  => -0.04388282011689339
  ├─ VIndex(3, :busbar₊u_r)  => 0.9677395909465843
  └─ VIndex(3, :busbar₊u_i)  => -0.051117179883106614
 p = NWParameter([1.00676, 0.567646, 1.0, 1.0, 1.0, 0.1, 1.0, 1.0, -0.321207, 0.428276  …  1.0, 0.01, 0.1, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0])
 t = nothing

At the highest level, we see the values of the 7 states in the network and their symbolic indices. Those indices can be used to access values directly:

s0[VIndex(1, :droop₊Pfilt)]

1.0067646497521268

Often, you want to access observables or parameters instead of states. There is a whole filtering and access mechanism you can use for that. For example, in PD.jl, each bus has the states :busbar₊P and :busbar₊Q. We can inspect them on all vertices using:

s0.v(:, [:busbar₊P, :busbar₊Q])

FilteringProxy for NWState()
  Component filter: VIndex(Colon())
  State filter:     [:busbar₊P, :busbar₊Q]
  Types:  states ✓  parameters ✓  inputs ✓  outputs ✓  observables ✓
Matching Indices:
  ╭ VIndex(1, :busbar₊P)   1.0067646  :bus
  ╰ VIndex(1, :busbar₊Q)   0.5676465  
  ╭ VIndex(2, :busbar₊P)  -0.4        :loadbus
  ╰ VIndex(2, :busbar₊Q)  -0.3        
  ╭ VIndex(3, :busbar₊P)  -0.6        :loadbus
  ╰ VIndex(3, :busbar₊Q)  -0.2        

In the output, we clearly see how the load buses draw exactly the amount of power we specified in the power flow models. On the inverter bus, however, we inject slightly more power than the loads demand to compensate for the line losses.

Similarly, we can access all node parameters at initial state using

s0.v.p

FilteringProxy for NWState()
  Component filter: AllVertices() <- filter further by obj[1], obj["name"], ...
  State filter:     none
  Types:  states ✗  parameters ✓  inputs ✗  outputs ✗  observables ✗
Matching Indices:
  ╭ VIndex(1, :droop₊Pset)   1.0067646   :bus
  │ VIndex(1, :droop₊Qset)   0.5676465   
  │ VIndex(1, :droop₊Vset)   1           
  │ VIndex(1, :droop₊ω₀)     1           
  │ VIndex(1, :droop₊Kp)     1           
  │ VIndex(1, :droop₊Kq)     0.1         
  │ VIndex(1, :droop₊τ_p)    1           
  ╰ VIndex(1, :droop₊τ_q)    1           
  ╭ VIndex(2, :load₊B)      -0.321207    :loadbus
  ╰ VIndex(2, :load₊G)       0.42827601  
  ╭ VIndex(3, :load₊B)      -0.21296241  :loadbus
  ╰ VIndex(3, :load₊G)       0.63888722  

Here we see how $P_{set}$ and $Q_{set}$ of the inverter were initialized in a way, that they match the powerflow solution.

There is a lot more functionality in the NWState objects, see the Symbolic Indexing docs of ND.jl and especially the FilteringProxy for more details.

Stage III: Time Domain Simulation

With the initialized state, we can finally simulate the system in time domain.

Perturbing the System

Since we start from an equilibrium point, we expect the system to stay there if we don't perturb it. Therefore, to get interesting results, we need to perturb the system.

The simplest way to perturb the system is to change a parameter. For example, let's increase the admittance at bus 2 by 10% after 0.1 seconds.

For that, we define a so-called "callback function", more specifically a preset time callback, which is triggered at a specific simulation time and modifies the parameters. The simulation then continues. General information on callbacks in Differential Equations can be found in the DiffEq.jl docs. Specific extensions for NetworkDynamics.jl can be found in the NetworkDynamics.jl callback docs.

We define the callback and attach it to bus 2 (our first load) like this:

affect = ComponentAffect([], [:load₊G, :load₊B]) do u, p, ctx
    @info "Increase load admittance Y by 10% at t=$(ctx.t)"
    p[:load₊G] = p[:load₊G] * 1.1
    p[:load₊B] = p[:load₊B] * 1.1
end
cb = PresetTimeComponentCallback(0.1, affect)
set_callback!(bus2, cb)

VertexModel :loadbus PureStateMap() @ Vertex 2
 ├─ 2 inputs:   [busbar₊i_r, busbar₊i_i]
 ├─ 2 states:   [busbar₊u_r=1, busbar₊u_i=0]
 |    with diagonal mass matrix [0, 0]
 ├─ 2 outputs:  [busbar₊u_r=1, busbar₊u_i=0]
 ├─ 2 params:   [load₊B≈1, load₊G≈1]
 └─ 1 callback: (:load₊G, :load₊B) affected at t=0.1
Powerflow model :pqbus with [pq₊P=-0.4, pq₊Q=-0.3]

With the callback defined, we can finally create and solve the ODEProblem:

prob = ODEProblem(nw, uflat(s0), (0.0, 5.0), copy(pflat(s0)), callback=get_callbacks(nw))
sol = solve(prob, Rodas5P());

[ Info: Increase load admittance Y by 10% at t=0.1

Stage IV: Postprocessing and Visualization

Once we have the solution object, we can use it like any other solution from DifferentialEquations.jl. Most importantly, we can use symbolic indices to access states, parameters and observables.

Plotting Results

For example, we can quickly plot the frequency response of the droop using

lines(sol, idxs=VIndex(1, :droop₊ω); axis=(;xlabel="Time [s]", ylabel="Frequency ω [pu]"))

We clearly see how the increased active power leads to a drop in frequency. This drop is then compensated by the droop control (i.e., we stabilize at a lower frequency).

Of course, we can also create more complex plots, such as this one showing the active and reactive power at each bus:

let
    fig = Figure(size=(1000,600))
    ax = Axis(fig[1, 1]; xlabel="Time [s]", ylabel="Active Power Load [pu]")
    for i in 2:3
        lines!(ax, sol, idxs=VIndex(i, :busbar₊P), color=Cycled(i))
    end
    axislegend(ax)
    ax = Axis(fig[1,2]; xlabel="Time [s]", ylabel="Reactive Power Load [pu]")
    for i in 2:3
        lines!(ax, sol, idxs=VIndex(i, :busbar₊Q), color=Cycled(i))
    end
    axislegend(ax)
    ax = Axis(fig[2,1]; xlabel="Time [s]", ylabel="Active Power Injection [pu]")
    lines!(ax, sol, idxs=VIndex(1, :busbar₊P))
    axislegend(ax)
    ax = Axis(fig[2,2]; xlabel="Time [s]", ylabel="Reactive Power Injection [pu]")
    lines!(ax, sol, idxs=VIndex(1, :busbar₊Q))
    axislegend(ax)
    fig
end

Here we see that the active and reactive power demand shoots up in the beginning after we increase Y. The power demand then slowly decreases again. This is probably due to a drop in voltage, which leads to lower power demand on constant Y loads.

Let's plot the voltage to verify this:

let
    fig = Figure()
    ax = Axis(fig[1, 1]; xlabel="Time [s]", ylabel="Voltage [pu]")
    for i in 1:3
        lines!(ax, sol, idxs=VIndex(i, :busbar₊u_mag), color=Cycled(i))
    end
    axislegend(ax; position=:rc)
    fig
end

Programmatic Access to Variables

Instead of plotting, we can also always use the solution interpolation to access values programmatically. For example

sol(1.0, idxs=VIndex(1,:droop₊ω))

0.9796177657887593

gives us the frequency of the droop inverter at time t=1.0s.

Similarly, we can extract time series by passing a vector of time points rather than a single point:

sol([0.1,0.2,0.3], idxs=VIndex(1:3, :busbar₊u_mag))

t: 3-element Vector{Float64}:
 0.1
 0.2
 0.3
u: 3-element Vector{Vector{Float64}}:
 [1.0000000000000007, 0.9664249287416306, 0.9690886863258026]
 [0.9996829425894866, 0.963758247336889, 0.9675903419318874]
 [0.9993994847401563, 0.963484975411273, 0.9673159834227099]

This code gives us the voltage magnitude at all three buses at the time points 0.1s, 0.2s, and 0.3s.

To deeply inspect a single point, we can also construct a NWState object from the solution for a specific time:

s095 = NWState(sol, 0.95)

NWState{Vector{Float64}} of Network (3 vertices, 3 edges)
  ├─ VIndex(1, :droop₊δ)     => -0.009782234715488626
  ├─ VIndex(1, :droop₊Qfilt) => 0.5860197401082406
  ├─ VIndex(1, :droop₊Pfilt) => 1.0264978211890774
  ├─ VIndex(2, :busbar₊u_r)  => 0.9606923355541416
  ├─ VIndex(2, :busbar₊u_i)  => -0.05547322393118105
  ├─ VIndex(3, :busbar₊u_r)  => 0.9641569090056149
  └─ VIndex(3, :busbar₊u_i)  => -0.06153953104316798
 p = NWParameter([1.00676, 0.567646, 1.0, 1.0, 1.0, 0.1, 1.0, 1.0, -0.353328, 0.471104  …  1.0, 0.01, 0.1, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0])
 t = 0.95

gives us the state at t=0.95s. We can use the state object for inspection as we did before. For example, we can inspect the power at the destination end of all lines at this point in time:

s095.e(:, :dst₊P)

FilteringProxy for NWState()
  Component filter: EIndex(Colon())
  State filter:     :dst₊P
  Types:  states ✓  parameters ✓  inputs ✓  outputs ✓  observables ✓
Matching Indices:
  ● EIndex(1, :dst₊P)   0.48841531  :line
  ● EIndex(2, :dst₊P)   0.54420649  :line
  ● EIndex(3, :dst₊P)   0.05212171  :line

This page was generated using Literate.jl.