Tutorial on Stepwise Initialization of a Complex Model
This example demonstrates how to initialize a complex network model with both static and dynamic components. The models are closely related to the ones used in the gas network example, but greatly simplified for the sake of this tutorial. We'll create a gas network model with three nodes and pipes connecting them, and show how to:
- Create static models for initialization
- Find a steady-state solution
- Create corresponding dynamic models
- Initialize the dynamic models with the steady-state solution
- Simulate the system with dynamic behavior
This script can be downloaded as a normal Julia script here.
First, let's import the necessary packages:
using NetworkDynamics
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEqTsit5
using CairoMakieNode Models
We'll start by defining our node models using ModelingToolkit. First, let's create a template for common states and equations in all gas nodes:
@mtkmodel GasNode begin
@variables begin
p(t), [description="Pressure"] # node output
q̃_nw(t), [description="aggregated flow from pipes into node"] # node input
q̃_inj(t), [description="flow injected into the network"]
end
@equations begin
q̃_inj ~ -q̃_nw
end
endNow we'll define three specific node types:
A) A constant pressure node that forces pressure to maintain a specific value
@mtkmodel ConstantPressureNode begin
@extend GasNode()
@parameters begin
p_set, [description="Constant pressure setpoint"]
end
@equations begin
p ~ p_set
end
endB) A static prosumer node which forces a certain flow (pressure is fully implicit)
@mtkmodel StaticProsumerNode begin
@extend GasNode()
@parameters begin
q̃_prosumer, [description="flow injected by prosumer"]
end
@equations begin
-q̃_nw ~ q̃_prosumer
end
endC) A dynamic prosumer node with compliance, which introduces dynamic behavior to the pressure state
@mtkmodel DynamicProsumerNode begin
@extend GasNode()
@parameters begin
q̃_prosumer, [description="flow injected by prosumer"]
C=0.1, [description="Compliance"]
end
@equations begin
C*D(p) ~ q̃_prosumer + q̃_nw
end
endD) A pressure control node that tries to maintain a set pressure by adjusting its injection
@mtkmodel PressureControlNode begin
@extend GasNode()
@parameters begin
p_set, [description="Pressure setpoint", guess=1]
K_p=1, [description="Proportional gain"]
K_i=1, [description="Integral gain"]
C=0.1, [description="Compliance"]
end
@variables begin
Δp(t), [description="Pressure error"]
ξ(t), [description="Integral state", guess=0]
q̃_prosumer(t), [description="flow injected by producer"]
end
@equations begin
Δp ~ p_set - p
D(ξ) ~ Δp
q̃_prosumer ~ K_p*Δp + K_i*ξ
C*D(p) ~ q̃_prosumer + q̃_nw
end
endEdge Models
Now we'll define our edge models, starting with a template for the pipe:
@mtkmodel GasPipe begin
@variables begin
q̃(t), [description="flow through pipe"] #output
p_src(t), [description="pressure at start of pipe"] #input
p_dst(t), [description="pressure at end of pipe"] #input
end
endNext, we define a dynamic pipe with inertia (a simple delayed model):
@mtkmodel DynamicPipe begin
@extend GasPipe()
@parameters begin
R=0.1, [description="Resistance"]
M=0.1, [description="Inertia"]
end
@equations begin
M*D(q̃) ~ (p_src - p_dst)/R - q̃ # some simple delayed model
end
endAnd finally, a quasistatic pipe model for initialization purposes. This model equals the dynamic model at steady state, making it ideal for finding initial conditions:
@mtkmodel QuasistaticPipe begin
@extend GasPipe()
@parameters begin
R=0.1, [description="Resistance"]
end
@equations begin
q̃ ~ (p_src - p_dst)/R
end
endDefining a Static Model for Initialization
Our first step is to define a static model that we'll use to find the steady-state solution. This is a crucial step for initializing complex dynamic models.
Step 1: Define all the components of our static model First, node 1 is our producer which will later be a controlled producer. For initialization, we use a static model:
@named v1_mod_static = ConstantPressureNode(p_set=1)
v1_static = VertexModel(v1_mod_static, [:q̃_nw], [:p], vidx=1)
# Nodes 2 and 3 are consumers. For them, we'll use static prosumer models:
@named v2_mod_static = StaticProsumerNode(q̃_prosumer=-0.6) # consumer
v2_static = VertexModel(v2_mod_static, [:q̃_nw], [:p], vidx=2)
@named v3_mod_static = StaticProsumerNode(q̃_prosumer=-0.4) # consumer
v3_static = VertexModel(v3_mod_static, [:q̃_nw], [:p], vidx=3)Now we define the static pipe models connecting our nodes:
@named p_mod_static = QuasistaticPipe()
p12_static = EdgeModel(p_mod_static, [:p_src], [:p_dst], AntiSymmetric([:q̃]), src=1, dst=2)
p13_static = EdgeModel(p_mod_static, [:p_src], [:p_dst], AntiSymmetric([:q̃]), src=1, dst=3)
p23_static = EdgeModel(p_mod_static, [:p_src], [:p_dst], AntiSymmetric([:q̃]), src=2, dst=3)Assemble all components into a static network:
nw_static = Network([v1_static, v2_static, v3_static], [p12_static, p13_static, p23_static])Network with 2 states and 6 parameters
├─ 3 vertices (2 unique types)
└─ 3 edges (1 unique type)
Edge-Aggregation using SequentialAggregator(+)Create an initial guess for the steady state and modify it with reasonable values:
u_static_guess = NWState(nw_static)
u_static_guess.v[2, :p] = 1.0
u_static_guess.v[3, :p] = 1.0Find the steady-state solution using our initial guess:
u_static = find_fixpoint(nw_static, u_static_guess)NWState{Vector{Float64}} of Network (3 vertices, 3 edges)
├─ VIndex(2, :p) => 0.9466666666666667
└─ VIndex(3, :p) => 0.9533333333333334
p = NWParameter([1.0, -0.6, -0.4, 0.1, 0.1, 0.1])
t = nothingDefining a Dynamic Model
Now we'll define our dynamic model using more complex components:
@named v1_mod_dyn = PressureControlNode()
v1_dyn = VertexModel(v1_mod_dyn, [:q̃_nw], [:p], vidx=1)
@named v2_mod_dyn = DynamicProsumerNode(q̃_prosumer=-0.6)
v2_dyn = VertexModel(v2_mod_dyn, [:q̃_nw], [:p], vidx=2)
@named v3_mod_dyn = DynamicProsumerNode(q̃_prosumer=-0.4)
v3_dyn = VertexModel(v3_mod_dyn, [:q̃_nw], [:p], vidx=3)Create dynamic pipe models with inertia:
@named p_mod_dyn = DynamicPipe()
p12_dyn = EdgeModel(p_mod_dyn, [:p_src], [:p_dst], AntiSymmetric([:q̃]), src=1, dst=2)
p13_dyn = EdgeModel(p_mod_dyn, [:p_src], [:p_dst], AntiSymmetric([:q̃]), src=1, dst=3)
p23_dyn = EdgeModel(p_mod_dyn, [:p_src], [:p_dst], AntiSymmetric([:q̃]), src=2, dst=3)Assemble the dynamic network:
nw_dyn = Network([v1_dyn, v2_dyn, v3_dyn], [p12_dyn, p13_dyn, p23_dyn])Network with 7 states and 14 parameters
├─ 3 vertices (2 unique types)
└─ 3 edges (1 unique type)
Edge-Aggregation using SequentialAggregator(+)Initializing the Dynamic Model with the Static Solution
Now comes the important part: we need to initialize the dynamic model with the results from the static model. To do so, we need to make use of two functions:
interface_values: Extracts all interface values (inputs and outputs) from a network stateinitialize_componentwise!: Initializes all components in a network one by one
First, we extract all interface values from our static solution:
interface_vals = interface_values(u_static)OrderedCollections.OrderedDict{NetworkDynamics.SymbolicStateIndex{Int64, Symbol}, Float64} with 18 entries:
VIndex(1, :q̃_nw) => -1.0
VIndex(1, :p) => 1.0
VIndex(2, :q̃_nw) => 0.6
VIndex(2, :p) => 0.946667
VIndex(3, :q̃_nw) => 0.4
VIndex(3, :p) => 0.953333
EIndex(1, :p_src) => 1.0
EIndex(1, :p_dst) => 0.946667
EIndex(1, :₋q̃) => -0.533333
EIndex(1, :q̃) => 0.533333
EIndex(2, :p_src) => 1.0
EIndex(2, :p_dst) => 0.953333
EIndex(2, :₋q̃) => -0.466667
EIndex(2, :q̃) => 0.466667
EIndex(3, :p_src) => 0.946667
EIndex(3, :p_dst) => 0.953333
EIndex(3, :₋q̃) => 0.0666667
EIndex(3, :q̃) => -0.0666667Next, we initialize the dynamic model using these interface values as defaults:
u0_dyn = initialize_componentwise!(nw_dyn, default_overrides=interface_vals, verbose=true)NWState{Vector{Float64}} of Network (3 vertices, 3 edges)
├─ VIndex(1, :p) => 1.0
├─ VIndex(1, :ξ) => 0.9999999999999997
├─ VIndex(2, :p) => 0.9466666666666667
├─ VIndex(3, :p) => 0.9533333333333334
├─ EIndex(1, :q̃) => 0.5333333333333334
├─ EIndex(2, :q̃) => 0.46666666666666634
└─ EIndex(3, :q̃) => -0.0666666666666671
p = NWParameter([1.0, 0.1, 1.0, 1.0, 0.1, -0.6, 0.1, -0.4, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1])
t = nothingInternally, this function uses initialize_component! on every single component. For each component, it overwrites the defaults to be consistent with the interface values of the static model. Therefore, we ensure that the dynamic model is initialized near the steady state of the static model.
We can inspect individual components if needed:
dump_initial_state(nw_dyn[VIndex(1)])Inputs:
q̃_nw = -1
States:
p = 1
ξ = 1 (guess 0)
Outputs:
p = 1
Parameters:
C = 0.1
K_i = 1
K_p = 1
p_set = 1 (guess 1)
Observed:
q̃_inj = 1
q̃_prosumer = 1
Δp = 0Let's verify that our initialization is correct by checking that the derivatives are close to zero:
du = ones(dim(nw_dyn))
nw_dyn(du, uflat(u0_dyn), pflat(u0_dyn), 0.0)
extrema(du .- zeros(dim(nw_dyn))) # very close to zero, confirming we have a steady state!(-7.771561172376096e-15, 5.551115123125783e-15)Simulating the Dynamic Model
With our properly initialized model, we can now simulate the system to observe its behavior. To make the simulation more interesting, we'll introduce a disturbance to see how the system responds from its steady state.
We'll use a callback to increase consumer demand at a specific time. For more information on callbacks, see the documentation on Callbacks.
affect = ComponentAffect([], [:q̃_prosumer]) do u, p, ctx
@info "Increase consumer demand at t=$(ctx.t)"
p[:q̃_prosumer] -= 0.1
end
cb = PresetTimeComponentCallback([1.0], affect)
set_callback!(nw_dyn[VIndex(2)], cb) # attach disturbance to second nodeCreate and solve the ODE problem with the callback. Note that we're using the flat representation of our initialized state (via uflat and pflat) as input to the ODE solver:
prob = ODEProblem(nw_dyn, copy(uflat(u0_dyn)), (0, 7), copy(pflat(u0_dyn));
callback=get_callbacks(nw_dyn))
sol = solve(prob, Tsit5())[ Info: Increase consumer demand at t=1.0Visualizing the Results
Finally, let's visualize the results of our simulation. The plots show how our gas network responds to the increased consumer demand at t=1:
Pressure at nodes: We see a pressure drop at all nodes after the disturbance before the pressure is stabilized by the controller.
Injection by producer: Node 1 increases its injection to compensate for the higher demand.
Consumption by consumers: The solid lines show the actual flows at nodes 2 and 3, while the dashed lines show the set consumer demands. At t=1, we see the step change in consumer demand at node 2.
Flows through pipes: Shows how the flows in all pipes adjust to the new demand pattern.
let
fig = Figure(size=(1000,1000))
ax = Axis(fig[1, 1]; title="Pressure at nodes")
for i in 1:3
lines!(ax, sol; idxs=VIndex(i, :p), label="Node $i", color=Cycled(i))
end
ax = Axis(fig[2, 1]; title="Injection by producer")
lines!(ax, sol; idxs=VIndex(1, :q̃_inj), label="Node 1", color=Cycled(1))
ax = Axis(fig[3, 1]; title="Consumption by consumers")
for i in 2:3
lines!(ax, sol; idxs=@obsex(-1*VIndex(i, :q̃_inj)), label="Node $i", color=Cycled(i))
lines!(ax, sol; idxs=@obsex(-1*VIndex(i, :q̃_prosumer)), label="Node $i", linestyle=:dash, color=Cycled(i))
end
ax = Axis(fig[4, 1]; title="Flows through pipes")
for i in 1:3
lines!(ax, sol; idxs=@obsex(abs(EIndex(i, :q̃))), label="Pipe $i", color=Cycled(i))
end
fig
end
This page was generated using Literate.jl.