Modeling Concepts

In general, PowerDynamics models power grids as a set of dynamical systems for both nodes and edges on a graph. Check out the Mathematical Model documentation of NetworkDynamics for the underlying concepts.

The simulation happens entirely in a synchronous dq-frame. Due to their conceptual similarity to complex phasors, variables in this global dq-frame are referenced by subscripts _r and _i (for real and imaginary). This helps distinguish the variables from local dq frames, e.g. a generator model might transform u_r and u_i into u_d and u_q.

Both edge and node models are so-called input-output-systems: the edges receive the voltage of adjacent nodes as an input, the nodes receive the currents on adjacent edges as an input. In general, this leads to the following structure of a bus/node model:

\[\begin{aligned} M_{\mathrm v}\,\frac{\mathrm{d}}{\mathrm{d}t}x_{\mathrm v} &= f^{\mathrm v}\left(x^{\mathrm v}, \sum_k\begin{bmatrix}i^k_r\\ i^k_i\end{bmatrix}, p_{\mathrm v}, t\right)\\ \begin{bmatrix}u_r\\ u_i\end{bmatrix} &= g^{\mathrm v}(x^\mathrm{v},p_{\mathrm v}, t) \end{aligned}\]

where $M_{\mathrm v}$ is the (possibly singular) mass-matrix, $x_{\mathrm v}$ are the internal states and $p_{\mathrm v}$ are parameters. Function $f_{\mathrm v}$ describes the time evolution of the internal states while output equation $g_{\mathrm v}$ defines the output voltage. The input for the system is the sum of all inflowing currents from adjacent lines $k$. Note how vertices are modeled as one-port systems, i.e. they receive the accumulated current from all connected lines, they can't distinguish which line provides which current. For special cases, this limitation might be mitigated using External Inputs.

The edge model on the other hand looks like this:

\[\begin{aligned} M_{\mathrm e}\,\frac{\mathrm{d}}{\mathrm{d}t}x_{\mathrm e} &= f_{\mathrm e}\left(x_{\mathrm e}, \begin{bmatrix} u_r^\mathrm{src}\\u_i^\mathrm{src}\end{bmatrix}, \begin{bmatrix} u_r^\mathrm{dst}\\u_i^\mathrm{dst}\end{bmatrix},p_\mathrm{e}, t\right)\\ \begin{bmatrix}i_r^\mathrm{src}\\i_i^\mathrm{src}\end{bmatrix} &= g^\mathrm{src}_{\mathrm e}\left(x_{\mathrm e}, \begin{bmatrix} u_r^\mathrm{src}\\u_i^\mathrm{src}\end{bmatrix}, \begin{bmatrix} u_r^\mathrm{dst}\\u_i^\mathrm{dst}\end{bmatrix}, p_\mathrm{e}, t\right)\\ \begin{bmatrix}i_r^\mathrm{dst}\\i_i^\mathrm{dst}\end{bmatrix} &= g^\mathrm{dst}_{\mathrm e}\left(x_{\mathrm e}, \begin{bmatrix} u_r^\mathrm{src}\\u_i^\mathrm{src}\end{bmatrix}, \begin{bmatrix} u_r^\mathrm{dst}\\u_i^\mathrm{dst}\end{bmatrix}, p_\mathrm{e}, t\right)\\ \end{aligned}\]

There are a few notable differences: Edges are two port systems, they have two distinct inputs and two distinct outputs. Namely, they receive the dq voltage from both source and destination end and define the current for both ends separately. In very simple systems without losses, those output currents might be just antisymmetric, in general cases however the current on both ends can differ drastically.

The above descriptions are important to understand what's happening inside the package. However, since we use ModelingToolkit to define the individual models a lot of this complexity is hidden from the user. In the following, we'll go through the most important concepts when designing models using ModelingToolkit.

Relationship between ModelingToolkit and NetworkDynamics

A crucial part of using this Library is understanding the relationship between ModelingToolkit models and NetworkDynamics.

In a nutshell, ModelingToolkit models are symbolic models, i.e. they consist of symbolic equations which are not yet "compiled" for use as a numeric model. The modeling in MTK is very flexible and similar to the Modelica language. What we need in the end is models in the structure defined in the equations above. For that, we need the MTK Models to have a specific structure. Then we can use the Bus and Line function to compile the MTK models and create EdgeModel and VertexModel objects from them. Those objects are not symbolic anymore but compiled numeric versions of the symbolically created systems.

ModelingToolkit Models

Terminal

The Terminal─Connector is an important building block for every model. It represents a connection point with constant voltage in dq─coordinates u_r and u_i and enforces the Kirchhoff constraints sum(i_r)=0 and sum(i_i)=0.

Modeling of Buses

Model class Injector

An injector is a class of components with a single Terminal() (called :terminal). Examples for injectors might be Generators, Shunts, Loads.

      ┌───────────┐
(t)   │           │
 o←───┤  Injector │
      │           │
      └───────────┘

The current for injectors is always in injector convention, i.e. positive currents flow out of the injector towards the terminal.

using PowerDynamics, PowerDynamics.Library, ModelingToolkit
using ModelingToolkit: D_nounits as Dt, t_nounits as t
@mtkmodel MyPQLoad begin
    @components begin
        terminal = Terminal()
    end
    @parameters begin
        Pset, [description="Active Power demand"]
        Qset, [description="Reactive Power demand"]
    end
    @variables begin
        P(t), [description="Active Power"]
        Q(t), [description="Reactive Power"]
    end
    @equations begin
        P ~ terminal.u_r*terminal.i_r + terminal.u_i*terminal.i_i
        Q ~ terminal.u_i*terminal.i_r - terminal.u_r*terminal.i_i
        # if possible, it's better for the solver to explicitly provide algebraic equations for the current
        terminal.i_r ~ (Pset*terminal.u_r + Qset*terminal.u_i)/(terminal.u_r^2 + terminal.u_i^2)
        terminal.i_i ~ (Pset*terminal.u_i - Qset*terminal.u_r)/(terminal.u_r^2 + terminal.u_i^2)
    end
end

Model class MTKBus

A MTKBus is a class of models, which are used to describe the dynamic behavior of a full bus in a power grid. Each MTKBus must contain a predefined model of type BusBar() (named :busbar). This busbar represents the connection point to the grid. Optionally, it may contain various injectors. If there are no injectors, the model just describes a junction bus, i.e. a Bus that just satisfies the Kirchhoff constraint for the flows of connected lines.

 ┌───────────────────────────────────┐
 │ MTKBus             ┌───────────┐  │
 │  ┌──────────┐   ┌──┤ Generator │  │
 │  │          │   │  └───────────┘  │
 │  │  BusBar  ├───o                 │
 │  │          │   │  ┌───────────┐  │
 │  └──────────┘   └──┤ Load      │  │
 │                    └───────────┘  │
 └───────────────────────────────────┘

Sometimes it is not possible to connect all injectors directly but instead one needs or wants Branches between the busbar and injector terminal. As long as the :busbar is present at the toplevel, there are few limitations on the overall model complexity.

For simple models (direct connections of a few injectors) it is possible to use the convenience method MTKBus(injectors...) to create the composite model based on provided injector models.

using PowerDynamics, PowerDynamics.Library, ModelingToolkit
@mtkmodel MyMTKBus begin
    @components begin
        busbar = BusBar()
        swing = Swing()
        load = PQLoad()
    end
    @equations begin
        connect(busbar.terminal, swing.terminal)
        connect(busbar.terminal, load.terminal)
    end
end

mybus = MTKBus(Swing(;name=:swing), PQLoad(;name=:load))

Line Modeling

Model class Branch

A branch is the two-port equivalent to an injector. It needs to have two Terminal()s, one is called :src, the other :dst.

Examples for branches are: PI─Model branches, dynamic RL branches or transformers.

      ┌───────────┐
(src) │           │ (dst)
  o←──┤  Branch   ├──→o
      │           │
      └───────────┘

Both ends follow the injector interface, i.e. current leaving the device towards the terminals is always positive.

using PowerDynamics, PowerDynamics.Library, ModelingToolkit
@mtkmodel MyRLine begin
    @components begin
        src = Terminal()
        dst = Terminal()
    end
    @parameters begin
        R=0, [description="Resistance"]
    end
    @equations begin
        dst.i_r ~ (dst.u_r - src.u_r)/R
        dst.i_i ~ (dst.u_i - src.u_i)/R
        src.i_r ~ -dst.i_r
        src.i_i ~ -dst.i_i
    end
end

Model class: MTKLine

Similar to the MTKBus, a MTKLine is a model class which represents a transmission line in the network.

It must contain two LineEnd() instances, one called :src, one called :dst.

 ┌────────────────────────────────────────────────┐
 │ MTKLine          ┌──────────┐                  │
 │  ┌─────────┐  ┌──┤ Branch A ├──┐  ┌─────────┐  │
 │  │ LineEnd │  │  └──────────┘  │  │ LineEnd │  │
 │  │  :src   ├──o                o──┤  :dst   │  │
 │  │         │  │  ┌──────────┐  │  │         │  │
 │  └─────────┘  └──┤ Branch B ├──┘  └─────────┘  │
 │                  └──────────┘                  │
 └────────────────────────────────────────────────┘

Simple line models, which consist only of valid Branch models can be instantiated using the MTKLine(branches...) constructor.

More complex models can be created manually. For example if you want to chain multiple branches between the LineEnds, for example something like

LineEnd(:src) ──o── Transformer ──o── Pi─Line ──o── LineEnd(:dst)

using PowerDynamics, PowerDynamics.Library, ModelingToolkit
@mtkmodel MyMTKLine begin
    @components begin
        src = LineEnd()
        dst = LineEnd()
        branch1 = PiLine()
        branch2 = PiLine()
    end
    @equations begin
        connect(src.terminal, branch1.src)
        connect(src.terminal, branch2.src)
        connect(dst.terminal, branch1.dst)
        connect(dst.terminal, branch2.dst)
    end
end

line = MTKLine(PiLine(;name=:branch1), PiLine(;name=:branch2))

From MTK Models to NetworkDynamics

Both MTKLine and MTKBus are still purely symbolic ModelingToolkit models. However they have an important property: they possess the correct input-output-structure and variable names to be compiled into VertexModel and EdgeModel models. To do so, PowerDynamics.jl provides the Line and Bus constructors.

Putting the knowledge from this document together, we can start a short simulation of an example network:

using PowerDynamics, PowerDynamics.Library, ModelingToolkit
using Graphs, NetworkDynamics
using OrdinaryDiffEqRosenbrock, OrdinaryDiffEqNonlinearSolve
using CairoMakie

First, we define an MTKBus consisting of two predefined injector models from the Library: a Swing generator model and a PQLoad. To do so, we use the MTKBus(injectors...) constructor.

@named swing = Swing(; Pm=1, V=1, D=0.1)
@named load = PQLoad(; Pset=-.5, Qset=0)
bus1mtk = MTKBus(swing, load; name=:swingbus)

Model swingbus:
Subsystems (3): see hierarchy(swingbus)
  busbar
  swing
  load
Equations (25):
  19 standard: see equations(swingbus)
  6 connecting: see equations(expand_connections(swingbus))
Unknowns (27): see unknowns(swingbus)
  busbar₊u_r(t) [defaults to 1]: bus d-voltage
  busbar₊u_i(t) [defaults to 0]: bus q-voltage
  busbar₊i_r(t): bus d-current (flowing into bus)
  busbar₊i_i(t): bus d-current (flowing into bus)
  busbar₊P(t): bus active power (flowing into network)
  busbar₊Q(t): bus reactive power (flowing into network)
  busbar₊u_mag(t): bus voltage magnitude
  busbar₊u_arg(t): bus voltage argument
  busbar₊i_mag(t): bus current magnitude
  busbar₊i_arg(t): bus current argument
  busbar₊terminal₊u_r(t): d-voltage
  busbar₊terminal₊u_i(t): q-voltage
  busbar₊terminal₊i_r(t): d-current
  busbar₊terminal₊i_i(t): q-current
  swing₊ω(t): Rotor frequency
  swing₊θ(t): Rotor angle
  swing₊Pel(t): Electrical Power injected into the grid
  swing₊terminal₊u_r(t): d-voltage
  swing₊terminal₊u_i(t): q-voltage
  swing₊terminal₊i_r(t): d-current
  swing₊terminal₊i_i(t): q-current
  load₊P(t): Active Power
  load₊Q(t): Reactive Power
  load₊terminal₊u_r(t): d-voltage
  load₊terminal₊u_i(t): q-voltage
  load₊terminal₊i_r(t): d-current
  load₊terminal₊i_i(t): q-current
Parameters (7): see parameters(swingbus)
  swing₊ω_ref [defaults to 1]: Reference frequency
  swing₊Pm [defaults to 1]: Mechanical Power
  swing₊M [defaults to 0.005]: Inertia
  swing₊D [defaults to 0.1]: Damping
  swing₊V [defaults to 1]: Voltage magnitude
  load₊Pset [defaults to -0.5]: Active Power demand
  load₊Qset [defaults to 0]: Reactive Power demand

This results in an MTK model, which fulfills the MTKBus interface and thus can be compiled into an actual VertexModel for simulation:

vertex1f = Bus(bus1mtk) # extract component function

VertexModel :swingbus NoFeedForward()
 ├─ 2 inputs:  [busbar₊i_r, busbar₊i_i]
 ├─ 2 states:  [swing₊ω≈0, swing₊θ≈0]
 ├─ 2 outputs: [busbar₊u_r=1, busbar₊u_i=0]
 └─ 7 params:  [swing₊ω_ref=1, swing₊Pm=1, swing₊M=0.005, swing₊D=0.1, swing₊V=1, load₊Pset=-0.5, load₊Qset=0]

As a second bus in this example, we use a SlackDifferential from the Library. This model is not an Injector but an MTKBus directly, as it does not make sense to connect anything else to a slack bus.

bus2mtk = SlackDifferential(; name=:slackbus)
vertex2f = Bus(bus2mtk) # extract component function

VertexModel :slackbus PureStateMap()
 ├─ 2 inputs:  [busbar₊i_r, busbar₊i_i]
 ├─ 2 states:  [busbar₊u_i=0, busbar₊u_r=1]
 ├─ 2 outputs: [busbar₊u_r=1, busbar₊u_i=0]
 └─ 2 params:  [u_init_r=1, u_init_i=0]

For the connecting line, we instantiate two PiLine from the library. Each PiLine fulfills the Branch interface. Therefore we can define a MTKLine model by putting both Branches in parallel:

@named branch1 = PiLine()
@named branch2 = PiLine()
linemtk = MTKLine(branch1, branch2; name=:powerline)

Model swingbus:
Subsystems (3): see hierarchy(swingbus)
  busbar
  swing
  load
Equations (25):
  19 standard: see equations(swingbus)
  6 connecting: see equations(expand_connections(swingbus))
Unknowns (27): see unknowns(swingbus)
  busbar₊u_r(t) [defaults to 1]: bus d-voltage
  busbar₊u_i(t) [defaults to 0]: bus q-voltage
  busbar₊i_r(t): bus d-current (flowing into bus)
  busbar₊i_i(t): bus d-current (flowing into bus)
  busbar₊P(t): bus active power (flowing into network)
  busbar₊Q(t): bus reactive power (flowing into network)
  busbar₊u_mag(t): bus voltage magnitude
  busbar₊u_arg(t): bus voltage argument
  busbar₊i_mag(t): bus current magnitude
  busbar₊i_arg(t): bus current argument
  busbar₊terminal₊u_r(t): d-voltage
  busbar₊terminal₊u_i(t): q-voltage
  busbar₊terminal₊i_r(t): d-current
  busbar₊terminal₊i_i(t): q-current
  swing₊ω(t): Rotor frequency
  swing₊θ(t): Rotor angle
  swing₊Pel(t): Electrical Power injected into the grid
  swing₊terminal₊u_r(t): d-voltage
  swing₊terminal₊u_i(t): q-voltage
  swing₊terminal₊i_r(t): d-current
  swing₊terminal₊i_i(t): q-current
  load₊P(t): Active Power
  load₊Q(t): Reactive Power
  load₊terminal₊u_r(t): d-voltage
  load₊terminal₊u_i(t): q-voltage
  load₊terminal₊i_r(t): d-current
  load₊terminal₊i_i(t): q-current
Parameters (7): see parameters(swingbus)
  swing₊ω_ref [defaults to 1]: Reference frequency
  swing₊Pm [defaults to 1]: Mechanical Power
  swing₊M [defaults to 0.005]: Inertia
  swing₊D [defaults to 0.1]: Damping
  swing₊V [defaults to 1]: Voltage magnitude
  load₊Pset [defaults to -0.5]: Active Power demand
  load₊Qset [defaults to 0]: Reactive Power demand

Similar to before, we need to compile the MTKModel by calling Line.

edgef = Line(linemtk) # extract component function

EdgeModel :powerline PureFeedForward()
 ├─ 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]
 └─  18 params:  [branch1₊R=0, branch1₊X=0.1, branch1₊G_src=0, branch1₊B_src=0, branch1₊G_dst=0, branch1₊B_dst=0, branch1₊r_src=1, branch1₊r_dst=1, branch1₊active=1, branch2₊R=0, branch2₊X=0.1, branch2₊G_src=0, branch2₊B_src=0, branch2₊G_dst=0, branch2₊B_dst=0, branch2₊r_src=1, branch2₊r_dst=1, branch2₊active=1]

To simulate the system, we place both components on a graph and connect them. We define both graph topology as well as the models for the individual components.

g = complete_graph(2)
nw = Network(g, [vertex1f, vertex2f], edgef)
u0 = NWState(nw) # extract parameters and state from models
u0.v[1, :swing₊θ] = 0 # set missing initial conditions
u0.v[1, :swing₊ω] = 1

Then we can solve the problem

prob = ODEProblem(nw, uflat(u0), (0,1), pflat(u0))
sol = solve(prob, Rodas5P())

And finally we can plot the solution:

fig = Figure();
ax = Axis(fig[1,1])
lines!(ax, sol; idxs=VIndex(1,:busbar₊P), label="Power injection Bus", color=Cycled(1))
lines!(ax, sol; idxs=VIndex(1,:swing₊Pel), label="Power injection Swing", color=Cycled(2))
lines!(ax, sol; idxs=VIndex(1,:load₊P), label="Power injection load", color=Cycled(3))
axislegend(ax)

ax = Axis(fig[2,1])
lines!(ax, sol; idxs=VIndex(1,:busbar₊u_arg), label="swing bus voltage angle", color=Cycled(1))
lines!(ax, sol; idxs=VIndex(2,:busbar₊u_arg), label="slack bus voltage angle", color=Cycled(2))
axislegend(ax)

Internals

Internally, we use different input/output conventions for bus and line models. The predefined models BusBar() and LineEnd() are defined in the following way:

Model: BusBar()

A busbar is a concrete model used in bus modeling. It represents the physical connection within a bus, the thing where all injectors and lines attach.

           ┌──────────┐
i_lines ──→│          │  (t)
           │  Busbar  ├───o
  u_bus ←──│          │
           └──────────┘

It receives the sum of all line currents as an input and equals that to the currents flowing into the terminal. As an output, it forwards the terminal voltage to the backend.

Model: LineEnd()

A LineEnd model is very similar to the BusBar model. It represents one end of a transmission line.

          ┌───────────┐
 u_bus ──→│           │  (t)
          │  LineEnd  ├───o
i_line ←──│           │
          └───────────┘

It has special input/output connectors which handle the network interconnection. The main difference being the different input/output conventions for the network interface.