IEEE39 Bus Tutorial - Part II: Initialization
This is the second part of a four-part tutorial series for the IEEE 39-bus test system:
- Part I: Model Creation - Build the network structure with buses, lines, and components
- Part II: Initialization (this tutorial) - Perform power flow calculations and dynamic initialization
- Part III: Dynamic Simulation - Run time-domain simulations and analyze system behavior
- Part IV: Advanced Modeling & Parameter Optimization - Create custom components and optimize system parameters
The goal of this tutorial is to get an understanding of the initialization process in PowerDynamics.jl.
For comprehensive documentation on initialization, see:
If you're looking for the practical initialization approach without diving into implementation details, jump directly to the Initialize all Components section at the end.
This tutorial goes deep into the initialization internals for educational purposes. In practice, you'll typically use the high-level functions shown at the end rather than the detailed step-by-step process demonstrated here.
As a prerequisite, we load part I of the tutorial, which contains the network model:
using PowerDynamics
EXAMPLEDIR = joinpath(pkgdir(PowerDynamics), "docs", "examples")
include(joinpath(EXAMPLEDIR, "ieee39_part1.jl"))
nw # nw object now availableNetwork with 201 states and 1306 parameters
├─ 39 vertices (5 unique types)
└─ 46 edges (1 unique type)
Edge-Aggregation using SequentialAggregator(+)The initialization process in PowerDynamics.jl is a two-step process: first we solve the power flow, then we use the power flow results to initialize the individual network components.
There are shortcut functions to do this (as shown later), but we will go through the steps in detail for educational purposes.
Power Flow
To solve the power flow, we first need to get the power flow model. We can use the function powerflow_model.
pfnw = powerflow_model(nw)Network with 76 states and 814 parameters
├─ 39 vertices (3 unique types)
└─ 46 edges (1 unique type)
Edge-Aggregation using SequentialAggregator(+)The power flow model is a Network object like the original network. It is built from the original network, by calling powerflow_model on the individual components. For example, we have this rather complex dynamic model at bus 30:
nw[VIndex(30)]VertexModel :bus30 PureStateMap() @ Vertex 30
├─ 2 inputs: [busbar₊i_r, busbar₊i_i]
├─ 15 states: [ctrld_gen₊gov₊τ_m₊u≈0, ctrld_gen₊gov₊xg1≈1, ctrld_gen₊gov₊xg2≈0, ctrld_gen₊avr₊vfout≈1, ctrld_gen₊avr₊v_fb≈0, ctrld_gen₊avr₊vr≈0, ctrld_gen₊avr₊vm≈1, ctrld_gen₊machine₊ψ″_q≈0, ctrld_gen₊machine₊ψ″_d≈1, ctrld_gen₊machine₊E′_d≈1, ctrld_gen₊machine₊E′_q≈0, ctrld_gen₊machine₊ω≈1, ctrld_gen₊machine₊δ≈0, busbar₊u_r, busbar₊u_i]
| with diagonal mass matrix [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0]
├─ 2 outputs: [busbar₊u_r, busbar₊u_i]
└─ 42 params: [ctrld_gen₊machine₊R_s=0, ctrld_gen₊machine₊X_d=1, ctrld_gen₊machine₊X_q=0.69, ctrld_gen₊machine₊X′_d=0.31, ctrld_gen₊machine₊X′_q=0.5, ctrld_gen₊machine₊X″_d=0.25, ctrld_gen₊machine₊X″_q=0.25, ctrld_gen₊machine₊X_ls=0.125, ctrld_gen₊machine₊T′_d0=10.2, ctrld_gen₊machine₊T″_d0=0.05, ctrld_gen₊machine₊T′_q0=2, ctrld_gen₊machine₊T″_q0=0.035, ctrld_gen₊machine₊H=4.2, ctrld_gen₊machine₊D=0, ctrld_gen₊machine₊S_b=100, ctrld_gen₊machine₊V_b=16.5, ctrld_gen₊machine₊ω_b=376.99, ctrld_gen₊machine₊Sn=1000, ctrld_gen₊machine₊Vn=16.5, ctrld_gen₊avr₊Tr=0.01, ctrld_gen₊avr₊vref≈1, ctrld_gen₊avr₊Ka=5, ctrld_gen₊avr₊Ke=-0.0485, ctrld_gen₊avr₊Kf=0.04, ctrld_gen₊avr₊Ta=0.06, ctrld_gen₊avr₊Tf=1, ctrld_gen₊avr₊Te=0.25, ctrld_gen₊avr₊vr_min=-1, ctrld_gen₊avr₊vr_max=1, ctrld_gen₊avr₊E1=3.5461, ctrld_gen₊avr₊E2=4.7281, ctrld_gen₊avr₊Se1=0.08, ctrld_gen₊avr₊Se2=0.26, ctrld_gen₊gov₊ω_ref=1, ctrld_gen₊gov₊p_ref≈1, ctrld_gen₊gov₊V_min=0, ctrld_gen₊gov₊V_max=1, ctrld_gen₊gov₊R=0.05, ctrld_gen₊gov₊T1=0.5, ctrld_gen₊gov₊T2=2.1, ctrld_gen₊gov₊T3=7.2, ctrld_gen₊gov₊DT=0]
Powerflow model :pvbus with [pv₊P=2.5, pv₊V=1.0475]From the printout, we can see that a power flow model :pvbus is attached to the generator model.
We can extract the attached PV power flow model by calling powerflow_model
powerflow_model(nw[VIndex(30)])VertexModel :pvbus PureStateMap()
├─ 2 inputs: [busbar₊i_r, busbar₊i_i]
├─ 2 states: [busbar₊u_r=1.0475, busbar₊u_i=0]
| with diagonal mass matrix [0, 0]
├─ 2 outputs: [busbar₊u_r=1.0475, busbar₊u_i=0]
└─ 2 params: [pv₊P=2.5, pv₊V=1.0475]The function powerflow_model checks if there is a power flow model attached (it checks the :pfmodel metadata).
Per component, the powerflow_model function will do the following:
- If the model has the
:pfmodelmetadata set (see Metadata), it will return theVertexModelstored in the metadata. In Part I we set the:pfmodelmetadata using thepfkeyword to the [Bus]. - If the model does not have the
:pfmodelmetadata set, PowerDynamics will check if the model itself is a valid power flow model. If so, it'll just use the dynamic model as the power flow model.
What is a valid power flow model?: A valid power flow model is a model that has no internal dynamics, i.e. it either has dim(model) == 0 OR it has a zero mass matrix (i.e. only constraints).
For example, the PiLine models are completely static:
nw[EIndex(1)]EdgeModel :piline_template PureFeedForward() @ Edge 1=>2
├─ 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]
└─ 16 params: [piline₊R=0.0035, piline₊X=0.0411, piline₊G_src=0, piline₊B_src=0.34935, piline₊G_dst=0, piline₊B_dst=0.34935, piline₊r_src=1, piline₊r_dst=1, piline₊R_fault=0, piline₊X_fault=0, piline₊G_fault=0, piline₊B_fault=0, piline₊pos=0.5, piline₊active=1, piline₊shortcircuit=0, piline₊faultimp=0]This model has no internal states and no :pfmodel metadata.
@assert ispfmodel(nw[EIndex(1)]) # is pf model itself
powerflow_model(nw[EIndex(1)]) === nw[EIndex(1)]trueAs a result, when we call powerflow_model(nw), we get a completely static network, i.e. only constraints, no dynamics.
all(iszero, pfnw.mass_matrix)trueThe fully static network has the form
\[\dot{x} = 0 = f_{\mathrm{nw}}(x, p, t)\]
Where x are the network states (mainly voltages $u_r$ and $u_i$ at the buses) and p are all the parameters such as $P$, $V$ and $Q$ values for bus models and line parameters for branch models.
To solve this root-finding problem, we first need to find an initial guess. The initial guess is prefilled with all the default values from the components. In our case, all parameters and states have default values attached, so the initial state is fully determined.
pfs0 = NWState(pfnw)NWState{Vector{Float64}} of Network (39 vertices, 46 edges)
├─ VIndex(1, :busbar₊u_i) => 0.0
├─ VIndex(1, :busbar₊u_r) => 1.0
├─ VIndex(2, :busbar₊u_i) => 0.0
├─ VIndex(2, :busbar₊u_r) => 1.0
├─ VIndex(3, :busbar₊u_i) => 0.0
├─ VIndex(3, :busbar₊u_r) => 1.0
├─ VIndex(4, :busbar₊u_i) => 0.0
├─ VIndex(4, :busbar₊u_r) => 1.0
⋮
├─ VIndex(36, :busbar₊u_r) => 1.0635000467300415
├─ VIndex(36, :busbar₊u_i) => 0.0
├─ VIndex(37, :busbar₊u_r) => 1.0277999639511108
├─ VIndex(37, :busbar₊u_i) => 0.0
├─ VIndex(38, :busbar₊u_r) => 1.0264999866485596
├─ VIndex(38, :busbar₊u_i) => 0.0
├─ VIndex(39, :busbar₊u_r) => 1.0299999713897705
└─ VIndex(39, :busbar₊u_i) => 0.0
p = NWParameter([0.0, 0.0, 0.0, 0.0, -3.22, -0.024, -5.0, -1.84, 0.0, 0.0 … 0.97561, 1.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0])
t = nothingWith the default state, we can call find_fixpoint, which keeps $p$ constant and tries to find an $x$ such that the root-finding problem stated above is fulfilled:
pfs = find_fixpoint(pfnw, pfs0)NWState{Vector{Float64}} of Network (39 vertices, 46 edges)
├─ VIndex(1, :busbar₊u_i) => -0.1537005168018229
├─ VIndex(1, :busbar₊u_r) => 1.0360171229045203
├─ VIndex(2, :busbar₊u_i) => -0.10513936683885863
├─ VIndex(2, :busbar₊u_r) => 1.0434526576051297
├─ VIndex(3, :busbar₊u_i) => -0.15402146042175108
├─ VIndex(3, :busbar₊u_r) => 1.0185938884459298
├─ VIndex(4, :busbar₊u_i) => -0.16752811347984345
├─ VIndex(4, :busbar₊u_r) => 0.989785207542011
⋮
├─ VIndex(36, :busbar₊u_r) => 1.0522992315769326
├─ VIndex(36, :busbar₊u_i) => 0.15394374497652633
├─ VIndex(37, :busbar₊u_r) => 1.0268824837830115
├─ VIndex(37, :busbar₊u_i) => 0.043418088368074
├─ VIndex(38, :busbar₊u_r) => 1.0169838480457412
├─ VIndex(38, :busbar₊u_i) => 0.13944918574007406
├─ VIndex(39, :busbar₊u_r) => 1.01418619676775
└─ VIndex(39, :busbar₊u_i) => -0.17979515941397015
p = NWParameter([0.0, 0.0, 0.0, 0.0, -3.22, -0.024, -5.0, -1.84, 0.0, 0.0 … 0.97561, 1.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0])
t = nothingAs a result, we get a NWState object which contains the full state for the power flow model.
Since power flow model and dynamic model share the same topology and the same network interface (i.e. nodes create voltages, edges create currents), we can extract the interface values from the power flow model state and apply them to the dynamic model:
interf = interface_values(pfs)OrderedCollections.OrderedDict{NetworkDynamics.SymbolicStateIndex{Int64, Symbol}, Float64} with 524 entries:
VIndex(1, :busbar₊i_r) => 2.22045e-15
VIndex(1, :busbar₊i_i) => -1.58762e-14
VIndex(1, :busbar₊u_r) => 1.03602
VIndex(1, :busbar₊u_i) => -0.153701
VIndex(2, :busbar₊i_r) => -2.22045e-15
VIndex(2, :busbar₊i_i) => -9.54792e-15
VIndex(2, :busbar₊u_r) => 1.04345
VIndex(2, :busbar₊u_i) => -0.105139
VIndex(3, :busbar₊i_r) => 3.08707
VIndex(3, :busbar₊i_i) => -0.490358
VIndex(3, :busbar₊u_r) => 1.01859
VIndex(3, :busbar₊u_i) => -0.154021
VIndex(4, :busbar₊i_r) => 4.60503
VIndex(4, :busbar₊i_i) => -2.63842
VIndex(4, :busbar₊u_r) => 0.989785
VIndex(4, :busbar₊u_i) => -0.167528
VIndex(5, :busbar₊i_r) => 0.0
VIndex(5, :busbar₊i_i) => -2.22045e-14
VIndex(5, :busbar₊u_r) => 0.993976
⋮ => ⋮The interface values give us the inputs and outputs for every component in the network, i.e. for all buses we get values for
busbar₊i_randbusbar₊i_i: current inputbusbar₊u_randbusbar₊u_i: voltage output
For all branches we get:
src₊u_randsrc₊u_i: source side voltage inputdst₊u_randdst₊u_i: destination side voltage inputsrc₊i_randsrc₊i_i: source side current outputdst₊i_randdst₊i_i: destination side current output
With those interface values fixed, we can go over to the second step of the initialization process: the initialization of the dynamic component.
Initialization of Dynamic Components
Initialization of a bus model means, we want to find a steady state of the dynamics given the interface values from the power flow model.
Therefore, the equations for the bus models (recall from Modeling Concepts) become
\[\begin{aligned} M_{\mathrm v}\,\frac{\mathrm{d}}{\mathrm{d}t}x_{\mathrm v} = \color{red}{0} &= f^{\mathrm v}\left(x^{\mathrm v}, \color{red}{\sum_k\begin{bmatrix}i^k_r\\ i^k_i\end{bmatrix}}, p_{\mathrm v}, t\right)\\ \color{red}{\begin{bmatrix}u_r\\ u_i\end{bmatrix}} &= g^{\mathrm v}(x^\mathrm{v},p_{\mathrm v}, t) \end{aligned}\]
where red symbols are fixed by either the power flow solution or our steady state condition (i.e. $\dot{x}=0$). This leaves us with a system of $N=\mathrm{dim}(x) + \mathrm{dim}(u)$ equations – we can solve for $N$ unknowns.
Edge Initialization Equations
\[\begin{aligned} M_{\mathrm e}\,\frac{\mathrm{d}}{\mathrm{d}t}x_{\mathrm e} = \color{red}{0} &= f_{\mathrm e}\left(x_{\mathrm e}, \color{red}{\begin{bmatrix} u_r^\mathrm{src}\\u_i^\mathrm{src}\end{bmatrix}}, \color{red}{\begin{bmatrix} u_r^\mathrm{dst}\\u_i^\mathrm{dst}\end{bmatrix}},p_\mathrm{e}, t\right)\\ \color{red}{\begin{bmatrix}i_r^\mathrm{src}\\i_i^\mathrm{src}\end{bmatrix}} &= g^\mathrm{src}_{\mathrm e}\left(x_{\mathrm e},\color{red}{ \begin{bmatrix} u_r^\mathrm{src}\\u_i^\mathrm{src}\end{bmatrix}}, \color{red}{\begin{bmatrix} u_r^\mathrm{dst}\\u_i^\mathrm{dst}\end{bmatrix}}, p_\mathrm{e}, t\right)\\ \color{red}{\begin{bmatrix}i_r^\mathrm{dst}\\i_i^\mathrm{dst}\end{bmatrix}} &= g^\mathrm{dst}_{\mathrm e}\left(x_{\mathrm e}, \color{red}{\begin{bmatrix} u_r^\mathrm{src}\\u_i^\mathrm{src}\end{bmatrix}}, \color{red}{\begin{bmatrix} u_r^\mathrm{dst}\\u_i^\mathrm{dst}\end{bmatrix}}, p_\mathrm{e}, t\right)\\ \end{aligned}\]
Notably, the unknowns can come from either the set of states or the set of parameters. We divide the set of symbols into two sets:
- fixed symbols have a
defaultmetadata set. They are considered fixed in the solution of the nonlinear system. - free symbols only have a
guessmetadata set. They are considered free in the solution of the nonlinear system.
Let's take a look at the bus model at bus 30:
gen = nw[VIndex(30)]VertexModel :bus30 PureStateMap() @ Vertex 30
├─ 2 inputs: [busbar₊i_r, busbar₊i_i]
├─ 15 states: [ctrld_gen₊gov₊τ_m₊u≈0, ctrld_gen₊gov₊xg1≈1, ctrld_gen₊gov₊xg2≈0, ctrld_gen₊avr₊vfout≈1, ctrld_gen₊avr₊v_fb≈0, ctrld_gen₊avr₊vr≈0, ctrld_gen₊avr₊vm≈1, ctrld_gen₊machine₊ψ″_q≈0, ctrld_gen₊machine₊ψ″_d≈1, ctrld_gen₊machine₊E′_d≈1, ctrld_gen₊machine₊E′_q≈0, ctrld_gen₊machine₊ω≈1, ctrld_gen₊machine₊δ≈0, busbar₊u_r, busbar₊u_i]
| with diagonal mass matrix [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0]
├─ 2 outputs: [busbar₊u_r, busbar₊u_i]
└─ 42 params: [ctrld_gen₊machine₊R_s=0, ctrld_gen₊machine₊X_d=1, ctrld_gen₊machine₊X_q=0.69, ctrld_gen₊machine₊X′_d=0.31, ctrld_gen₊machine₊X′_q=0.5, ctrld_gen₊machine₊X″_d=0.25, ctrld_gen₊machine₊X″_q=0.25, ctrld_gen₊machine₊X_ls=0.125, ctrld_gen₊machine₊T′_d0=10.2, ctrld_gen₊machine₊T″_d0=0.05, ctrld_gen₊machine₊T′_q0=2, ctrld_gen₊machine₊T″_q0=0.035, ctrld_gen₊machine₊H=4.2, ctrld_gen₊machine₊D=0, ctrld_gen₊machine₊S_b=100, ctrld_gen₊machine₊V_b=16.5, ctrld_gen₊machine₊ω_b=376.99, ctrld_gen₊machine₊Sn=1000, ctrld_gen₊machine₊Vn=16.5, ctrld_gen₊avr₊Tr=0.01, ctrld_gen₊avr₊vref≈1, ctrld_gen₊avr₊Ka=5, ctrld_gen₊avr₊Ke=-0.0485, ctrld_gen₊avr₊Kf=0.04, ctrld_gen₊avr₊Ta=0.06, ctrld_gen₊avr₊Tf=1, ctrld_gen₊avr₊Te=0.25, ctrld_gen₊avr₊vr_min=-1, ctrld_gen₊avr₊vr_max=1, ctrld_gen₊avr₊E1=3.5461, ctrld_gen₊avr₊E2=4.7281, ctrld_gen₊avr₊Se1=0.08, ctrld_gen₊avr₊Se2=0.26, ctrld_gen₊gov₊ω_ref=1, ctrld_gen₊gov₊p_ref≈1, ctrld_gen₊gov₊V_min=0, ctrld_gen₊gov₊V_max=1, ctrld_gen₊gov₊R=0.05, ctrld_gen₊gov₊T1=0.5, ctrld_gen₊gov₊T2=2.1, ctrld_gen₊gov₊T3=7.2, ctrld_gen₊gov₊DT=0]
Powerflow model :pvbus with [pv₊P=2.5, pv₊V=1.0475]We have a system with 15 States and 42 parameters. Of those, most are fixed, i.e. have a default metadata set. In the VertexModel printout, defaults are shown with = while guesses are shown with ≈. We can use dump_initial_state to get an overview of the free and set states:
dump_initial_state(gen; obs=false)Inputs:
busbar₊i_i = uninitialized
busbar₊i_r = uninitialized
States:
busbar₊u_i = uninitialized
busbar₊u_r = uninitialized
ctrld_gen₊avr₊v_fb = uninitialized (guess 0)
ctrld_gen₊avr₊vfout = uninitialized (guess 1) (bounds 0.0..Inf)
ctrld_gen₊avr₊vm = uninitialized (guess 1)
ctrld_gen₊avr₊vr = uninitialized (guess 0)
ctrld_gen₊gov₊xg1 = uninitialized (guess 1)
ctrld_gen₊gov₊xg2 = uninitialized (guess 0)
ctrld_gen₊gov₊τ_m₊u = uninitialized (guess 0)
ctrld_gen₊machine₊E′_d = uninitialized (guess 1)
ctrld_gen₊machine₊E′_q = uninitialized (guess 0)
ctrld_gen₊machine₊δ = uninitialized (guess 0)
ctrld_gen₊machine₊ψ″_d = uninitialized (guess 1)
ctrld_gen₊machine₊ψ″_q = uninitialized (guess 0)
ctrld_gen₊machine₊ω = uninitialized (guess 1)
Outputs:
busbar₊u_i = uninitialized
busbar₊u_r = uninitialized
Parameters:
ctrld_gen₊avr₊E1 = 3.5461
ctrld_gen₊avr₊E2 = 4.7281
ctrld_gen₊avr₊Ka = 5
ctrld_gen₊avr₊Ke = -0.0485
ctrld_gen₊avr₊Kf = 0.04
ctrld_gen₊avr₊Se1 = 0.08
ctrld_gen₊avr₊Se2 = 0.26
ctrld_gen₊avr₊Ta = 0.06
ctrld_gen₊avr₊Te = 0.25
ctrld_gen₊avr₊Tf = 1
ctrld_gen₊avr₊Tr = 0.01
ctrld_gen₊avr₊vr_max = 1
ctrld_gen₊avr₊vr_min = -1
ctrld_gen₊avr₊vref = uninitialized (guess 1)
ctrld_gen₊gov₊DT = 0
ctrld_gen₊gov₊R = 0.05
ctrld_gen₊gov₊T1 = 0.5
ctrld_gen₊gov₊T2 = 2.1
ctrld_gen₊gov₊T3 = 7.2
ctrld_gen₊gov₊V_max = 1
ctrld_gen₊gov₊V_min = 0
ctrld_gen₊gov₊p_ref = uninitialized (guess 1)
ctrld_gen₊gov₊ω_ref = 1
ctrld_gen₊machine₊D = 0
ctrld_gen₊machine₊H = 4.2
ctrld_gen₊machine₊R_s = 0
ctrld_gen₊machine₊S_b = 100
ctrld_gen₊machine₊Sn = 1000
ctrld_gen₊machine₊T′_d0 = 10.2
ctrld_gen₊machine₊T′_q0 = 2
ctrld_gen₊machine₊T″_d0 = 0.05
ctrld_gen₊machine₊T″_q0 = 0.035
ctrld_gen₊machine₊V_b = 16.5
ctrld_gen₊machine₊Vn = 16.5
ctrld_gen₊machine₊X_d = 1
ctrld_gen₊machine₊X_ls = 0.125
ctrld_gen₊machine₊X_q = 0.69
ctrld_gen₊machine₊X′_d = 0.31
ctrld_gen₊machine₊X′_q = 0.5
ctrld_gen₊machine₊X″_d = 0.25
ctrld_gen₊machine₊X″_q = 0.25
ctrld_gen₊machine₊ω_b = 376.99Right now, we have 2 free parameters, 2 free inputs, 2 free outputs and 15 free states. We can use initialize_component to find values for the free symbols:
initialize_component(gen)┌ Error: Initialization problem underconstraint. 17 Equations for 21 free variables: [:busbar₊u_r, :busbar₊u_i, :ctrld_gen₊gov₊τ_m₊u, :ctrld_gen₊gov₊xg1, :ctrld_gen₊gov₊xg2, :ctrld_gen₊avr₊vfout, :ctrld_gen₊avr₊v_fb, :ctrld_gen₊avr₊vr, :ctrld_gen₊avr₊vm, :ctrld_gen₊machine₊ψ″_q, :ctrld_gen₊machine₊ψ″_d, :ctrld_gen₊machine₊E′_d, :ctrld_gen₊machine₊E′_q, :ctrld_gen₊machine₊ω, :ctrld_gen₊machine₊δ, :busbar₊u_r, :busbar₊u_i, :busbar₊i_r, :busbar₊i_i, :ctrld_gen₊avr₊vref, :ctrld_gen₊gov₊p_ref]. Consider passing additional constraints using `InitConstraint`.
└ @ Main ieee39_part2.md:188Wait! The initialization failed! Why? Well, we need to apply additional defaults, so called default_overwrites for the inputs/outputs to make the system solvable.
interf_v30 = Dict( # manually define interface values for demonstration
:busbar₊u_r => 1.04573,
:busbar₊u_i => -0.0609188,
:busbar₊i_i => 1.53174,
:busbar₊i_r => -2.30145,
)
initialize_component(gen; default_overrides=interf_v30)[ Info: Apply positivity/negativity conserving variable transformation on [:ctrld_gen₊avr₊vfout] to satisfy bounds.
[ Info: Initialization problem is overconstrained (15 vars for 17 equations). Create NonlinearLeastSquaresProblem for [:ctrld_gen₊gov₊τ_m₊u, :ctrld_gen₊gov₊xg1, :ctrld_gen₊gov₊xg2, :ctrld_gen₊avr₊vfout, :ctrld_gen₊avr₊v_fb, :ctrld_gen₊avr₊vr, :ctrld_gen₊avr₊vm, :ctrld_gen₊machine₊ψ″_q, :ctrld_gen₊machine₊ψ″_d, :ctrld_gen₊machine₊E′_d, :ctrld_gen₊machine₊E′_q, :ctrld_gen₊machine₊ω, :ctrld_gen₊machine₊δ, :ctrld_gen₊avr₊vref, :ctrld_gen₊gov₊p_ref].
[ Info: Initialization successful with residual 2.2358126609832515e-14The non mutating initialize_component returns a dictionary containing a full initialized component state. That can be useful for certain purposes, often it is easier to work with the mutating version of initialization function, which will write the initialized values back to the component metadata (i.e setting the :init property for the symbols).
Let's call the mutating initialize_component! function to write the initialized values back to the component and inspect the initial state using dump_initial_state:
initialize_component!(gen; default_overrides=interf_v30)
dump_initial_state(gen; obs=false)Set additional default for busbar₊u_r: 1.04573
Set additional default for busbar₊u_i: -0.0609188
Set additional default for busbar₊i_i: 1.53174
Set additional default for busbar₊i_r: -2.30145
[ Info: Apply positivity/negativity conserving variable transformation on [:ctrld_gen₊avr₊vfout] to satisfy bounds.
[ Info: Initialization problem is overconstrained (15 vars for 17 equations). Create NonlinearLeastSquaresProblem for [:ctrld_gen₊gov₊τ_m₊u, :ctrld_gen₊gov₊xg1, :ctrld_gen₊gov₊xg2, :ctrld_gen₊avr₊vfout, :ctrld_gen₊avr₊v_fb, :ctrld_gen₊avr₊vr, :ctrld_gen₊avr₊vm, :ctrld_gen₊machine₊ψ″_q, :ctrld_gen₊machine₊ψ″_d, :ctrld_gen₊machine₊E′_d, :ctrld_gen₊machine₊E′_q, :ctrld_gen₊machine₊ω, :ctrld_gen₊machine₊δ, :ctrld_gen₊avr₊vref, :ctrld_gen₊gov₊p_ref].
[ Info: Initialization successful with residual 2.2358126609832515e-14
Inputs:
busbar₊i_i = 1.5317
busbar₊i_r = -2.3014
States:
busbar₊u_i = -0.060919
busbar₊u_r = 1.0457
ctrld_gen₊avr₊v_fb = 7.2327e-23 (guess 0)
ctrld_gen₊avr₊vfout = 1.2089 (guess 1) (bounds 0.0..Inf)
ctrld_gen₊avr₊vm = 1.0475 (guess 1)
ctrld_gen₊avr₊vr = -0.058633 (guess 0)
ctrld_gen₊gov₊xg1 = 0.25 (guess 1)
ctrld_gen₊gov₊xg2 = 0.25 (guess 0)
ctrld_gen₊gov₊τ_m₊u = 0.25 (guess 0)
ctrld_gen₊machine₊E′_d = 0.041105 (guess 1)
ctrld_gen₊machine₊E′_q = 1.0902 (guess 0)
ctrld_gen₊machine₊δ = 0.084805 (guess 0)
ctrld_gen₊machine₊ψ″_d = 1.0583 (guess 1)
ctrld_gen₊machine₊ψ″_q = -0.12223 (guess 0)
ctrld_gen₊machine₊ω = 1 (guess 1)
Outputs:
busbar₊u_i = -0.060919
busbar₊u_r = 1.0457
Parameters:
ctrld_gen₊avr₊E1 = 3.5461
ctrld_gen₊avr₊E2 = 4.7281
ctrld_gen₊avr₊Ka = 5
ctrld_gen₊avr₊Ke = -0.0485
ctrld_gen₊avr₊Kf = 0.04
ctrld_gen₊avr₊Se1 = 0.08
ctrld_gen₊avr₊Se2 = 0.26
ctrld_gen₊avr₊Ta = 0.06
ctrld_gen₊avr₊Te = 0.25
ctrld_gen₊avr₊Tf = 1
ctrld_gen₊avr₊Tr = 0.01
ctrld_gen₊avr₊vr_max = 1
ctrld_gen₊avr₊vr_min = -1
ctrld_gen₊avr₊vref = 1.0358 (guess 1)
ctrld_gen₊gov₊DT = 0
ctrld_gen₊gov₊R = 0.05
ctrld_gen₊gov₊T1 = 0.5
ctrld_gen₊gov₊T2 = 2.1
ctrld_gen₊gov₊T3 = 7.2
ctrld_gen₊gov₊V_max = 1
ctrld_gen₊gov₊V_min = 0
ctrld_gen₊gov₊p_ref = 0.0125 (guess 1)
ctrld_gen₊gov₊ω_ref = 1
ctrld_gen₊machine₊D = 0
ctrld_gen₊machine₊H = 4.2
ctrld_gen₊machine₊R_s = 0
ctrld_gen₊machine₊S_b = 100
ctrld_gen₊machine₊Sn = 1000
ctrld_gen₊machine₊T′_d0 = 10.2
ctrld_gen₊machine₊T′_q0 = 2
ctrld_gen₊machine₊T″_d0 = 0.05
ctrld_gen₊machine₊T″_q0 = 0.035
ctrld_gen₊machine₊V_b = 16.5
ctrld_gen₊machine₊Vn = 16.5
ctrld_gen₊machine₊X_d = 1
ctrld_gen₊machine₊X_ls = 0.125
ctrld_gen₊machine₊X_q = 0.69
ctrld_gen₊machine₊X′_d = 0.31
ctrld_gen₊machine₊X′_q = 0.5
ctrld_gen₊machine₊X″_d = 0.25
ctrld_gen₊machine₊X″_q = 0.25
ctrld_gen₊machine₊ω_b = 376.99In the state dump we see how the initialization successfully set all the previously unknown values, including the control parameters avr₊vref and gov₊p_ref. We have our first initialized component!
However, in practice it's not always so easy.
Handling Structurally Underconstrained Components
Recalling from Part 1, we have an uncontrolled machine together with a load on bus 39.
╔════════════════════════════════╗
║ Unctr. Ma. Load Bus (compiled) ║
║ ┌────────────────────────┐ ║
Network ║ │MTKBus ┌─────────┐ │ ║
interface ║ │ ┌─┤ Machine │ │ ║
current ────→│ ┌──────┐ │ └─────────┘ │ ║
║ │ │BusBar├─o │ ║
voltage ←────│ └──────┘ │ ┌──────┐ │ ║
║ │ └─┤ Load │ │ ║
║ │ └──────┘ │ ║
║ └────────────────────────┘ ║
╚════════════════════════════════╝If we try to initialize this component as before, we run into a problem:
interf_v39 = Dict(
:busbar₊u_r => 1.01419,
:busbar₊u_i => -0.179795,
:busbar₊i_i => -1.72223,
:busbar₊i_r => 0.720135,
)
initialize_component!(nw[VIndex(39)]; default_overrides=interf_v39)Set additional default for busbar₊u_r: 1.01419
Set additional default for busbar₊u_i: -0.179795
Set additional default for busbar₊i_i: -1.72223
Set additional default for busbar₊i_r: 0.720135
┌ Error: Initialization problem underconstraint. 10 Equations for 11 free variables: [:machine₊terminal₊i_i, :machine₊terminal₊i_r, :machine₊ψ″_q, :machine₊ψ″_d, :machine₊E′_d, :machine₊E′_q, :machine₊ω, :machine₊δ, :machine₊vf_set, :machine₊τ_m_set, :ZIPLoad₊Vset]. Consider passing additional constraints using `InitConstraint`.
└ @ Main ieee39_part2.md:251Even though we set the interface values, the problem is still underconstrained! Let's check the free symbols:
println("free u: ", free_u(nw[VIndex(39)]))
println("free p: ", free_p(nw[VIndex(39)]))free u: [:machine₊terminal₊i_i, :machine₊terminal₊i_r, :machine₊ψ″_q, :machine₊ψ″_d, :machine₊E′_d, :machine₊E′_q, :machine₊ω, :machine₊δ]
free p: [:machine₊vf_set, :machine₊τ_m_set, :ZIPLoad₊Vset]We see 8 free states and 3 free parameters, however we only have 8 state + 2 output equations:
VertexModel :bus39 NoFeedForward() @ Vertex 39
├─ 2 inputs: [busbar₊i_r=0.72013, busbar₊i_i=-1.7222]
├─ 8 states: [machine₊terminal₊i_i≈0, machine₊terminal₊i_r≈0, machine₊ψ″_q≈0, machine₊ψ″_d≈1, machine₊E′_d≈1, machine₊E′_q≈0, machine₊ω≈1, machine₊δ≈0]
| with diagonal mass matrix [0, 0, 1, 1, 1, 1, 1, 1]
├─ 2 outputs: [busbar₊u_r=1.0142, busbar₊u_i=-0.1798]
└─ 30 params: [machine₊vf_set≈1, machine₊τ_m_set≈1, machine₊R_s=0, machine₊X_d=2, machine₊X_q=1.9, machine₊X′_d=0.6, machine₊X′_q=0.8, machine₊X″_d=0.4, machine₊X″_q=0.4, machine₊X_ls=0.3, machine₊T′_d0=7, machine₊T″_d0=0.05, machine₊T′_q0=0.7, machine₊T″_q0=0.035, machine₊H=5, machine₊D=0, machine₊S_b=100, machine₊V_b=345, machine₊ω_b=376.99, machine₊Sn=10000, machine₊Vn=345, ZIPLoad₊Pset=-11.04, ZIPLoad₊Qset=-2.5, ZIPLoad₊Vset≈1, ZIPLoad₊KpZ=1, ZIPLoad₊KqZ=1, ZIPLoad₊KpI=0, ZIPLoad₊KqI=0, ZIPLoad₊KpC=0, ZIPLoad₊KqC=0]
Powerflow model :pvbus with [pv₊P=-1.04, pv₊V=1.03]Even though we have enough set parameters to initialize machine and load on its own, we cannot do it simultaneously. Intuitively speaking, it's just not clear for the solver which of the two components provides how much power.
To solve this, we have essentially 3 methods:
Method 1: Manual setting of defaults
The simplest solution is to manually set more defaults. For example, we know that we want to initialize the ZIP load around the initialization point, i.e. $V_\mathrm{set}$ should be the same as the bus voltage magnitude.
vm_manual = copy(nw[VIndex(39)])
u_r = get_initial_state(vm_manual, :busbar₊u_r)
u_i = get_initial_state(vm_manual, :busbar₊u_i)
set_default!(vm_manual, :ZIPLoad₊Vset, sqrt(u_r^2 + u_i^2))
initialize_component!(vm_manual)[ Info: Apply positivity/negativity conserving variable transformation on [:machine₊vf_set, :machine₊τ_m_set] to satisfy bounds.
[ Info: Initialization problem is fully constrained. Created NonlinearLeastSquaresProblem for [:machine₊terminal₊i_i, :machine₊terminal₊i_r, :machine₊ψ″_q, :machine₊ψ″_d, :machine₊E′_d, :machine₊E′_q, :machine₊ω, :machine₊δ, :machine₊vf_set, :machine₊τ_m_set]
[ Info: Initialization successful with residual 1.2820961574625452e-15The initialization succeeded now!
Note how we can skip the default_overrides keyword argument, since the first (failing) call of initialize_component! already "burned in" the default overrides! Mutating state is a powerful tool, but it needs care!
Method 2: Adding an init_formula
The problem with the previous method is that it is quite manual. In reality, we would never go through this very manual initialization process. The fact that the model is structurally underconstrained is a property of the model and should therefore be handled by the model. To do so, NetworkDynamics.jl provides the InitFormula mechanism.
An InitFormula is a symbolic formula that is evaluated during the initialization process. It is attached to the VertexModel so it can be evaluated automatically during the initialization process. The "formula" we want to apply is simply the equation
\[V_\mathrm{set} = \sqrt{u_r^2 + u_i^2}\]
vm_formula = copy(nw[VIndex(39)])
formula = @initformula :ZIPLoad₊Vset = sqrt(:busbar₊u_r^2 + :busbar₊u_i^2)
set_initformula!(vm_formula, formula)VertexModel :bus39 NoFeedForward() @ Vertex 39
├─ 2 inputs: [busbar₊i_r=0.72013, busbar₊i_i=-1.7222]
├─ 8 states: [machine₊terminal₊i_i≈0, machine₊terminal₊i_r≈0, machine₊ψ″_q≈0, machine₊ψ″_d≈1, machine₊E′_d≈1, machine₊E′_q≈0, machine₊ω≈1, machine₊δ≈0]
| with diagonal mass matrix [0, 0, 1, 1, 1, 1, 1, 1]
├─ 2 outputs: [busbar₊u_r=1.0142, busbar₊u_i=-0.1798]
├─ 30 params: [machine₊vf_set≈1, machine₊τ_m_set≈1, machine₊R_s=0, machine₊X_d=2, machine₊X_q=1.9, machine₊X′_d=0.6, machine₊X′_q=0.8, machine₊X″_d=0.4, machine₊X″_q=0.4, machine₊X_ls=0.3, machine₊T′_d0=7, machine₊T″_d0=0.05, machine₊T′_q0=0.7, machine₊T″_q0=0.035, machine₊H=5, machine₊D=0, machine₊S_b=100, machine₊V_b=345, machine₊ω_b=376.99, machine₊Sn=10000, machine₊Vn=345, ZIPLoad₊Pset=-11.04, ZIPLoad₊Qset=-2.5, ZIPLoad₊Vset≈1, ZIPLoad₊KpZ=1, ZIPLoad₊KqZ=1, ZIPLoad₊KpI=0, ZIPLoad₊KqI=0, ZIPLoad₊KpC=0, ZIPLoad₊KqC=0]
└─ 1 add. init eq. from 1 formula setting [:ZIPLoad₊Vset]
Powerflow model :pvbus with [pv₊P=-1.04, pv₊V=1.03]The printout shows 1 additional initialization equation was attached to the model.
The initialization works now:
initialize_component!(vm_formula)InitFomula: setting default for :ZIPLoad₊Vset to 1.0300036884035901
[ Info: Apply positivity/negativity conserving variable transformation on [:machine₊vf_set, :machine₊τ_m_set] to satisfy bounds.
[ Info: Initialization problem is fully constrained. Created NonlinearLeastSquaresProblem for [:machine₊terminal₊i_i, :machine₊terminal₊i_r, :machine₊ψ″_q, :machine₊ψ″_d, :machine₊E′_d, :machine₊E′_q, :machine₊ω, :machine₊δ, :machine₊vf_set, :machine₊τ_m_set]
[ Info: Initialization successful with residual 1.2820961574625452e-15The init formula is applied early in the initialization process, essentially writing a new default for ZIPLoad₊Vset based on the other defaults.
This reduced the number of free variables to 10, thus the system was solvable.
Method 3: Using an InitConstraint
Sometimes, your additional initialization needs are more complicated. Similar to defining a formula, which is evaluated before the actual initialization, NetworkDynamics provides a mechanism for injecting additional constraints into the initialization process.
In contrast to the formula, the constraint does not need to be explicitly solvable, as it defines a residual equation
\[0 = c(x) = V_\mathrm{set} - \sqrt{u_r^2 + u_i^2}\]
vm_constraint = copy(nw[VIndex(39)])
constraint = @initconstraint begin
:ZIPLoad₊Vset - sqrt(:busbar₊u_r^2 + :busbar₊u_i^2)
end
set_initconstraint!(vm_constraint, constraint)VertexModel :bus39 NoFeedForward() @ Vertex 39
├─ 2 inputs: [busbar₊i_r=0.72013, busbar₊i_i=-1.7222]
├─ 8 states: [machine₊terminal₊i_i≈0, machine₊terminal₊i_r≈0, machine₊ψ″_q≈0, machine₊ψ″_d≈1, machine₊E′_d≈1, machine₊E′_q≈0, machine₊ω≈1, machine₊δ≈0]
| with diagonal mass matrix [0, 0, 1, 1, 1, 1, 1, 1]
├─ 2 outputs: [busbar₊u_r=1.0142, busbar₊u_i=-0.1798]
├─ 30 params: [machine₊vf_set≈1, machine₊τ_m_set≈1, machine₊R_s=0, machine₊X_d=2, machine₊X_q=1.9, machine₊X′_d=0.6, machine₊X′_q=0.8, machine₊X″_d=0.4, machine₊X″_q=0.4, machine₊X_ls=0.3, machine₊T′_d0=7, machine₊T″_d0=0.05, machine₊T′_q0=0.7, machine₊T″_q0=0.035, machine₊H=5, machine₊D=0, machine₊S_b=100, machine₊V_b=345, machine₊ω_b=376.99, machine₊Sn=10000, machine₊Vn=345, ZIPLoad₊Pset=-11.04, ZIPLoad₊Qset=-2.5, ZIPLoad₊Vset≈1, ZIPLoad₊KpZ=1, ZIPLoad₊KqZ=1, ZIPLoad₊KpI=0, ZIPLoad₊KqI=0, ZIPLoad₊KpC=0, ZIPLoad₊KqC=0]
└─ 1 add. init eq. from 1 constraint for [:ZIPLoad₊Vset, :busbar₊u_r, :busbar₊u_i]
Powerflow model :pvbus with [pv₊P=-1.04, pv₊V=1.03]With this added constraint, the initialization process is solvable again, since we now have 11 equations for the 11 free variables.
initialize_component!(vm_constraint)[ Info: Apply positivity/negativity conserving variable transformation on [:machine₊vf_set, :machine₊τ_m_set] to satisfy bounds.
[ Info: Initialization problem is fully constrained. Created NonlinearLeastSquaresProblem for [:machine₊terminal₊i_i, :machine₊terminal₊i_r, :machine₊ψ″_q, :machine₊ψ″_d, :machine₊E′_d, :machine₊E′_q, :machine₊ω, :machine₊δ, :machine₊vf_set, :machine₊τ_m_set, :ZIPLoad₊Vset]
[ Info: Initialization successful with residual 1.3086704628040893e-15For this particular case, method (2) is the way to go. However there are cases where the constraint is more complex and cannot be expressed as a formula.
See NetworkDynamics docs on Advanced Component Initialization: Formulas and Constraints and the PowerDynamics specific extension Advanced Component Initialization for more information on method 2 and 3.
Automatic Initialization of Full Network
Let's return from our excursion into individual component initialization and focus on the whole network again. As we've just seen, we have structurally underconstrained components in the network. Let's define the init formulas for the two buses which have loads and machines:
formula = @initformula :ZIPLoad₊Vset = sqrt(:busbar₊u_r^2 + :busbar₊u_i^2)
set_initformula!(nw[VIndex(31)], formula)
set_initformula!(nw[VIndex(39)], formula)With that, the componentwise initialization of the whole network is possible:
initialize_componentwise!(nw; default_overrides=interf)Even shorter, we can just use initialize_from_pf! to do everything from exporting the power flow model, finding the fixpoint and initializing all components:
s0 = initialize_from_pf!(nw)NWState{Vector{Float64}} of Network (39 vertices, 46 edges)
├─ VIndex(1, :busbar₊u_i) => -0.1537005168018229
├─ VIndex(1, :busbar₊u_r) => 1.0360171229045203
├─ VIndex(2, :busbar₊u_i) => -0.10513936683885863
├─ VIndex(2, :busbar₊u_r) => 1.0434526576051297
├─ VIndex(5, :busbar₊u_i) => -0.15053530762365022
├─ VIndex(5, :busbar₊u_r) => 0.993976368877831
├─ VIndex(6, :busbar₊u_i) => -0.1393645217247641
├─ VIndex(6, :busbar₊u_r) => 0.9979886162106123
⋮
├─ VIndex(39, :machine₊terminal₊i_i) => -2.538683831099698
├─ VIndex(39, :machine₊terminal₊i_r) => 9.410063917554226
├─ VIndex(39, :machine₊ψ″_q) => -0.1505993365374145
├─ VIndex(39, :machine₊ψ″_d) => 1.0219450150609855
├─ VIndex(39, :machine₊E′_d) => 0.10353704204613715
├─ VIndex(39, :machine₊E′_q) => 1.029534387321863
├─ VIndex(39, :machine₊ω) => 1.0
└─ VIndex(39, :machine₊δ) => -0.0009449613162717396
p = NWParameter([-3.22, -0.024, 1.03017, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, -5.0 … 0.97561, 1.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0])
t = nothingThis page was generated using Literate.jl.