Merge pull request #195 from JuliaDynamics/hw/callbacks

Component Callbacks

Author: hexaeder

Project.toml

@@ -6,6 +6,7 @@ version = "0.9.7"
 [deps]
 ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197"
 Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458"
+DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -48,6 +49,7 @@ Adapt = "4.0.4"
 ArgCheck = "2.3.0"
 Atomix = "1"
 CUDA = "5.5.2"
+DiffEqCallbacks = "4.2.2"
 DocStringExtensions = "0.9.3"
 FastClosures = "0.3.2"
 ForwardDiff = "0.10.36"

docs/Project.toml

@@ -5,6 +5,7 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
 DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"

docs/examples/cascading_failure.jl

@@ -8,6 +8,12 @@ This script can be dowloaded as a normal Julia script [here](@__NAME__.jl). #md
 
 > [1] Schäfer, B., Witthaut, D., Timme, M., & Latora, V. (2018). Dynamically induced cascading failures in power grids. Nature communications, 9(1), 1-13. https://www.nature.com/articles/s41467-018-04287-5
 
+This example has three subchaperts:
+ - first we [define the network model](#define-the-model),
+ - secondly, we implement [component based callbacks](#component-based-callbacks) and
+ - thirdly we solve the problem using [systemwide callbacks](#system-wide-callbacks).
+
+
 The system is modeled using swing equation and active power edges. The nodes are
 characterized by the voltage angle `δ`, the active power on each line is symmetric
 and a function of the difference between source and destination angle `δ_src - δ_dst`.
@@ -22,6 +28,8 @@ using Test #hide
 import SymbolicIndexingInterface as SII
 
 #=
+## Defining the Model
+
 For the nodes we define the swing equation. State `v[1] = δ`, `v[2] = ω`.
 The swing equation has three parameters: `p = (P_ref, I, γ)` where `P_ref`
 is the power setpopint, `I` is the inertia and `γ` is the droop or damping coeficcient.
@@ -58,25 +66,102 @@ g = SimpleGraph([0 1 1 0 1;
                  1 1 0 1 0;
                  0 1 1 0 1;
                  1 0 0 1 0])
-swing_network = Network(g, vertex, edge)
+nw = Network(g, vertex, edge; dealias=true)
 
 #=
-For the parameters, we create the `NWParameter` object prefilled with default p values
+Note that we used `dealias=true` to automaticially generate separate
+`ComponentModels` for each vertex/edge. Doing so allows us to later
+set different metadata (callbacks, default values, etc.) for each vertex/edge.
+
+We proceed by setting the default reference power for the nodes:
 =#
-p = NWParameter(swing_network)
-## vertices 1, 3 and 4 act as loads
-p.v[(1,3,4), :P_ref] .= -1
-## vertices 2 and 5 act as generators
-p.v[(2,5), :P_ref] .= 1.5
+set_default!(nw, VIndex(1, :P_ref), -1.0) # load
+set_default!(nw, VIndex(2, :P_ref),  1.5) # generator
+set_default!(nw, VIndex(3, :P_ref), -1.0) # load
+set_default!(nw, VIndex(4, :P_ref), -1.0) # load
+set_default!(nw, VIndex(5, :P_ref),  1.5) # generator
 nothing #hide #md
 
 #=
-We can use `find_fixpoint` to find a valid initial condition of the network
+We can use `find_fixpoint` to find a valid initial condition of the network.
+We also use `set_defaults!` to overwirte all the default values for states and parameters
+with the one of the fixpoint, this means that we can allways re-extract this setpoint by
+using `u0 = NWState(nw)`.
+=#
+u0 = find_fixpoint(nw)
+set_defaults!(nw, u0)
+nothing #hide #md
+
+#=
+## Component-based Callbacks
+
+For the component based callback we need to define a condtion and an affect.
+Both functions take three inputs:
+  - the actual function `f`
+  - the states which to be accessed `sym`
+  - the parameters to be accessed `psym`
+=#
+cond = ComponentCondition([:P], [:limit]) do u, p, t
+    abs(u[:P]) - p[:limit]
+end
+affect = ComponentAffect([], [:K]) do u, p, ctx
+    println("Line $(ctx.eidx) tripped at t=$(ctx.integrator.t)")
+    p[:K] = 0
+end
+edge_cb = ContinousComponentCallback(cond, affect)
+#=
+To enable the callback in simulation, we need to attach them to the individual
+edgemodels/vertexmodels.
 =#
-u0 = find_fixpoint(swing_network, p)
+for i in 1:ne(g)
+    edgemodel = nw[EIndex(i)]
+    set_callback!(edgemodel, edge_cb)
+end
 nothing #hide #md
 
 #=
+The system starts at a steady state.
+In order to see any dynamic, we need to fail a first line intentionally.
+For that we use a [`PresetTimeComponentCallback`](@ref), which triggers an
+`ComponentAffect` at a given time. We can reuse the previously defined component
+affect for that and just add it to line number 5 at time 1.0.
+=#
+trip_first_cb = PresetTimeComponentCallback(1.0, affect)
+add_callback!(nw[EIndex(5)], trip_first_cb)
+nothing #hide #md
+#=
+When we inspect the edge model for 5 no, we see that we've registered 2 callbacks:
+=#
+nw[EIndex(5)]
+
+#=
+Now we can simulate the network. We use [`get_callbacks(::Network)`](@ref)
+to generate a callback set for the whole network which represents all of the individual
+component callbacks.
+=#
+u0 = NWState(nw)
+network_cb = get_callbacks(nw)
+prob = ODEProblem(nw, uflat(u0), (0, 6), pflat(u0); callback=network_cb)
+sol = solve(prob, Tsit5());
+nothing #hide #md
+
+#=
+Lastly we plot the power flow on all lines using the [`eidxs`](@ref) function to generate the
+symbolic indices for the states of interest:
+=#
+plot(sol; idxs=eidxs(sol, :, :P))
+
+#=
+## System wide Callbacks
+
+The above solution relies on the `ComponentCallback` features of
+NetworkDyanmics. The "low-level" API would be to use `VectorContinousCallback`
+and `PresetTimeCallback` directly to achieve the same effect, essentially doing
+manually what [`get_callbacks(::Network)`](@ref) is doing for us.
+
+While not necessary in this case, this method offers more flexiblity then the
+component based appraoch.
+
 In order to implement the line failures, we need to create a `VectorContinousCallback`.
 In the callback, we compare the current flow on the line with the limit. If the limit is reached,
 the coupling `K` is set to 0.
@@ -125,8 +210,8 @@ cache accessors to the internal states.
 This still isn't ideal beacuse both `getlim` and `getflow` getters will create arrays
 within the callback. But is far better then resolving the flat state indices every time.
 =#
-condition = let getlim = SII.getp(swing_network, epidxs(swing_network, :, :limit)),
-                getflow = SII.getu(swing_network, eidxs(swing_network, :, :P))
+condition = let getlim = SII.getp(nw, epidxs(nw, :, :limit)),
+                getflow = SII.getu(nw, eidxs(nw, :, :P))
     function (out, u, t, integrator)
         ## careful,  u != integrator.u
         ## therefore construct nwstate with Network info from integrator but u
@@ -143,28 +228,26 @@ We can combine affect and condition to form the callback.
 trip_cb = VectorContinuousCallback(condition, affect!, ne(g));
 
 #=
-However, there is another component missing. If we look at the powerflow on the
-lines in the initial steady state
-=#
-u0.e[:, :P]
-#=
+Similarily to before, we need to generate a initial perturbation by failing one line
+using a `PresetTimeCallback`.
 We see that every flow is below the trip value 1.0. Therefor we need to add a distrubance
 to the network. We do this by manually disabeling line 5 at time 1.
 =#
 trip_first_cb = PresetTimeCallback(1.0, integrator->affect!(integrator, 5));
 
 #=
-With those components, we can create the problem and solve it.
+Now we are set for solving the system again. This time we create our own callback
+set by combining both Callbacks manually.
 =#
-
-prob = ODEProblem(swing_network, uflat(u0), (0,6), copy(pflat(p));
-                  callback=CallbackSet(trip_cb, trip_first_cb))
+u0 = NWState(nw)
+cbset = CallbackSet(trip_cb, trip_first_cb)
+prob = ODEProblem(nw, uflat(u0), (0,6), pflat(u0); callback=cbset)
 Main.test_execution_styles(prob) # testing all ex styles #src
-sol = solve(prob, Tsit5());
+sol2 = solve(prob, Tsit5());
 ## we want to test the reconstruction of the observables # hide
-@test all(!iszero, sol(sol.t; idxs=eidxs(sol,:,:P))[begin]) # hide
-@test all(iszero, sol(sol.t; idxs=eidxs(sol,:,:P))[end][[1:5...,7]]) # hide
+@test all(!iszero, sol2(sol2.t; idxs=eidxs(sol2,:,:P))[begin]) # hide
+@test all(iszero, sol2(sol2.t; idxs=eidxs(sol2,:,:P))[end][[1:5...,7]]) # hide
 nothing #hide
 
-# Through the magic of symbolic indexing we can plot the power flows on all lines:
-plot(sol; idxs=eidxs(sol,:,:P))
+# Then again we plot the solution:
+plot(sol2; idxs=eidxs(sol2,:,:P))

docs/make.jl

@@ -5,6 +5,14 @@ using NetworkDynamics
 using SciMLBase
 using Literate
 using ModelingToolkit
+using DocumenterInterLinks
+
+links = InterLinks(
+    "DiffEq" => "https://docs.sciml.ai/DiffEqDocs/stable/",
+    "MTK" => "https://docs.sciml.ai/ModelingToolkit/stable/",
+    "SymbolicIndexingInterface" => "https://docs.sciml.ai/SymbolicIndexingInterface/stable/",
+    "DiffEqCallbacks" => "https://docs.sciml.ai/DiffEqCallbacks/stable/",
+)
 
 # generate examples
 example_dir = joinpath(@__DIR__, "examples")
@@ -17,16 +25,14 @@ for example in filter(contains(r".jl$"), readdir(example_dir, join=true))
     Literate.script(example, outdir; keep_comments=true)
 end
 
-# TODO: doc on steady state solve https://docs.sciml.ai/NonlinearSolve/stable/native/steadystatediffeq/#SteadyStateDiffEq.SSRootfind
-# TODO: doc on parameter & state handling? -> symbolic indexing
-
 mtkext = Base.get_extension(NetworkDynamics, :MTKExt)
-makedocs(;
+kwargs = (;
     root=joinpath(pkgdir(NetworkDynamics), "docs"),
     sitename="NetworkDynamics",
     modules=[NetworkDynamics, mtkext],
     linkcheck=true, # checks if external links resolve
     pagesonly=true,
+    plugins=[links],
     pages=[
         "General" => "index.md",
         "mathematical_model.md",
@@ -35,6 +41,7 @@ makedocs(;
             "symbolic_indexing.md",
             "metadata.md",
             "initialization.md",
+            "callbacks.md",
             "mtk_integration.md",
             "external_inputs.md",
         ],
@@ -52,12 +59,30 @@ makedocs(;
     ],
     draft=false,
     format = Documenter.HTML(ansicolor = true),
-    warnonly=[:cross_references, :missing_docs, :docs_block],
-    # warnonly=true,
+    warnonly=[:missing_docs],
 )
+kwargs_warnonly = (; kwargs..., warnonly=true)
+
+if haskey(ENV,"GITHUB_ACTIONS")
+    success = true
+    thrown_ex = nothing
+    try
+        makedocs(; kwargs...)
+    catch e
+        @info "Strict doc build failed, try again with warnonly=true"
+        global success = false
+        global thrown_ex = e
+        makedocs(; kwargs_warnonly...)
+    end
+
+    deploydocs(; repo="github.com/JuliaDynamics/NetworkDynamics.jl.git",
+            devbranch="main", push_preview=true)
+
+    success || throw(thrown_ex)
+else # local build
+    makedocs(; kwargs_warnonly...)
+end
 
-deploydocs(; repo="github.com/JuliaDynamics/NetworkDynamics.jl.git",
-           devbranch="main", push_preview=true)
 
 # warnonly options
 # :autodocs_block

docs/src/API.md

@@ -124,13 +124,35 @@ get_bounds
 set_bounds!
 ```
 
+## Callbacks API
+### Define Callbacks
+```@docs
+NetworkDynamics.ComponentCallback
+ContinousComponentCallback
+VectorContinousComponentCallback
+DiscreteComponentCallback
+PresetTimeComponentCallback
+ComponentCondition
+ComponentAffect
+SymbolicView
+get_callbacks(::NetworkDynamics.Network)
+```
+### Attach Callbacks to Edge/VertexModels
+```@docs
+has_callback
+get_callbacks(::NetworkDynamics.ComponentModel)
+set_callback!
+add_callback!
+```
+
 ## Initialization
 ```@docs
 find_fixpoint
 initialize_component!
 init_residual
 get_initial_state
 dump_initial_state
+set_defaults!
 ```
 
 ## Execution Types