IEEE39 Bus Tutorial - Part II: Initialization
This tutorial can be downloaded as a normal Julia script here.
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 192 states and 1662 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 1200 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]
├─ 14 states: [ctrld_gen₊gov₊xg2≈0, ctrld_gen₊gov₊xg1≈1, ctrld_gen₊avr₊v_fb≈0, ctrld_gen₊avr₊vfout≈1, 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≈0, ctrld_gen₊machine₊E′_q≈1, ctrld_gen₊machine₊ω≈1, ctrld_gen₊machine₊δ≈0, busbar₊u_r=1, busbar₊u_i=0]
| with diagonal mass matrix [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0]
├─ 2 outputs: [busbar₊u_r=1, busbar₊u_i=0]
├─ 43 params: [busbar₊Vbase=16.5, systembase₊Sbase=100, systembase₊ωbase=376.99, systembase₊ωframe=1, 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₊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]
└─ 2 InitFormula clusters:
explicitly initializes 1/16 variables from seeds [1 already set param]
explicitly initializes 1/16 variables from seeds [1 already set param]
Powerflow model :pvbus with [busbar₊Vbase=1, systembase₊Sbase=100, systembase₊ωbase=376.99, systembase₊ωframe=1, 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=-2.3866, busbar₊i_i=0]
├─ 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]
└─ 6 params: [busbar₊Vbase=1, systembase₊Sbase=100, systembase₊ωbase=376.99, systembase₊ωframe=1, 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 [compile_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]
├─ 21 params: [src₊Vbase, dst₊Vbase, systembase₊Sbase=100, systembase₊ωbase=376.99, systembase₊ωframe=1, 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]
└─ 2 params default_from: :src₊Vbase ← src busbar₊Vbase, :dst₊Vbase ← dst busbar₊VbaseThis 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_r) => 1.0
├─ VIndex(1, :busbar₊u_i) => 0.0
├─ VIndex(2, :busbar₊u_r) => 1.0
├─ VIndex(2, :busbar₊u_i) => 0.0
├─ VIndex(3, :busbar₊u_r) => 1.0
├─ VIndex(3, :busbar₊u_i) => 0.0
├─ VIndex(4, :busbar₊u_r) => 1.0
├─ VIndex(4, :busbar₊u_i) => 0.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([1.0, 100.0, 376.991, 1.0, 0.0, 0.0, 1.0, 100.0, 376.991, 1.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_r) => 1.0360171229045212
├─ VIndex(1, :busbar₊u_i) => -0.1537005168018203
├─ VIndex(2, :busbar₊u_r) => 1.0434526576051306
├─ VIndex(2, :busbar₊u_i) => -0.10513936683885611
├─ VIndex(3, :busbar₊u_r) => 1.018593888445931
├─ VIndex(3, :busbar₊u_i) => -0.1540214604217488
├─ VIndex(4, :busbar₊u_r) => 0.9897852075420123
├─ VIndex(4, :busbar₊u_i) => -0.16752811347984162
⋮
├─ VIndex(36, :busbar₊u_r) => 1.0522992315769322
├─ VIndex(36, :busbar₊u_i) => 0.15394374497652874
├─ VIndex(37, :busbar₊u_r) => 1.0268824837830113
├─ VIndex(37, :busbar₊u_i) => 0.04341808836807661
├─ VIndex(38, :busbar₊u_r) => 1.016983848045741
├─ VIndex(38, :busbar₊u_i) => 0.13944918574007684
├─ VIndex(39, :busbar₊u_r) => 1.0141861967677503
└─ VIndex(39, :busbar₊u_i) => -0.1797951594139676
p = NWParameter([1.0, 100.0, 376.991, 1.0, 0.0, 0.0, 1.0, 100.0, 376.991, 1.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.SymbolicIndex{Int64, Symbol}, Float64} with 524 entries:
VIndex(1, :busbar₊i_r) => -4.44089e-16
VIndex(1, :busbar₊i_i) => 8.32667e-16
VIndex(1, :busbar₊u_r) => 1.03602
VIndex(1, :busbar₊u_i) => -0.153701
VIndex(2, :busbar₊i_r) => 3.9968e-15
VIndex(2, :busbar₊i_i) => 1.53211e-14
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) => 8.43769e-15
VIndex(5, :busbar₊i_i) => 1.4877e-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]
├─ 14 states: [ctrld_gen₊gov₊xg2≈0, ctrld_gen₊gov₊xg1≈1, ctrld_gen₊avr₊v_fb≈0, ctrld_gen₊avr₊vfout≈1, 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≈0, ctrld_gen₊machine₊E′_q≈1, ctrld_gen₊machine₊ω≈1, ctrld_gen₊machine₊δ≈0, busbar₊u_r=1, busbar₊u_i=0]
| with diagonal mass matrix [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0]
├─ 2 outputs: [busbar₊u_r=1, busbar₊u_i=0]
├─ 43 params: [busbar₊Vbase=16.5, systembase₊Sbase=100, systembase₊ωbase=376.99, systembase₊ωframe=1, 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₊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]
└─ 2 InitFormula clusters:
explicitly initializes 1/16 variables from seeds [1 already set param]
explicitly initializes 1/16 variables from seeds [1 already set param]
Powerflow model :pvbus with [busbar₊Vbase=1, systembase₊Sbase=100, systembase₊ωbase=376.99, systembase₊ωframe=1, 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 = 0
busbar₊u_r = 1
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₊machine₊E′_d = uninitialized (guess 0)
ctrld_gen₊machine₊E′_q = uninitialized (guess 1)
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 = 0
busbar₊u_r = 1
Parameters:
busbar₊Vbase = 16.5
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₊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₊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
systembase₊Sbase = 100
systembase₊ωbase = 376.99
systembase₊ωframe = 1Right 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)Dict{Symbol, Float64} with 59 entries:
:systembase₊ωframe => 1.0
:ctrld_gen₊machine₊E′_q => 0.999701
:ctrld_gen₊avr₊Tr => 0.01
:ctrld_gen₊gov₊V_min => 0.0
:ctrld_gen₊machine₊D => 0.0
:ctrld_gen₊gov₊V_max => 1.0
:ctrld_gen₊avr₊Se2 => 0.26
:ctrld_gen₊machine₊T″_d0 => 0.05
:systembase₊ωbase => 376.991
:ctrld_gen₊avr₊Ka => 5.0
:ctrld_gen₊machine₊X_ls => 0.125
:ctrld_gen₊avr₊E2 => 4.72813
:ctrld_gen₊avr₊Tf => 1.0
:ctrld_gen₊avr₊Ke => -0.0485
:ctrld_gen₊gov₊xg1 => -0.0519183
:ctrld_gen₊gov₊xg2 => -0.0519183
:ctrld_gen₊gov₊p_ref => -0.00259591
:ctrld_gen₊avr₊E1 => 3.5461
:ctrld_gen₊machine₊R_s => 0.0
⋮ => ⋮Wait! 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) - De-aliased guesses (2 conflicting dropped):
- :busbar₊terminal₊i_i ⇒ :busbar₊i_i ( 0)
- :busbar₊terminal₊i_r ⇒ :busbar₊i_r ( 0)
- :ctrld_gen₊avr₊v_mag₊u ⇒ :ctrld_gen₊machine₊v_mag ( 0)
- :ctrld_gen₊gov₊τ_m₊u ⇒ :ctrld_gen₊machine₊τ_m ( 0)
- :ctrld_gen₊gov₊ω_meas₊u ⇒ :ctrld_gen₊machine₊ω ( 0) dropped, class holds 1
- :ctrld_gen₊machine₊Pout₊u ⇒ :ctrld_gen₊machine₊P ( 0)
- :ctrld_gen₊machine₊Qout₊u ⇒ :ctrld_gen₊machine₊Q ( 0)
- :ctrld_gen₊machine₊ωout₊u ⇒ :ctrld_gen₊machine₊ω ( 0) dropped, class holds 1
- InitFormula: ctrld_gen₊machine₊Sn = ctrld_gen₊machine₊Sbase yields to existing default :ctrld_gen₊machine₊Sn = 1000.0
- InitFormula: ctrld_gen₊machine₊Vn = ctrld_gen₊machine₊Vbase yields to existing default :ctrld_gen₊machine₊Vn = 16.5
- Apply positivity/negativity conserving variable transformation on [:ctrld_gen₊avr₊vfout] to satisfy bounds.
- Initialization problem is overconstrained (14 vars for 16 equations). Created NonlinearLeastSquaresProblem for:
- :ctrld_gen₊gov₊xg2 (guess= 0)
- :ctrld_gen₊gov₊xg1 (guess= 1)
- :ctrld_gen₊avr₊v_fb (guess= 0)
- :ctrld_gen₊avr₊vfout (guess= 1)
- :ctrld_gen₊avr₊vr (guess= 0)
- :ctrld_gen₊avr₊vm (guess= 1)
- :ctrld_gen₊machine₊ψ″_q (guess= 0)
- :ctrld_gen₊machine₊ψ″_d (guess= 1)
- :ctrld_gen₊machine₊E′_d (guess= 0)
- :ctrld_gen₊machine₊E′_q (guess= 1)
- :ctrld_gen₊machine₊ω (guess= 1)
- :ctrld_gen₊machine₊δ (guess= 0)
- :ctrld_gen₊avr₊vref (guess= 1)
- :ctrld_gen₊gov₊p_ref (guess= 1)
- Initialization successful with residual 1.5855369605350334e-14
- :ctrld_gen₊gov₊xg2 => 0.25
- :ctrld_gen₊gov₊xg1 => 0.25
- :ctrld_gen₊avr₊v_fb => -1.7638e-24
- :ctrld_gen₊avr₊vfout => 1.2089
- :ctrld_gen₊avr₊vr => -0.058633
- :ctrld_gen₊avr₊vm => 1.0475
- :ctrld_gen₊machine₊ψ″_q => -0.12223
- :ctrld_gen₊machine₊ψ″_d => 1.0583
- :ctrld_gen₊machine₊E′_d => 0.041105
- :ctrld_gen₊machine₊E′_q => 1.0902
- :ctrld_gen₊machine₊ω => 1
- :ctrld_gen₊machine₊δ => 0.084805
- :ctrld_gen₊avr₊vref => 1.0358
- :ctrld_gen₊gov₊p_ref => 0.0125The 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) - Additional defaults:
- :busbar₊u_r = 1.0457
- :busbar₊u_i = -0.060919
- :busbar₊i_i = 1.5317
- :busbar₊i_r = -2.3014
- De-aliased guesses (2 conflicting dropped):
- :busbar₊terminal₊i_i ⇒ :busbar₊i_i ( 0)
- :busbar₊terminal₊i_r ⇒ :busbar₊i_r ( 0)
- :ctrld_gen₊avr₊v_mag₊u ⇒ :ctrld_gen₊machine₊v_mag ( 0)
- :ctrld_gen₊gov₊τ_m₊u ⇒ :ctrld_gen₊machine₊τ_m ( 0)
- :ctrld_gen₊gov₊ω_meas₊u ⇒ :ctrld_gen₊machine₊ω ( 0) dropped, class holds 1
- :ctrld_gen₊machine₊Pout₊u ⇒ :ctrld_gen₊machine₊P ( 0)
- :ctrld_gen₊machine₊Qout₊u ⇒ :ctrld_gen₊machine₊Q ( 0)
- :ctrld_gen₊machine₊ωout₊u ⇒ :ctrld_gen₊machine₊ω ( 0) dropped, class holds 1
- InitFormula: ctrld_gen₊machine₊Sn = ctrld_gen₊machine₊Sbase yields to existing default :ctrld_gen₊machine₊Sn = 1000.0
- InitFormula: ctrld_gen₊machine₊Vn = ctrld_gen₊machine₊Vbase yields to existing default :ctrld_gen₊machine₊Vn = 16.5
- Apply positivity/negativity conserving variable transformation on [:ctrld_gen₊avr₊vfout] to satisfy bounds.
- Initialization problem is overconstrained (14 vars for 16 equations). Created NonlinearLeastSquaresProblem for:
- :ctrld_gen₊gov₊xg2 (guess= 0)
- :ctrld_gen₊gov₊xg1 (guess= 1)
- :ctrld_gen₊avr₊v_fb (guess= 0)
- :ctrld_gen₊avr₊vfout (guess= 1)
- :ctrld_gen₊avr₊vr (guess= 0)
- :ctrld_gen₊avr₊vm (guess= 1)
- :ctrld_gen₊machine₊ψ″_q (guess= 0)
- :ctrld_gen₊machine₊ψ″_d (guess= 1)
- :ctrld_gen₊machine₊E′_d (guess= 0)
- :ctrld_gen₊machine₊E′_q (guess= 1)
- :ctrld_gen₊machine₊ω (guess= 1)
- :ctrld_gen₊machine₊δ (guess= 0)
- :ctrld_gen₊avr₊vref (guess= 1)
- :ctrld_gen₊gov₊p_ref (guess= 1)
- Initialization successful with residual 1.5855369605350334e-14
- :ctrld_gen₊gov₊xg2 => 0.25
- :ctrld_gen₊gov₊xg1 => 0.25
- :ctrld_gen₊avr₊v_fb => -1.7638e-24
- :ctrld_gen₊avr₊vfout => 1.2089
- :ctrld_gen₊avr₊vr => -0.058633
- :ctrld_gen₊avr₊vm => 1.0475
- :ctrld_gen₊machine₊ψ″_q => -0.12223
- :ctrld_gen₊machine₊ψ″_d => 1.0583
- :ctrld_gen₊machine₊E′_d => 0.041105
- :ctrld_gen₊machine₊E′_q => 1.0902
- :ctrld_gen₊machine₊ω => 1
- :ctrld_gen₊machine₊δ => 0.084805
- :ctrld_gen₊avr₊vref => 1.0358
- :ctrld_gen₊gov₊p_ref => 0.0125
Inputs:
busbar₊i_i = 1.5317 (guess 0)
busbar₊i_r = -2.3014 (guess 0)
States:
busbar₊u_i = -0.060919
busbar₊u_r = 1.0457
ctrld_gen₊avr₊v_fb = -1.7638e-24 (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₊machine₊E′_d = 0.041105 (guess 0)
ctrld_gen₊machine₊E′_q = 1.0902 (guess 1)
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:
busbar₊Vbase = 16.5
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₊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₊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
systembase₊Sbase = 100
systembase₊ωbase = 376.99
systembase₊ωframe = 1In 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)VertexModel :bus39 PureStateMap() @ Vertex 39
├─ 2 inputs: [busbar₊i_r=0.72013, busbar₊i_i=-1.7222]
├─ 8 states: [machine₊ψ″_q=-0.042639, machine₊ψ″_d=1.0269, machine₊E′_d=0.029315, machine₊E′_q=1.0251, machine₊ω=1, machine₊δ=-0.12628, busbar₊u_r=1.0142, busbar₊u_i=-0.1798]
| with diagonal mass matrix [1, 1, 1, 1, 1, 1, 0, 0]
├─ 2 outputs: [busbar₊u_r=1.0142, busbar₊u_i=-0.1798]
├─ 31 params: [busbar₊Vbase=345, systembase₊Sbase=100, systembase₊ωbase=376.99, systembase₊ωframe=1, machine₊vf_set=1.0165, machine₊τ_m_set=0.027104, 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₊Sn=10000, machine₊Vn=345, ZIPLoad₊Pset=-11.04, ZIPLoad₊Qset=-2.5, ZIPLoad₊Vset=1.7672, ZIPLoad₊KpZ=1, ZIPLoad₊KqZ=1, ZIPLoad₊KpI=0, ZIPLoad₊KqI=0, ZIPLoad₊KpC=0, ZIPLoad₊KqC=0]
└─ 2 InitFormula clusters:
explicitly initializes 1/2 variables from seeds [1 already set param]
explicitly initializes 1/2 variables from seeds [1 already set param]
Powerflow model :pvbus with [busbar₊Vbase=1, systembase₊Sbase=100, systembase₊ωbase=376.99, systembase₊ωframe=1, pv₊P=-1.04, pv₊V=1.03]Even 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₊ψ″_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 PureStateMap() @ Vertex 39
├─ 2 inputs: [busbar₊i_r=0.72013, busbar₊i_i=-1.7222]
├─ 8 states: [machine₊ψ″_q=-0.042639, machine₊ψ″_d=1.0269, machine₊E′_d=0.029315, machine₊E′_q=1.0251, machine₊ω=1, machine₊δ=-0.12628, busbar₊u_r=1.0142, busbar₊u_i=-0.1798]
| with diagonal mass matrix [1, 1, 1, 1, 1, 1, 0, 0]
├─ 2 outputs: [busbar₊u_r=1.0142, busbar₊u_i=-0.1798]
├─ 31 params: [busbar₊Vbase=345, systembase₊Sbase=100, systembase₊ωbase=376.99, systembase₊ωframe=1, machine₊vf_set=1.0165, machine₊τ_m_set=0.027104, 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₊Sn=10000, machine₊Vn=345, ZIPLoad₊Pset=-11.04, ZIPLoad₊Qset=-2.5, ZIPLoad₊Vset=1.7672, ZIPLoad₊KpZ=1, ZIPLoad₊KqZ=1, ZIPLoad₊KpI=0, ZIPLoad₊KqI=0, ZIPLoad₊KpC=0, ZIPLoad₊KqC=0]
└─ 2 InitFormula clusters:
explicitly initializes 1/2 variables from seeds [1 already set param]
explicitly initializes 1/2 variables from seeds [1 already set param]
Powerflow model :pvbus with [busbar₊Vbase=1, systembase₊Sbase=100, systembase₊ωbase=376.99, systembase₊ωframe=1, 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) - De-aliased guesses (1 conflicting dropped):
- :machine₊ωout₊u ⇒ :machine₊ω ( 0) dropped, class holds 1
- InitFormula: machine₊Sn = machine₊Sbase yields to existing default :machine₊Sn = 10000.0
- InitFormula: machine₊Vn = machine₊Vbase yields to existing default :machine₊Vn = 345.0
- Apply positivity/negativity conserving variable transformation on [:machine₊vf_set, :machine₊τ_m_set] to satisfy bounds.
- Initialization problem is overconstrained (8 vars for 10 equations). Created NonlinearLeastSquaresProblem for:
- :machine₊ψ″_q (guess= 0)
- :machine₊ψ″_d (guess= 1)
- :machine₊E′_d (guess= 0)
- :machine₊E′_q (guess= 1)
- :machine₊ω (guess= 1)
- :machine₊δ (guess= 0)
- :machine₊vf_set (guess= 1)
- :machine₊τ_m_set (guess= 1)
- Initialization successful with residual 1.2953324773051273e-15
- :machine₊ψ″_q => -0.1506
- :machine₊ψ″_d => 1.0219
- :machine₊E′_d => 0.10354
- :machine₊E′_q => 1.0295
- :machine₊ω => 1
- :machine₊δ => -0.0009454
- :machine₊vf_set => 1.065
- :machine₊τ_m_set => 0.1The 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 PureStateMap() @ Vertex 39
├─ 2 inputs: [busbar₊i_r=0.72013, busbar₊i_i=-1.7222]
├─ 8 states: [machine₊ψ″_q=-0.042639, machine₊ψ″_d=1.0269, machine₊E′_d=0.029315, machine₊E′_q=1.0251, machine₊ω=1, machine₊δ=-0.12628, busbar₊u_r=1.0142, busbar₊u_i=-0.1798]
| with diagonal mass matrix [1, 1, 1, 1, 1, 1, 0, 0]
├─ 2 outputs: [busbar₊u_r=1.0142, busbar₊u_i=-0.1798]
├─ 31 params: [busbar₊Vbase=345, systembase₊Sbase=100, systembase₊ωbase=376.99, systembase₊ωframe=1, machine₊vf_set=1.0165, machine₊τ_m_set=0.027104, 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₊Sn=10000, machine₊Vn=345, ZIPLoad₊Pset=-11.04, ZIPLoad₊Qset=-2.5, ZIPLoad₊Vset=1.7672, ZIPLoad₊KpZ=1, ZIPLoad₊KqZ=1, ZIPLoad₊KpI=0, ZIPLoad₊KqI=0, ZIPLoad₊KpC=0, ZIPLoad₊KqC=0]
└─ 1 InitFormula cluster:
explicitly initializes 1/1 variables from seeds [busbar₊u_i, busbar₊u_r]
Powerflow model :pvbus with [busbar₊Vbase=1, systembase₊Sbase=100, systembase₊ωbase=376.99, systembase₊ωframe=1, 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) - De-aliased guesses (1 conflicting dropped):
- :machine₊ωout₊u ⇒ :machine₊ω ( 0) dropped, class holds 1
- InitFormulas set:
- :ZIPLoad₊Vset = 1.03 (via ZIPLoad₊Vset = sqrt(busbar₊u_r ^ 2 + busbar₊u_i ^ 2))
- Apply positivity/negativity conserving variable transformation on [:machine₊vf_set, :machine₊τ_m_set] to satisfy bounds.
- Initialization problem is overconstrained (8 vars for 10 equations). Created NonlinearLeastSquaresProblem for:
- :machine₊ψ″_q (guess= 0)
- :machine₊ψ″_d (guess= 1)
- :machine₊E′_d (guess= 0)
- :machine₊E′_q (guess= 1)
- :machine₊ω (guess= 1)
- :machine₊δ (guess= 0)
- :machine₊vf_set (guess= 1)
- :machine₊τ_m_set (guess= 1)
- Initialization successful with residual 1.2953324773051273e-15
- :machine₊ψ″_q => -0.1506
- :machine₊ψ″_d => 1.0219
- :machine₊E′_d => 0.10354
- :machine₊E′_q => 1.0295
- :machine₊ω => 1
- :machine₊δ => -0.0009454
- :machine₊vf_set => 1.065
- :machine₊τ_m_set => 0.1The 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 PureStateMap() @ Vertex 39
├─ 2 inputs: [busbar₊i_r=0.72013, busbar₊i_i=-1.7222]
├─ 8 states: [machine₊ψ″_q=-0.042639, machine₊ψ″_d=1.0269, machine₊E′_d=0.029315, machine₊E′_q=1.0251, machine₊ω=1, machine₊δ=-0.12628, busbar₊u_r=1.0142, busbar₊u_i=-0.1798]
| with diagonal mass matrix [1, 1, 1, 1, 1, 1, 0, 0]
├─ 2 outputs: [busbar₊u_r=1.0142, busbar₊u_i=-0.1798]
├─ 31 params: [busbar₊Vbase=345, systembase₊Sbase=100, systembase₊ωbase=376.99, systembase₊ωframe=1, machine₊vf_set=1.0165, machine₊τ_m_set=0.027104, 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₊Sn=10000, machine₊Vn=345, ZIPLoad₊Pset=-11.04, ZIPLoad₊Qset=-2.5, ZIPLoad₊Vset=1.7672, ZIPLoad₊KpZ=1, ZIPLoad₊KqZ=1, ZIPLoad₊KpI=0, ZIPLoad₊KqI=0, ZIPLoad₊KpC=0, ZIPLoad₊KqC=0]
├─ 2 InitFormula clusters:
| explicitly initializes 1/2 variables from seeds [1 already set param]
| explicitly initializes 1/2 variables from seeds [1 already set param]
└─ 1 add. init eq. from 1 constraint for [:ZIPLoad₊Vset, :busbar₊u_r, :busbar₊u_i]
Powerflow model :pvbus with [busbar₊Vbase=1, systembase₊Sbase=100, systembase₊ωbase=376.99, systembase₊ωframe=1, 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) - De-aliased guesses (1 conflicting dropped):
- :machine₊ωout₊u ⇒ :machine₊ω ( 0) dropped, class holds 1
- InitFormula: machine₊Sn = machine₊Sbase yields to existing default :machine₊Sn = 10000.0
- InitFormula: machine₊Vn = machine₊Vbase yields to existing default :machine₊Vn = 345.0
- Apply positivity/negativity conserving variable transformation on [:machine₊vf_set, :machine₊τ_m_set] to satisfy bounds.
- Initialization problem is overconstrained (9 vars for 11 equations). Created NonlinearLeastSquaresProblem for:
- :machine₊ψ″_q (guess= 0)
- :machine₊ψ″_d (guess= 1)
- :machine₊E′_d (guess= 0)
- :machine₊E′_q (guess= 1)
- :machine₊ω (guess= 1)
- :machine₊δ (guess= 0)
- :machine₊vf_set (guess= 1)
- :machine₊τ_m_set (guess= 1)
- :ZIPLoad₊Vset (guess= 1)
- Initialization successful with residual 1.2953395488600261e-15
- :machine₊ψ″_q => -0.1506
- :machine₊ψ″_d => 1.0219
- :machine₊E′_d => 0.10354
- :machine₊E′_q => 1.0295
- :machine₊ω => 1
- :machine₊δ => -0.0009454
- :machine₊vf_set => 1.065
- :machine₊τ_m_set => 0.1
- :ZIPLoad₊Vset => 1.03For 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.
We start by with solve_powerflow, which extracts the powerflow models from each component, builds the (static) power flow network and solves the problem using a root-finding algorithm.
pfs = solve_powerflow(nw)Next, we get the voltages and the currents at the buses from the powerflow state. Those interface_values are the shared state between the power flow model and the full dynamic model.
interf = interface_values(pfs)Lastly, we initialize the components one by one around the state prescribed by the power flow solution (i.e. we initialize the components, such that their steady state reproduces the bus voltages and currents from the power flow solution). For that we use initialize_componentwise!.
initialize_componentwise!(nw; default_overrides=interf)Since this init flow is so typical, there is a power dynamics specific wrapper around it. We can just use initialize_from_pf! to do everything from exporting the power flow model, finding the fixpoint and initializing all components based on the PF solution:
s0 = initialize_from_pf!(nw)NWState{Vector{Float64}} of Network (39 vertices, 46 edges)
├─ VIndex(1, :busbar₊u_r) => 1.0360171229045212
├─ VIndex(1, :busbar₊u_i) => -0.1537005168018203
├─ VIndex(2, :busbar₊u_r) => 1.0434526576051306
├─ VIndex(2, :busbar₊u_i) => -0.10513936683885611
├─ VIndex(5, :busbar₊u_r) => 0.9939763688778325
├─ VIndex(5, :busbar₊u_i) => -0.15053530762364867
├─ VIndex(6, :busbar₊u_r) => 0.9979886162106136
├─ VIndex(6, :busbar₊u_i) => -0.13936452172476257
⋮
├─ VIndex(39, :machine₊ψ″_q) => -0.15059933653741472
├─ VIndex(39, :machine₊ψ″_d) => 1.0219450150609855
├─ VIndex(39, :machine₊E′_d) => 0.1035370420461373
├─ VIndex(39, :machine₊E′_q) => 1.029534387321863
├─ VIndex(39, :machine₊ω) => 1.0
├─ VIndex(39, :machine₊δ) => -0.000944961316268967
├─ VIndex(39, :busbar₊u_r) => 1.0141861967677503
└─ VIndex(39, :busbar₊u_i) => -0.1797951594139676
p = NWParameter([345.0, 100.0, 376.991, 1.0, 345.0, 100.0, 376.991, 1.0, 345.0, 100.0 … 0.97561, 1.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0])
t = nothingVoltage Bases of the Lines
Initialization is also where the lines learn their voltage bases. In Part I we set busbar₊Vbase on every bus from bus.csv, but we never touched the lines. A LineEnd does not carry a value of its own: it declares that its Vbase is inherited from the bus on that side of the line, and the inheritance is resolved here, when the line is initialized in the context of its neighbours.
The effect is easiest to see on a transformer, whose two ends sit at genuinely different voltage levels — branch 41 connects the 345 kV grid (bus 25) to a 16.5 kV generator bus (37). Lets check the voltage bases for both ends of that line:
s0.e(25=>37, r"Vbase$")FilteringProxy for NWState()
Component filter: EIndex(25 => 37)
State filter: r"Vbase$"
Types: states ✓ parameters ✓ inputs ✓ outputs ✓ observables ✓
Matching Indices:
╭ EIndex(41, :src₊Vbase) 345 :piline_template
╰ EIndex(41, :dst₊Vbase) 16.5 To cross check, lets check the voltage base of the two buses themselves:
s0.v([25,37], :busbar₊Vbase)FilteringProxy for NWState()
Component filter: VIndex([25, 37])
State filter: :busbar₊Vbase
Types: states ✓ parameters ✓ inputs ✓ outputs ✓ observables ✓
Matching Indices:
● VIndex(25, :busbar₊Vbase) 345 :bus25
● VIndex(37, :busbar₊Vbase) 16.5 :bus37Only the SI observables depend on this — the dynamics are formulated in per unit throughout — but it means u_kV, P_MW and friends report the right numbers on both sides without anyone having to keep a table of line voltage levels in sync with the buses.
This page was generated using Literate.jl.