EMT Toy Model Example

This example demonstrates an electromagnetic transient (EMT) simulation of a simple two-bus system using PowerDynamics.jl. The system consists of a slack bus connected to a load bus through an RL transmission line, with the load bus having both a dynamic PQ load and a capacitive shunt element.

We compare our simulation results with PowerFactory reference data to validate the EMT modeling approach.

Pedagogical Example

This is a pedagogical example that demonstrates the modeling concepts in PowerDynamics.jl are generally compatible with EMT simulations. However, this is far from being an actual interesting simulation study. The way we want to handle EMT simulations in PowerDynamics.jl is not yet fully clear and remains an active area of development.

The example serves to illustrate the flexibility of the modeling framework rather than provide a production-ready EMT simulation tool.

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

System Description

The test system includes:

Bus 1: Slack bus (infinite bus with constant voltage)
Bus 2: Load bus with PQ load and shunt capacitor
Transmission line: RL model with distributed capacitance
Disturbance: Load disconnection at t=0.1s

using PowerDynamics
using PowerDynamics.Library
using NetworkDynamics
using ModelingToolkit
using ModelingToolkit: D_nounits as Dt, t_nounits as t
using CSV
using SteadyStateDiffEq
using OrdinaryDiffEqRosenbrock
using DataFrames
using CairoMakie

System Parameters

First, we define the base system parameters and component values.

ω0    = 2π*50    ## Nominal frequency [rad/s]
Sbase = 300      ## Base power [MW]
Vbase = 110      ## Base voltage [kV]

Rline = 1        ## Line resistance [Ω]
Lline = (1/100π) ## Line inductance [H]
Cline = (2e-6)   ## Line capacitance [F]
Pload = -300     ## Load power (negative = consumption) [MW]

# Convert to per-unit values
Rline_pu = Rline / Zbase(Sbase, Vbase)
Lline_pu = Lline / Zbase(Sbase, Vbase)
Cline_pu = Cline / Ybase(Sbase, Vbase)
Pload_pu = Pload / Sbase

Bus Definitions

Slack Bus

The slack bus maintains constant voltage magnitude and angle, representing an infinite bus or strong grid connection.

slackbus = Bus(pfSlack(; V=1), vidx=1)

VertexModel :slackbus NoFeedForward() @ Vertex 1
 ├─ 2 inputs:  [busbar₊i_r, busbar₊i_i]
 ├─ 0 states:  []
 ├─ 2 outputs: [busbar₊u_r=1, busbar₊u_i=0]
 └─ 2 params:  [slack₊V=1, slack₊δ=0]

Dynamic Shunt Capacitor Model

The shunt capacitor is modeled as a dynamic component in the dq-frame. This allows us to observe the three-phase voltages (u_a, u_b, u_c) by transforming from the dq-frame back to abc coordinates.

The capacitor dynamics are given by:

\[\begin{aligned} \frac{du_r}{dt} &= \phantom{-}\omega_0 u_i + \frac{1}{C} i_r \\ \frac{du_i}{dt} &= -\omega_0 u_r + \frac{1}{C} i_i \end{aligned}\]

where the $\omega_0 u_i$ and $-\omega_0 u_r$ terms account for the rotating dq-frame.

@mtkmodel DynamicShunt begin
    @components begin
        terminal = Terminal()
    end
    @variables begin
        u_r(t), [guess=1, description="Real part of voltage"]
        u_i(t), [guess=0, description="Imaginary part of voltage"]
        # Three-phase voltages as observables
        u_a(t), [description="Voltage in a phase"]
        u_b(t), [description="Voltage in b phase"]
        u_c(t), [description="Voltage in c phase"]
    end
    @parameters begin
        C, [description="Capacitance"]
        ω0, [description="Angular frequency of dq Frame"]
    end
    begin
        # Transformation matrix from dq to abc coordinates
        Tdqinv(δ) = [cos(δ)       -sin(δ)
                     cos(δ-2pi/3) -sin(δ-2pi/3)
                     cos(δ+2pi/3) -sin(δ+2pi/3)]
    end
    @equations begin
        # Capacitor dynamics in rotating dq-frame
        Dt(u_r) ~  ω0*u_i + 1/C * terminal.i_r
        Dt(u_i) ~ -ω0*u_r + 1/C * terminal.i_i
        # Terminal connections
        terminal.u_r ~ u_r
        terminal.u_i ~ u_i
        # Transform to three-phase voltages
        [u_a, u_b, u_c] ~ Tdqinv(ω0*t) * [u_r, u_i]
    end
end

Load Bus Components

The load bus combines two components:

A PQ load consuming constant active power (injector model from Library)
A dynamic shunt capacitor representing line charging

@named load = PQLoad(Pset=-Pload_pu, Qset=0)
@named shunt = DynamicShunt(C=Cline_pu, ω0=ω0)
loadbus = Bus(
    MTKBus(load, shunt);
    vidx=2
)

VertexModel :bus PureStateMap() @ Vertex 2
 ├─ 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]
 └─ 4 params:  [load₊Pset=1, load₊Qset=0, shunt₊C=8.0667e-5, shunt₊ω0=314.16]

Transmission Line Model

Dynamic RL Branch

The transmission line is modeled as a dynamic RL branch in the dq-frame. The line current dynamics are given by:

\[\begin{aligned} \frac{di_r}{dt} &= \phantom{-}\omega_0 i_i - \frac{R}{L} i_r + \frac{1}{L}(u_{\text{dst},r} - u_{\text{src},r}) \\ \frac{di_i}{dt} &= -\omega_0 i_r - \frac{R}{L} i_i + \frac{1}{L}(u_{\text{dst},i} - u_{\text{src},i}) \end{aligned}\]

where the voltage difference drives the current through the line impedance.

@mtkmodel DynamicRLBranch begin
    @components begin
        src = Terminal()
        dst = Terminal()
    end
    @variables begin
        i_r(t)=0, [description="Current in real part"]
        i_i(t)=-1, [description="Current in imaginary part"]
    end
    @parameters begin
        R, [description="Resistance"]
        L, [description="Inductance"]
        ω0, [description="Angular frequency of dq Frame"]
    end
    @equations begin
        # RL line dynamics in rotating dq-frame
        Dt(i_r) ~  ω0 * i_i  - R/L * i_r + 1/L*(dst.u_r - src.u_r)
        Dt(i_i) ~ -ω0 * i_r  - R/L * i_i + 1/L*(dst.u_i - src.u_i)
        # Terminal current connections (KCL enforcement)
        src.i_r ~ -i_r  ## Current flows out of source
        src.i_i ~ -i_i
        dst.i_r ~ i_r   ## Current flows into destination
        dst.i_i ~ i_i
    end
end

@named branch = DynamicRLBranch(; R=Rline_pu, L=Lline_pu, ω0=ω0)
line_model = Line(
    MTKLine(branch);
    src=1, dst=2
)

EdgeModel :line PureStateMap() @ Edge 1=>2
 ├─ 2/2 inputs:  src=[src₊u_r, src₊u_i] dst=[dst₊u_r, dst₊u_i]
 ├─   2 states:  [branch₊i_i=-1, branch₊i_r=0]
 ├─ 2/2 outputs: src=[src₊i_r, src₊i_i] dst=[dst₊i_r, dst₊i_i]
 └─   3 params:  [branch₊R=0.024793, branch₊L=7.892e-5, branch₊ω0=314.16]

Network Assembly and Initialization

We assemble the complete network and attempt initialization.

nw = Network([slackbus, loadbus], line_model)
s0 = find_fixpoint(nw; alg=DynamicSS(Rodas5P()))

┌ Warning: At t=5.048047796807022e-5, dt was forced below floating point epsilon 6.776263578034403e-21, and step error estimate = 28375.542759452987. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).
└ @ SciMLBase ~/.julia/packages/SciMLBase/MEH0d/src/integrator_interface.jl:657
┌ Warning: Solver did not finish retcode = Failure (alg = SteadyStateDiffEq.DynamicSS{Rodas5P{0, ADTypes.AutoForwardDiff{nothing, Nothing}, Nothing, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Float64}(Rodas5P{0, ADTypes.AutoForwardDiff{nothing, Nothing}, Nothing, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}(nothing, OrdinaryDiffEqCore.DEFAULT_PRECS, OrdinaryDiffEqCore.trivial_limiter!, OrdinaryDiffEqCore.trivial_limiter!, ADTypes.AutoForwardDiff()), Inf))!
└ @ NetworkDynamics ~/.julia/packages/NetworkDynamics/x440l/src/initialization.jl:54
┌ Error: ErrorException("Could not find fixpoint, solver returned Failure! For debugging, it is advised to manually construct the steady state problem and try different solvers/arguments:\n\nprob = SteadyStateproblem(nw, uflat(nwstat), pflat(nwpara))\nsol = solve(prob, alg; kwargs...)\nx0 = NWState(nw, sol.u)\n\nFor detail see https://docs.sciml.ai/NonlinearSolve/stable/native/steadystatediffeq/\n")
└ @ Main emt_toymodel.md:204

Initialization Challenge

The direct initialization fails due to the stiffness of the PQ load model. When the load current is computed algebraically from $i = P \frac{u}{|u|^2}$, the system becomes numerically challenging to initialize.

To overcome this, we use a "less stiff" load model with dynamics that smooth out the algebraic singularity during initialization.

We create a "less stiff" version of the PQ load that introduces first-order dynamics with a fast time constant (1/1000 s). This smooths the algebraic relation and makes initialization more robust:

\[\begin{aligned} \frac{di_r}{dt} &= 1000 \left( P \frac{u_r}{u_r^2 + u_i^2} - i_r \right) \\ \frac{di_i}{dt} &= 1000 \left( P \frac{u_i}{u_r^2 + u_i^2} - i_i \right) \end{aligned}\]

This approaches the algebraic PQ load behavior but avoids initialization issues.

@mtkmodel LessStiffPQLoad begin
    @components begin
        terminal = Terminal()
    end
    @variables begin
        i_r(t)=0, [description="Current in real part"]
        i_i(t)=0, [description="Current in imaginary part"]
    end
    @parameters begin
        Pset, [description="Active Power demand"]
    end
    @equations begin
        # First-order dynamics with fast time constant
        Dt(i_r) ~ 1e3*(Pset * terminal.u_r/(terminal.u_r^2 + terminal.u_i^2) - i_r)
        Dt(i_i) ~ 1e3*(Pset * terminal.u_i/(terminal.u_r^2 + terminal.u_i^2) - i_i)
        terminal.i_r ~ i_r
        terminal.i_i ~ i_i
    end
end

@named less_stiff_load = LessStiffPQLoad(Pset=-Pload_pu)
less_stiff_loadbus = Bus(
    MTKBus(less_stiff_load, shunt);
    vidx=2
)
less_stiff_nw = Network([slackbus, less_stiff_loadbus], line_model)
less_stiff_s0 = find_fixpoint(less_stiff_nw; alg=DynamicSS(Rodas5P()))

NWState{Vector{Float64}} of Network (2 vertices, 1 edges)
  ├─ VIndex(2, :busbar₊u_i)          => -0.02539048768641113
  ├─ VIndex(2, :busbar₊u_r)          => 0.9745092560101434
  ├─ VIndex(2, :less_stiff_load₊i_i) => -0.02671802718435793
  ├─ VIndex(2, :less_stiff_load₊i_r) => 1.0254613901218117
  ├─ EIndex(1, :branch₊i_i)          => 0.0020218374630564407
  └─ EIndex(1, :branch₊i_r)          => -1.0261048404693773
 p = NWParameter([1.0, 0.0, 1.0, 8.06667e-5, 314.159, 0.0247934, 7.89198e-5, 314.159])
 t = nothing

Initialize Target System

Perfect! The less stiff load initialization worked. Now we use this solution as an initial guess for our target system with the algebraic PQ load.

s0guess = NWState(nw)
# Transfer key state variables from less stiff solution
s0guess[VIndex(2, :busbar₊u_i)] = less_stiff_s0[VIndex(2, :busbar₊u_i)]
s0guess[VIndex(2, :busbar₊u_r)] = less_stiff_s0[VIndex(2, :busbar₊u_r)]
s0guess[EIndex(1, :branch₊i_i)] = less_stiff_s0[EIndex(1, :branch₊i_i)]
s0guess[EIndex(1, :branch₊i_r)] = less_stiff_s0[EIndex(1, :branch₊i_r)]
s0 = find_fixpoint(nw, s0guess; alg=DynamicSS(Rodas5P()))

NWState{Vector{Float64}} of Network (2 vertices, 1 edges)
  ├─ VIndex(2, :busbar₊u_i) => -0.025390487677360737
  ├─ VIndex(2, :busbar₊u_r) => 0.9745092559179812
  ├─ EIndex(1, :branch₊i_i) => 0.002021837487855655
  └─ EIndex(1, :branch₊i_r) => -1.0261048404519089
 p = NWParameter([1.0, 0.0, 1.0, 0.0, 8.06667e-5, 314.159, 0.0247934, 7.89198e-5, 314.159])
 t = nothing

Disturbance Setup

Excellent! The initialization succeeded. Now we set up a disturbance to observe the system's transient response. We'll disable the load at t=0.1s to simulate a sudden load disconnection.

disable_load_affect = ComponentAffect([], [:load₊Pset]) do u, p, ctx
    println("Disabling load at time $(ctx.t)")
    p[:load₊Pset] = 0  ## Set load power to zero
end
set_callback!(loadbus, PresetTimeComponentCallback(0.1, disable_load_affect))

Dynamic Simulation

With the system properly initialized and the disturbance configured, we can now run the electromagnetic transient simulation.

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

Disabling load at time 0.1

Results and Validation

We compare our EMT simulation results with PowerFactory reference data to validate the modeling approach. The comparison focuses on the three-phase voltages at bus 2 during the load disconnection transient.

The thick gray lines show the PowerFactory reference, while our PowerDynamics.jl results are overlaid in color. The close agreement validates our EMT modeling approach.

fig = let
    fig = Figure()
    ax = Axis(fig[1,1];
        title="Three-Phase Voltage at Bus 2",
        xlabel="Time [s]",
        ylabel="Voltage [pu]")
    ts = range(0.09, 0.13; length=2000)

    # Load PowerFactory reference data
    df = CSV.read(
        joinpath(pkgdir(PowerDynamics),"docs","examples", "emt_data_minimal.csv.gz"),
        DataFrame
    )
    # Plot PowerFactory results (thick gray lines)
    lines!(ax, df.t, df.u_2_a; label="PowerFactory A", color=:lightgray, linewidth=5)
    lines!(ax, df.t, df.u_2_b; label="PowerFactory B", color=:lightgray, linewidth=5)
    lines!(ax, df.t, df.u_2_c; label="PowerFactory C", color=:lightgray, linewidth=5)

    # Extract and plot our simulation results
    a = sol(ts, idxs=VIndex(2, :shunt₊u_a)).u
    b = sol(ts, idxs=VIndex(2, :shunt₊u_b)).u
    c = sol(ts, idxs=VIndex(2, :shunt₊u_c)).u
    lines!(ax, ts, a, label="PowerDynamics A", color=Cycled(1))
    lines!(ax, ts, b, label="PowerDynamics B", color=Cycled(2))
    lines!(ax, ts, c, label="PowerDynamics C", color=Cycled(3))

    axislegend(ax, position=:rt)
    xlims!(ax, ts[begin], ts[end])
    fig
end

Detailed View of Transient Response

Let's zoom in on the critical period around the load disconnection to better observe the transient behavior and compare with the reference.

xlims!(0.0995, 0.105)
fig

This page was generated using Literate.jl.