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 adds dynamics 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 equals the dynamic model in 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 interface values (pressures and flows) of the dynamic model with the results from the static model.
We can do this manually:
# Vertex 1: output
set_default!(nw_dyn[VIndex(1)], :p, u_static.v[1, :p])
# Vertex 1: input
set_default!(nw_dyn[VIndex(1)], :q̃_nw, u_static.v[1, :q̃_nw])But there is also a built-in method set_interface_defaults! which we can use automatically:
set_interface_defaults!(nw_dyn, u_static; verbose=true)Setting the interface defaults:
Vertex 1 output: p => 1.0
Vertex 1 input: q̃_nw => -0.9999999999999998
Vertex 2 output: p => 0.9466666666666667
Vertex 2 input: q̃_nw => 0.6000000000000005
Vertex 3 output: p => 0.9533333333333334
Vertex 3 input: q̃_nw => 0.39999999999999925
Edge 1 output: ₋q̃ => -0.5333333333333334
Edge 1 output: q̃ => 0.5333333333333334
Edge 1 input: p_src => 1.0
Edge 1 input: p_dst => 0.9466666666666667
Edge 2 output: ₋q̃ => -0.46666666666666634
Edge 2 output: q̃ => 0.46666666666666634
Edge 2 input: p_src => 1.0
Edge 2 input: p_dst => 0.9533333333333334
Edge 3 output: ₋q̃ => 0.0666666666666671
Edge 3 output: q̃ => -0.0666666666666671
Edge 3 input: p_src => 0.9466666666666667
Edge 3 input: p_dst => 0.9533333333333334With the interfaces all set, we can "initialize" the internal states of the dynamic models.
For example, let's inspect the state of our first vertex:
nw_dyn[VIndex(1)]VertexModel :v1_mod_dyn PureStateMap() @ Vertex 1
├─ 1 input: [q̃_nw=-1]
├─ 2 states: [ξ≈0, p=1]
├─ 1 output: [p=1]
└─ 4 params: [K_i=1, C=0.1, K_p=1, p_set≈1]We observe that both the initial state ξ as well as the pressure setpoint p_set are left "free". Using initialize_component!, we can try to find values for the "free" states and parameters such that the interface constraints are fulfilled.
initialize_component!(nw_dyn[VIndex(1)])VertexModel :v1_mod_dyn PureStateMap() @ Vertex 1
├─ 1 input: [q̃_nw=-1]
├─ 2 states: [ξ=1, p=1]
├─ 1 output: [p=1]
└─ 4 params: [K_i=1, C=0.1, K_p=1, p_set=1]We may also use dump_initial_state to get a more detailed view of the state:
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 = 0We can also initialize the other two vertices, however it is unnecessary since their state is already completely determined by the fixed input/output:
initialize_component!(nw_dyn[VIndex(2)])
initialize_component!(nw_dyn[VIndex(3)])[ Info: No free variables! Residual 5.551115123125783e-15
[ Info: No free variables! Residual 7.771561172376096e-15Similarly, we can initialize the dynamic pipe models. However, since their dynamic state equals the output, once again there is nothing to initialize.
initialize_component!(nw_dyn[EIndex(1)])
initialize_component!(nw_dyn[EIndex(2)])
initialize_component!(nw_dyn[EIndex(3)])[ Info: No free variables! Residual 0.0
[ Info: No free variables! Residual 0.0
[ Info: No free variables! Residual 0.0Now, everything is initialized, which means every input, output, state and parameter either has a default metadata or an init metadata. When constructing the NWState for this network, it will be filled with all those values which should now correspond to a steady state of the system:
u0_dyn = NWState(nw_dyn)NWState{Vector{Float64}} of Network (3 vertices, 3 edges)
├─ VIndex(1, :ξ) => 0.9999999999999997
├─ VIndex(1, :p) => 1.0
├─ 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 = nothingLet'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
Now we can solve the dynamic model and add a disturbance to see how the system responds:
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:
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.
Draw 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="Draw 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.