add set_interface_defaults! method

Author: hexaeder

docs/examples/init_tutorial.jl

@@ -44,7 +44,7 @@ nothing #hide
 #=
 Now we'll define three specific node types:
 
-1. A constant pressure node that forces pressure to maintain a specific value
+**A) A constant pressure node that forces pressure to maintain a specific value**
 =#
 @mtkmodel ConstantPressureNode begin
     @extend GasNode()
@@ -58,7 +58,7 @@ end
 nothing #hide
 
 #=
-2. A static prosumer node which forces a certain flow (pressure is fully implicit)
+**B) A static prosumer node which forces a certain flow (pressure is fully implicit)**
 =#
 @mtkmodel StaticProsumerNode begin
     @extend GasNode()
@@ -72,7 +72,7 @@ end
 nothing #hide
 
 #=
-3. A dynamic prosumer node with compliance, which adds dynamics to the pressure state
+**C) A dynamic prosumer node with compliance, which adds dynamics to the pressure state**
 =#
 @mtkmodel DynamicProsumerNode begin
     @extend GasNode()
@@ -87,7 +87,7 @@ end
 nothing #hide
 
 #=
-4. A pressure control node that tries to maintain a set pressure by adjusting its injection
+**D) A pressure control node that tries to maintain a set pressure by adjusting its injection**
 =#
 @mtkmodel PressureControlNode begin
     @extend GasNode()
@@ -206,7 +206,7 @@ u_static = find_fixpoint(nw_static, u_static_guess)
 
 Now we'll define our dynamic model using more complex components:
 =#
-@named v1_mod_dyn = PressureControlNode(;p_set=1)
+@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)
@@ -236,37 +236,62 @@ nw_dyn = Network([v1_dyn, v2_dyn, v3_dyn], [p12_dyn, p13_dyn, p23_dyn])
 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.
 
-First, let's handle node 1 (the pressure control node):
+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])
-v1_dyn # hide
+nothing #hide
+
+#=
+But there is also a built in method [`set_interface_defaults!`](@ref) which we can use
+automaticially:
+=#
+set_interface_defaults!(nw_dyn, u_static; verbose=true)
+nothing #hide
+
 #=
-In the output, you can see that state ξ is "approx" 1 (guess value) while the rest is fixed.
-Now we can initialize the dynamic model's internal states:
+With 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)]
+
+#=
+We observe, that both the initial state `ξ` as well as the pressure setpoint `p_set`
+is left "free". Using [`initialize_component!`](@ref), we can try to find values for the
+"free" states and parameters such that the interface constraints are fulfilled.
 =#
-initialize_component!(v1_dyn)
-v1_dyn #hide
+initialize_component!(nw_dyn[VIndex(1)])
+#=
+We may also use [`dump_initial_state`](@ref) to get a more detailed view of the state:
+=#
+dump_initial_state(nw_dyn[VIndex(1)])
+nothing #hide
 
 #=
-For the other two vertices (which are simpler), we just need to set the default values:
+We can also initialize the other two vertices, however it is a bit useless
+since their state is allready completely determined by the fixed input/output:
 =#
-set_default!(nw_dyn[VIndex(2)], :p, u_static.v[2, :p])
-set_default!(nw_dyn[VIndex(3)], :p, u_static.v[3, :p])
+initialize_component!(nw_dyn[VIndex(2)])
+initialize_component!(nw_dyn[VIndex(3)])
 nothing #hide
 
 #=
-For the pipe models, we manually initialize the flow state of the dynamic line model:
+Similarily, we can initialize the dynamic pipe models. However since their dynamic state
+equals the output, once again there is nothing to initialize.
 =#
-set_default!(nw_dyn[EIndex(1)], :q̃, u_static.e[1, :q̃])
-set_default!(nw_dyn[EIndex(2)], :q̃, u_static.e[2, :q̃])
-set_default!(nw_dyn[EIndex(3)], :q̃, u_static.e[3, :q̃])
+initialize_component!(nw_dyn[EIndex(1)])
+initialize_component!(nw_dyn[EIndex(2)])
+initialize_component!(nw_dyn[EIndex(3)])
 nothing #hide
 
 #=
-Now that we've set all the "default" values for all the states, we can call `NWState` on the
-network to get a fully initialized state vector:
+Now, everything is initialized, which means every input, output, state and parameter
+either has a `default` metadata or a `init` metadate. 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)
 

docs/src/API.md

@@ -164,6 +164,7 @@ init_residual
 get_initial_state
 dump_initial_state
 set_defaults!
+set_interface_defaults!
 ```
 
 ## Execution Types

src/NetworkDynamics.jl

@@ -87,7 +87,7 @@ include("callbacks.jl")
 using NonlinearSolve: AbstractNonlinearSolveAlgorithm, NonlinearFunction
 using NonlinearSolve: NonlinearLeastSquaresProblem, NonlinearProblem
 using SteadyStateDiffEq: SteadyStateProblem, SteadyStateDiffEqAlgorithm, SSRootfind
-export find_fixpoint, initialize_component!, init_residual
+export find_fixpoint, initialize_component!, init_residual, set_interface_defaults!
 include("initialization.jl")
 
 include("show.jl")

src/initialization.jl

@@ -66,6 +66,10 @@ function initialization_problem(cf::T; t=NaN, apply_bound_transformation, verbos
 
     # count free variables and equations
     Nfree = mapreduce(sum, +, outfree_ms) + sum(ufree_m) + mapreduce(sum, +, infree_ms) + sum(pfree_m)
+
+    # return "fake NonlinearProblem" if no free variables
+    iszero(Nfree) && return ((; u0=[]), identity)
+
     Neqs = dim(cf)  + mapreduce(length, +, outfree_ms)
 
 
@@ -269,7 +273,11 @@ function initialize_component!(cf; verbose=true, apply_bound_transformation=true
         resid = sol.resid
     else
         resid = init_residual(cf; recalc=true)
+        if resid < 1e-10
         verbose && @info "No free variables! Residual $(LinearAlgebra.norm(resid))"
+        else
+            @warn "No free variables! However model does not appear to be initialized in steady state. Residual $(LinearAlgebra.norm(resid))"
+        end
     end
 
     set_metadata!(cf, :init_residual, resid)
@@ -321,7 +329,7 @@ function init_residual(cf::T; t=NaN, recalc=false) where {T<:ComponentModel}
 end
 
 function broken_bounds(cf)
-    allsyms = vcat(sym(cf), psym(cf), insym_all(cf), outsym_flat(cf), obssym(cf))
+    allsyms = vcat(sym(cf), psym(cf), insym_flat(cf), outsym_flat(cf), obssym(cf))
     boundsyms = filter(s -> has_bounds(cf, s), allsyms)
     bounds = get_bounds.(Ref(cf), boundsyms)
     vals = get_initial_state(cf, boundsyms)
@@ -349,3 +357,50 @@ function set_defaults!(nw::Network, s::NWState)
     end
     nw
 end
+
+"""
+    set_interface_defaults!(nw::Network, s::NWState; verbose=false)
+
+Sets the **interface** (i.e. node and edge inputs/outputs) defaults of a given
+network to the ones defined by the given state. Notably, while the graph
+topology and interface dimensions of the target network `nw` and the source
+network of `s` musst be identicaly, the systems may differ in the dynamical
+components.
+
+This is mainly inteded for initialization purposes: solve the interface values
+with a simpler -- possible static -- network and "transfer" the steady state
+interface values to the full network.
+"""
+function set_interface_defaults!(nw::Network, s::NWState; verbose=false)
+    @argcheck s.nw.im.g == nw.im.g "Graphs musst have the same structure!"
+    verbose && println("Setting the interface defaults:")
+    values = Dict{SymbolicIndex, Float64}()
+    for (idxs, _Index, comp) in ((1:nv(nw), VIndex, "Vertex"), (1:ne(nw), EIndex, "Edge"))
+        for i in idxs
+            target = nw[_Index(i)]
+            source = s.nw[_Index(i)]
+
+            tsyms = _Index(i, outsym_flat(target)) |> collect
+            svals = s[_Index(i, outsym_flat(source))]
+            @assert length(tsyms) == length(svals) "$comp $i: Source network has $(length(svals)) outputs, but target network has $(length(tsyms)) outputs!"
+            for (tsym, sval) in zip(tsyms, svals)
+                values[tsym] = sval
+                verbose && println("  $comp $(lpad(i, length(repr(idxs[end])))) \
+                                    output: $(tsym.subidx) => $sval")
+            end
+
+            tsyms = _Index(i, insym_flat(target)) |> collect
+            svals = s[_Index(i, insym_flat(source))]
+            @assert length(tsyms) == length(svals) "$comp $i: Source network has $(length(svals)) inputs, but target network has $(length(tsyms)) inputs!"
+            for (tsym, sval) in zip(tsyms, svals)
+                values[tsym] = sval
+                verbose && println("  $comp $(lpad(i, length(repr(idxs[end])))) \
+                                    input:  $(tsym.subidx) => $sval")
+            end
+        end
+    end
+    for (sym, val) in values
+        set_default!(nw, sym, val)
+    end
+    nw
+end