diff --git a/docs/Project.toml b/docs/Project.toml index 8d0f1667c9..24f00b4db2 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -15,6 +15,7 @@ OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226" +Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index f425fdce5b..23e1e6d7d1 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -378,3 +378,218 @@ sol.ps[c] # sol[c] will error, since `c` is not a timeseries value ``` It can be seen that the timeseries for `c` is not saved. + +## [(Experimental) Imperative affects](@id imp_affects) + +The `ImperativeAffect` can be used as an alternative to the aforementioned functional affect form. Note +that `ImperativeAffect` is still experimental; to emphasize this, we do not export it and it should be +included as `ModelingToolkit.ImperativeAffect`. `ImperativeAffect` aims to simplify the manipulation of +system state. + +We will use two examples to describe `ImperativeAffect`: a simple heater and a quadrature encoder. +These examples will also demonstrate advanced usage of `ModelingToolkit.SymbolicContinuousCallback`, +the low-level interface of the tuple form converts into that allows control over the SciMLBase-level +event that is generated for a continuous event. + +### [Heater](@id heater_events) + +Bang-bang control of a heater connected to a leaky plant requires hysteresis in order to prevent rapid control oscillation. + +```@example events +@variables temp(t) +params = @parameters furnace_on_threshold=0.5 furnace_off_threshold=0.7 furnace_power=1.0 leakage=0.1 furnace_on(t)::Bool=false +eqs = [ + D(temp) ~ furnace_on * furnace_power - temp^2 * leakage +] +``` + +Our plant is simple. We have a heater that's turned on and off by the time-indexed parameter `furnace_on` +which adds `furnace_power` forcing to the system when enabled. We then leak heat proportional to `leakage` +as a function of the square of the current temperature. + +We need a controller with hysteresis to control the plant. We wish the furnace to turn on when the temperature +is below `furnace_on_threshold` and off when above `furnace_off_threshold`, while maintaining its current state +in between. To do this, we create two continuous callbacks: + +```@example events +using Setfield +furnace_disable = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i + @set! x.furnace_on = false + end) +furnace_enable = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_on_threshold], + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i + @set! x.furnace_on = true + end) +``` + +We're using the explicit form of `SymbolicContinuousCallback` here, though +so far we aren't using anything that's not possible with the implicit interface. +You can also write + +```julia +[temp ~ furnace_off_threshold] => ModelingToolkit.ImperativeAffect(modified = (; + furnace_on)) do x, o, i, c + @set! x.furnace_on = false +end +``` + +and it would work the same. + +The `ImperativeAffect` is the larger change in this example. `ImperativeAffect` has the constructor signature + +```julia +ImperativeAffect(f::Function; modified::NamedTuple, observed::NamedTuple, ctx) +``` + +that accepts the function to call, a named tuple of both the names of and symbolic values representing +values in the system to be modified, a named tuple of the values that are merely observed (that is, used from +the system but not modified), and a context that's passed to the affect function. + +In our example, each event merely changes whether the furnace is on or off. Accordingly, we pass a `modified` tuple +`(; furnace_on)` (creating a `NamedTuple` equivalent to `(furnace_on = furnace_on)`). `ImperativeAffect` will then +evaluate this before calling our function to fill out all of the numerical values, then apply them back to the system +once our affect function returns. Furthermore, it will check that it is possible to do this assignment. + +The function given to `ImperativeAffect` needs to have the signature: + +```julia +f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple +``` + +The function `f` will be called with `observed` and `modified` `NamedTuple`s that are derived from their respective `NamedTuple` definitions. +In our example, if `furnace_on` is `false`, then the value of the `x` that's passed in as `modified` will be `(furnace_on = false)`. +The modified values should be passed out in the same format: to set `furnace_on` to `true` we need to return a tuple `(furnace_on = true)`. +The examples does this with Setfield, recreating the result tuple before returning it; the returned tuple may optionally be missing values as +well, in which case those values will not be written back to the problem. + +Accordingly, we can now interpret the `ImperativeAffect` definitions to mean that when `temp = furnace_off_threshold` we +will write `furnace_on = false` back to the system, and when `temp = furnace_on_threshold` we will write `furnace_on = true` back +to the system. + +```@example events +@named sys = ODESystem( + eqs, t, [temp], params; continuous_events = [furnace_disable, furnace_enable]) +ss = structural_simplify(sys) +prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 10.0)) +sol = solve(prob, Tsit5()) +plot(sol) +hline!([sol.ps[furnace_off_threshold], sol.ps[furnace_on_threshold]], + l = (:black, 1), primary = false) +``` + +Here we see exactly the desired hysteresis. The heater starts on until the temperature hits +`furnace_off_threshold`. The temperature then bleeds away until `furnace_on_threshold` at which +point the furnace turns on again until `furnace_off_threshold` and so on and so forth. The controller +is effectively regulating the temperature of the plant. + +### [Quadrature Encoder](@id quadrature) + +For a more complex application we'll look at modeling a quadrature encoder attached to a shaft spinning at a constant speed. +Traditionally, a quadrature encoder is built out of a code wheel that interrupts the sensors at constant intervals and two sensors slightly out of phase with one another. +A state machine can take the pattern of pulses produced by the two sensors and determine the number of steps that the shaft has spun. The state machine takes the new value +from each sensor and the old values and decodes them into the direction that the wheel has spun in this step. + +```@example events +@variables theta(t) omega(t) +params = @parameters qA=0 qB=0 hA=0 hB=0 cnt::Int=0 +eqs = [D(theta) ~ omega + omega ~ 1.0] +``` + +Our continuous-time system is extremely simple. We have two unknown variables `theta` for the angle of the shaft +and `omega` for the rate at which it's spinning. We then have parameters for the state machine `qA, qB, hA, hB` +(corresponding to the current quadrature of the A/B sensors and the historical ones) and a step count `cnt`. + +We'll then implement the decoder as a simple Julia function. + +```@example events +function decoder(oldA, oldB, newA, newB) + state = (oldA, oldB, newA, newB) + if state == (0, 0, 1, 0) || state == (1, 0, 1, 1) || state == (1, 1, 0, 1) || + state == (0, 1, 0, 0) + return 1 + elseif state == (0, 0, 0, 1) || state == (0, 1, 1, 1) || state == (1, 1, 1, 0) || + state == (1, 0, 0, 0) + return -1 + elseif state == (0, 0, 0, 0) || state == (0, 1, 0, 1) || state == (1, 0, 1, 0) || + state == (1, 1, 1, 1) + return 0 + else + return 0 # err is interpreted as no movement + end +end +``` + +Based on the current and old state, this function will return 1 if the wheel spun in the positive direction, +-1 if in the negative, and 0 otherwise. + +The encoder state advances when the occlusion begins or ends. We model the +code wheel as simply detecting when `cos(100*theta)` is 0; if we're at a positive +edge of the 0 crossing, then we interpret that as occlusion (so the discrete `qA` goes to 1). Otherwise, if `cos` is +going negative, we interpret that as lack of occlusion (so the discrete goes to 0). The decoder function is +then invoked to update the count with this new information. + +We can implement this in one of two ways: using edge sign detection or right root finding. For exposition, we +will implement each sensor differently. + +For sensor A, we're using the edge detection method. By providing a different affect to `SymbolicContinuousCallback`'s +`affect_neg` argument, we can specify different behaviour for the negative crossing vs. the positive crossing of the root. +In our encoder, we interpret this as occlusion or nonocclusion of the sensor, update the internal state, and tick the decoder. + +```@example events +qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0], + ModelingToolkit.ImperativeAffect((; qA, hA, hB, cnt), (; qB)) do x, o, c, i + @set! x.hA = x.qA + @set! x.hB = o.qB + @set! x.qA = 1 + @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) + x + end, + affect_neg = ModelingToolkit.ImperativeAffect( + (; qA, hA, hB, cnt), (; qB)) do x, o, c, i + @set! x.hA = x.qA + @set! x.hB = o.qB + @set! x.qA = 0 + @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) + x + end) +``` + +The other way we can implement a sensor is by changing the root find. +Normally, we use left root finding; the affect will be invoked instantaneously _before_ +the root is crossed. This makes it trickier to figure out what the new state is. +Instead, we can use right root finding: + +```@example events +qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π / 2) ~ 0], + ModelingToolkit.ImperativeAffect((; qB, hA, hB, cnt), (; qA, theta)) do x, o, c, i + @set! x.hA = o.qA + @set! x.hB = x.qB + @set! x.qB = clamp(sign(cos(100 * o.theta - π / 2)), 0.0, 1.0) + @set! x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) + x + end; rootfind = SciMLBase.RightRootFind) +``` + +Here, sensor B is located `π / 2` behind sensor A in angular space, so we're adjusting our +trigger function accordingly. We here ask for right root finding on the callback, so we know +that the value of said function will have the "new" sign rather than the old one. Thus, we can +determine the new state of the sensor from the sign of the indicator function evaluated at the +affect activation point, with -1 mapped to 0. + +We can now simulate the encoder. + +```@example events +@named sys = ODESystem( + eqs, t, [theta, omega], params; continuous_events = [qAevt, qBevt]) +ss = structural_simplify(sys) +prob = ODEProblem(ss, [theta => 0.0], (0.0, pi)) +sol = solve(prob, Tsit5(); dtmax = 0.01) +sol.ps[cnt] +``` + +`cos(100*theta)` will have 200 crossings in the half rotation we've gone through, so the encoder would notionally count 200 steps. +Our encoder counts 198 steps (it loses one step to initialization and one step due to the final state falling squarely on an edge). diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index de9a06d774..20b2ada8fa 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -146,6 +146,7 @@ include("systems/parameter_buffer.jl") include("systems/abstractsystem.jl") include("systems/model_parsing.jl") include("systems/connectors.jl") +include("systems/imperative_affect.jl") include("systems/callbacks.jl") include("systems/problem_utils.jl") diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index e5db9a71b6..eaf31a9c5f 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -62,8 +62,6 @@ function Base.hash(a::FunctionalAffect, s::UInt) hash(a.ctx, s) end -has_functional_affect(cb) = affects(cb) isa FunctionalAffect - namespace_affect(affect, s) = namespace_equation(affect, s) function namespace_affect(affect::FunctionalAffect, s) FunctionalAffect(func(affect), @@ -75,6 +73,10 @@ function namespace_affect(affect::FunctionalAffect, s) context(affect)) end +function has_functional_affect(cb) + (affects(cb) isa FunctionalAffect || affects(cb) isa ImperativeAffect) +end + #################################### continuous events ##################################### const NULL_AFFECT = Equation[] @@ -88,17 +90,26 @@ By default `affect_neg = affect`; to only get rising edges specify `affect_neg = Assume without loss of generality that the equation is of the form `c(u,p,t) ~ 0`; we denote the integrator state as `i.u`. For compactness, we define `prev_sign = sign(c(u[t-1], p[t-1], t-1))` and `cur_sign = sign(c(u[t], p[t], t))`. A condition edge will be detected and the callback will be invoked iff `prev_sign * cur_sign <= 0`. +The positive edge `affect` will be triggered iff an edge is detected and if `prev_sign < 0`; similarly, `affect_neg` will be +triggered iff an edge is detected and `prev_sign > 0`. + Inter-sample condition activation is not guaranteed; for example if we use the dirac delta function as `c` to insert a sharp discontinuity between integrator steps (which in this example would not normally be identified by adaptivity) then the condition is not guaranteed to be triggered. Once detected the integrator will "wind back" through a root-finding process to identify the point when the condition became active; the method used -is specified by `rootfind` from [`SciMLBase.RootfindOpt`](@ref). Multiple callbacks in the same system with different `rootfind` operations will be resolved -into separate VectorContinuousCallbacks in the enumeration order of `SciMLBase.RootfindOpt`, which may cause some callbacks to not fire if several become -active at the same instant. See the `SciMLBase` documentation for more information on the semantic rules. - -The positive edge `affect` will be triggered iff an edge is detected and if `prev_sign < 0`; similarly, `affect_neg` will be -triggered iff an edge is detected `prev_sign > 0`. +is specified by `rootfind` from [`SciMLBase.RootfindOpt`](@ref). If we denote the time when the condition becomes active as `tc`, +the value in the integrator after windback will be: +* `u[tc-epsilon], p[tc-epsilon], tc` if `LeftRootFind` is used, +* `u[tc+epsilon], p[tc+epsilon], tc` if `RightRootFind` is used, +* or `u[t], p[t], t` if `NoRootFind` is used. +For example, if we want to detect when an unknown variable `x` satisfies `x > 0` using the condition `x ~ 0` on a positive edge (that is, `D(x) > 0`), +then left root finding will get us `x=-epsilon`, right root finding `x=epsilon` and no root finding will produce whatever the next step of the integrator was after +it passed through 0. + +Multiple callbacks in the same system with different `rootfind` operations will be grouped +by their `rootfind` value into separate VectorContinuousCallbacks in the enumeration order of `SciMLBase.RootfindOpt`. This may cause some callbacks to not fire if several become +active at the same instant. See the `SciMLBase` documentation for more information on the semantic rules. Affects (i.e. `affect` and `affect_neg`) can be specified as either: * A list of equations that should be applied when the callback is triggered (e.g. `x ~ 3, y ~ 7`) which must be of the form `unknown ~ observed value` where each `unknown` appears only once. Equations will be applied in the order that they appear in the vector; parameters and state updates will become immediately visible to following equations. @@ -108,6 +119,7 @@ Affects (i.e. `affect` and `affect_neg`) can be specified as either: + `read_parameters` is a vector of the parameters that are *used* by `f!`. Their indices are passed to `f` in `p` similarly to the indices of `unknowns` passed in `u`. + `modified_parameters` is a vector of the parameters that are *modified* by `f!`. Note that a parameter will not appear in `p` if it only appears in `modified_parameters`; it must appear in both `parameters` and `modified_parameters` if it is used in the affect definition. + `ctx` is a user-defined context object passed to `f!` when invoked. This value is aliased for each problem. +* A [`ImperativeAffect`](@ref); refer to its documentation for details. DAEs will be reinitialized using `reinitializealg` (which defaults to `SciMLBase.CheckInit`) after callbacks are applied. This reinitialization algorithm ensures that the DAE is satisfied after the callback runs. The default value of `CheckInit` will simply validate @@ -119,10 +131,10 @@ will run as soon as the solver starts, while finalization affects will be execut """ struct SymbolicContinuousCallback eqs::Vector{Equation} - initialize::Union{Vector{Equation}, FunctionalAffect} - finalize::Union{Vector{Equation}, FunctionalAffect} - affect::Union{Vector{Equation}, FunctionalAffect} - affect_neg::Union{Vector{Equation}, FunctionalAffect, Nothing} + initialize::Union{Vector{Equation}, FunctionalAffect, ImperativeAffect} + finalize::Union{Vector{Equation}, FunctionalAffect, ImperativeAffect} + affect::Union{Vector{Equation}, FunctionalAffect, ImperativeAffect} + affect_neg::Union{Vector{Equation}, FunctionalAffect, ImperativeAffect, Nothing} rootfind::SciMLBase.RootfindOpt reinitializealg::SciMLBase.DAEInitializationAlgorithm function SymbolicContinuousCallback(; @@ -369,7 +381,7 @@ SymbolicDiscreteCallback(cb::SymbolicDiscreteCallback) = cb # passthrough function Base.show(io::IO, db::SymbolicDiscreteCallback) println(io, "condition: ", db.condition) println(io, "affects:") - if db.affects isa FunctionalAffect + if db.affects isa FunctionalAffect || db.affects isa ImperativeAffect # TODO println(io, " ", db.affects) else @@ -722,6 +734,7 @@ function generate_single_rootfinding_callback( else initfn = user_initfun end + return ContinuousCallback( cond, affect_function.affect, affect_function.affect_neg, rootfind = cb.rootfind, initialize = initfn, @@ -755,7 +768,6 @@ function generate_vector_rootfinding_callback( finalize::Union{Function, Nothing}}[ compile_affect_fn(cb, sys, dvs, ps, kwargs) for cb in cbs] - cond = function (out, u, t, integ) rf_ip(out, u, parameter_values(integ), t) end @@ -801,31 +813,21 @@ function generate_vector_rootfinding_callback( initialize = nothing if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing initialize = handle_optional_setup_fn( - map( - (cb, fn) -> begin - if (save_idxs = get(ic.callback_to_clocks, cb, nothing)) !== nothing - let save_idxs = save_idxs - if !isnothing(fn.initialize) - (i) -> begin - fn.initialize(i) - for idx in save_idxs - SciMLBase.save_discretes!(i, idx) - end - end - else - (i) -> begin - for idx in save_idxs - SciMLBase.save_discretes!(i, idx) - end - end + map(cbs, affect_functions) do cb, fn + if (save_idxs = get(ic.callback_to_clocks, cb, nothing)) !== nothing + let save_idxs = save_idxs + custom_init = fn.initialize + (i) -> begin + !isnothing(custom_init) && custom_init(i) + for idx in save_idxs + SciMLBase.save_discretes!(i, idx) end end - else - fn.initialize end - end, - cbs, - affect_functions), + else + fn.initialize + end + end, SciMLBase.INITIALIZE_DEFAULT) else @@ -847,6 +849,13 @@ function compile_affect_fn(cb, sys::AbstractTimeDependentSystem, dvs, ps, kwargs eq_aff = affects(cb) eq_neg_aff = affect_negs(cb) affect = compile_affect(eq_aff, cb, sys, dvs, ps; expression = Val{false}, kwargs...) + function compile_optional_affect(aff, default = nothing) + if isnothing(aff) || aff == default + return nothing + else + return compile_affect(aff, cb, sys, dvs, ps; expression = Val{false}, kwargs...) + end + end if eq_neg_aff === eq_aff affect_neg = affect else @@ -944,10 +953,48 @@ function compile_user_affect(affect::FunctionalAffect, cb, sys, dvs, ps; kwargs. end end +function invalid_variables(sys, expr) + filter(x -> !any(isequal(x), all_symbols(sys)), reduce(vcat, vars(expr); init = [])) +end +function unassignable_variables(sys, expr) + assignable_syms = reduce( + vcat, Symbolics.scalarize.(vcat(unknowns(sys), parameters(sys))); init = []) + written = reduce(vcat, Symbolics.scalarize.(vars(expr)); init = []) + return filter( + x -> !any(isequal(x), assignable_syms), written) +end + +@generated function _generated_writeback(integ, setters::NamedTuple{NS1, <:Tuple}, + values::NamedTuple{NS2, <:Tuple}) where {NS1, NS2} + setter_exprs = [] + for name in NS2 + if !(name in NS1) + missing_name = "Tried to write back to $name from affect; only declared states ($NS1) may be written to." + error(missing_name) + end + push!(setter_exprs, :(setters.$name(integ, values.$name))) + end + return :(begin + $(setter_exprs...) + end) +end + +function check_assignable(sys, sym) + if symbolic_type(sym) == ScalarSymbolic() + is_variable(sys, sym) || is_parameter(sys, sym) + elseif symbolic_type(sym) == ArraySymbolic() + is_variable(sys, sym) || is_parameter(sys, sym) || + all(x -> check_assignable(sys, x), collect(sym)) + elseif sym isa Union{AbstractArray, Tuple} + all(x -> check_assignable(sys, x), sym) + else + false + end +end + function compile_affect(affect::FunctionalAffect, cb, sys, dvs, ps; kwargs...) compile_user_affect(affect, cb, sys, dvs, ps; kwargs...) end - function _compile_optional_affect(default, aff, cb, sys, dvs, ps; kwargs...) if isnothing(aff) || aff == default return nothing diff --git a/src/systems/imperative_affect.jl b/src/systems/imperative_affect.jl new file mode 100644 index 0000000000..2f489913ad --- /dev/null +++ b/src/systems/imperative_affect.jl @@ -0,0 +1,218 @@ + +""" + ImperativeAffect(f::Function; modified::NamedTuple, observed::NamedTuple, ctx) + +`ImperativeAffect` is a helper for writing affect functions that will compute observed values and +ensure that modified values are correctly written back into the system. The affect function `f` needs to have +the signature + +``` + f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple +``` + +The function `f` will be called with `observed` and `modified` `NamedTuple`s that are derived from their respective `NamedTuple` definitions. +Each declaration`NamedTuple` should map an expression to a symbol; for example if we pass `observed=(; x = a + b)` this will alias the result of executing `a+b` in the system as `x` +so the value of `a + b` will be accessible as `observed.x` in `f`. `modified` currently restricts symbolic expressions to only bare variables, so only tuples of the form +`(; x = y)` or `(; x)` (which aliases `x` as itself) are allowed. + +The argument NamedTuples (for instance `(;x=y)`) will be populated with the declared values on function entry; if we require `(;x=y)` in `observed` and `y=2`, for example, +then the NamedTuple `(;x=2)` will be passed as `observed` to the affect function `f`. + +The NamedTuple returned from `f` includes the values to be written back to the system after `f` returns. For example, if we want to update the value of `x` to be the result of `x + y` we could write + + ImperativeAffect(observed=(; x_plus_y = x + y), modified=(; x)) do m, o + @set! m.x = o.x_plus_y + end + +Where we use Setfield to copy the tuple `m` with a new value for `x`, then return the modified value of `m`. All values updated by the tuple must have names originally declared in +`modified`; a runtime error will be produced if a value is written that does not appear in `modified`. The user can dynamically decide not to write a value back by not including it +in the returned tuple, in which case the associated field will not be updated. +""" +@kwdef struct ImperativeAffect + f::Any + obs::Vector + obs_syms::Vector{Symbol} + modified::Vector + mod_syms::Vector{Symbol} + ctx::Any + skip_checks::Bool +end + +function ImperativeAffect(f::Function; + observed::NamedTuple = NamedTuple{()}(()), + modified::NamedTuple = NamedTuple{()}(()), + ctx = nothing, + skip_checks = false) + ImperativeAffect(f, + collect(values(observed)), collect(keys(observed)), + collect(values(modified)), collect(keys(modified)), + ctx, skip_checks) +end +function ImperativeAffect(f::Function, modified::NamedTuple; + observed::NamedTuple = NamedTuple{()}(()), ctx = nothing, skip_checks = false) + ImperativeAffect( + f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) +end +function ImperativeAffect( + f::Function, modified::NamedTuple, observed::NamedTuple; ctx = nothing, skip_checks = false) + ImperativeAffect( + f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) +end +function ImperativeAffect( + f::Function, modified::NamedTuple, observed::NamedTuple, ctx; skip_checks = false) + ImperativeAffect( + f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) +end + +function Base.show(io::IO, mfa::ImperativeAffect) + obs_vals = join(map((ob, nm) -> "$ob => $nm", mfa.obs, mfa.obs_syms), ", ") + mod_vals = join(map((md, nm) -> "$md => $nm", mfa.modified, mfa.mod_syms), ", ") + affect = mfa.f + print(io, + "ImperativeAffect(observed: [$obs_vals], modified: [$mod_vals], affect:$affect)") +end +func(f::ImperativeAffect) = f.f +context(a::ImperativeAffect) = a.ctx +observed(a::ImperativeAffect) = a.obs +observed_syms(a::ImperativeAffect) = a.obs_syms +discretes(a::ImperativeAffect) = filter(ModelingToolkit.isparameter, a.modified) +modified(a::ImperativeAffect) = a.modified +modified_syms(a::ImperativeAffect) = a.mod_syms + +function Base.:(==)(a1::ImperativeAffect, a2::ImperativeAffect) + isequal(a1.f, a2.f) && isequal(a1.obs, a2.obs) && isequal(a1.modified, a2.modified) && + isequal(a1.obs_syms, a2.obs_syms) && isequal(a1.mod_syms, a2.mod_syms) && + isequal(a1.ctx, a2.ctx) +end + +function Base.hash(a::ImperativeAffect, s::UInt) + s = hash(a.f, s) + s = hash(a.obs, s) + s = hash(a.obs_syms, s) + s = hash(a.modified, s) + s = hash(a.mod_syms, s) + hash(a.ctx, s) +end + +namespace_affects(af::ImperativeAffect, s) = namespace_affect(af, s) +function namespace_affect(affect::ImperativeAffect, s) + ImperativeAffect(func(affect), + namespace_expr.(observed(affect), (s,)), + observed_syms(affect), + renamespace.((s,), modified(affect)), + modified_syms(affect), + context(affect), + affect.skip_checks) +end + +function compile_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs...) + compile_user_affect(affect, cb, sys, dvs, ps; kwargs...) +end + +function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs...) + #= + Implementation sketch: + generate observed function (oop), should save to a component array under obs_syms + do the same stuff as the normal FA for pars_syms + call the affect method + unpack and apply the resulting values + =# + function check_dups(syms, exprs) # = (syms_dedup, exprs_dedup) + seen = Set{Symbol}() + syms_dedup = [] + exprs_dedup = [] + for (sym, exp) in Iterators.zip(syms, exprs) + if !in(sym, seen) + push!(syms_dedup, sym) + push!(exprs_dedup, exp) + push!(seen, sym) + elseif !affect.skip_checks + @warn "Expression $(expr) is aliased as $sym, which has already been used. The first definition will be used." + end + end + return (syms_dedup, exprs_dedup) + end + + obs_exprs = observed(affect) + if !affect.skip_checks + for oexpr in obs_exprs + invalid_vars = invalid_variables(sys, oexpr) + if length(invalid_vars) > 0 + error("Observed equation $(oexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing).") + end + end + end + obs_syms = observed_syms(affect) + obs_syms, obs_exprs = check_dups(obs_syms, obs_exprs) + + mod_exprs = modified(affect) + if !affect.skip_checks + for mexpr in mod_exprs + if !check_assignable(sys, mexpr) + @warn ("Expression $mexpr cannot be assigned to; currently only unknowns and parameters may be updated by an affect.") + end + invalid_vars = unassignable_variables(sys, mexpr) + if length(invalid_vars) > 0 + error("Modified equation $(mexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing) or they may have been reduced away.") + end + end + end + mod_syms = modified_syms(affect) + mod_syms, mod_exprs = check_dups(mod_syms, mod_exprs) + + overlapping_syms = intersect(mod_syms, obs_syms) + if length(overlapping_syms) > 0 && !affect.skip_checks + @warn "The symbols $overlapping_syms are declared as both observed and modified; this is a code smell because it becomes easy to confuse them and assign/not assign a value." + end + + # sanity checks done! now build the data and update function for observed values + mkzero(sz) = + if sz === () + 0.0 + else + zeros(sz) + end + obs_fun = build_explicit_observed_function( + sys, Symbolics.scalarize.(obs_exprs); + mkarray = (es, _) -> MakeTuple(es)) + obs_sym_tuple = (obs_syms...,) + + # okay so now to generate the stuff to assign it back into the system + mod_pairs = mod_exprs .=> mod_syms + mod_names = (mod_syms...,) + mod_og_val_fun = build_explicit_observed_function( + sys, Symbolics.scalarize.(first.(mod_pairs)); + mkarray = (es, _) -> MakeTuple(es)) + + upd_funs = NamedTuple{mod_names}((setu.((sys,), first.(mod_pairs))...,)) + + if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing + save_idxs = get(ic.callback_to_clocks, cb, Int[]) + else + save_idxs = Int[] + end + + let user_affect = func(affect), ctx = context(affect) + function (integ) + # update the to-be-mutated values; this ensures that if you do a no-op then nothing happens + modvals = mod_og_val_fun(integ.u, integ.p, integ.t) + upd_component_array = NamedTuple{mod_names}(modvals) + + # update the observed values + obs_component_array = NamedTuple{obs_sym_tuple}(obs_fun( + integ.u, integ.p, integ.t)) + + # let the user do their thing + upd_vals = user_affect(upd_component_array, obs_component_array, ctx, integ) + + # write the new values back to the integrator + _generated_writeback(integ, upd_funs, upd_vals) + + for idx in save_idxs + SciMLBase.save_discretes!(integ, idx) + end + end + end +end + +scalarize_affects(affects::ImperativeAffect) = affects diff --git a/src/systems/index_cache.jl b/src/systems/index_cache.jl index aca6a0547d..4fc2f18bfd 100644 --- a/src/systems/index_cache.jl +++ b/src/systems/index_cache.jl @@ -116,7 +116,7 @@ function IndexCache(sys::AbstractSystem) for affect in affs if affect isa Equation is_parameter(sys, affect.lhs) && push!(discs, affect.lhs) - elseif affect isa FunctionalAffect + elseif affect isa FunctionalAffect || affect isa ImperativeAffect union!(discs, unwrap.(discretes(affect))) else error("Unhandled affect type $(typeof(affect))") diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index 9d3aac4074..ccf0a17a40 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -8,6 +8,7 @@ using ModelingToolkit: SymbolicContinuousCallback, using StableRNGs import SciMLBase using SymbolicIndexingInterface +using Setfield rng = StableRNG(12345) @variables x(t) = 0 @@ -227,6 +228,119 @@ affect_neg = [x ~ 1] @test e[].affect == affect end +@testset "ImperativeAffect constructors" begin + fmfa(o, x, i, c) = nothing + m = ModelingToolkit.ImperativeAffect(fmfa) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test m.obs == [] + @test m.obs_syms == [] + @test m.modified == [] + @test m.mod_syms == [] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect(fmfa, (;)) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test m.obs == [] + @test m.obs_syms == [] + @test m.modified == [] + @test m.mod_syms == [] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect(fmfa, (; x)) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, []) + @test m.obs_syms == [] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:x] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect(fmfa, (; y = x)) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, []) + @test m.obs_syms == [] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:y] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect(fmfa; observed = (; y = x)) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:y] + @test m.modified == [] + @test m.mod_syms == [] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect(fmfa; modified = (; x)) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, []) + @test m.obs_syms == [] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:x] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect(fmfa; modified = (; y = x)) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, []) + @test m.obs_syms == [] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:y] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect(fmfa, (; x), (; x)) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:x] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:x] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect(fmfa, (; y = x), (; y = x)) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:y] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:y] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect( + fmfa; modified = (; y = x), observed = (; y = x)) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:y] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:y] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect( + fmfa; modified = (; y = x), observed = (; y = x), ctx = 3) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:y] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:y] + @test m.ctx === 3 + + m = ModelingToolkit.ImperativeAffect(fmfa, (; x), (; x), 3) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:x] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:x] + @test m.ctx === 3 +end + ## @named sys = ODESystem(eqs, t, continuous_events = [x ~ 1]) @@ -822,7 +936,7 @@ end @named trigsys = ODESystem(eqs, t; continuous_events = [evt1, evt2]) trigsys_ss = structural_simplify(trigsys) prob = ODEProblem(trigsys_ss, [], (0.0, 2π)) - sol = solve(prob, Tsit5()) + sol = solve(prob, Tsit5(); dtmax = 0.01) required_crossings_c1 = [π / 2, 3 * π / 2] required_crossings_c2 = [π / 6, π / 2, 5 * π / 6, 7 * π / 6, 3 * π / 2, 11 * π / 6] @test maximum(abs.(first.(cr1) .- required_crossings_c1)) < 1e-4 @@ -969,6 +1083,173 @@ end @test sol[b] == [2.0, 5.0, 5.0] @test sol[c] == [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0] end +@testset "Heater" begin + @variables temp(t) + params = @parameters furnace_on_threshold=0.5 furnace_off_threshold=0.7 furnace_power=1.0 leakage=0.1 furnace_on::Bool=false + eqs = [ + D(temp) ~ furnace_on * furnace_power - temp^2 * leakage + ] + + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i, c + @set! x.furnace_on = false + end) + furnace_enable = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_on_threshold], + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i, c + @set! x.furnace_on = true + end) + @named sys = ODESystem( + eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) + ss = structural_simplify(sys) + prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + sol = solve(prob, Tsit5(); dtmax = 0.01) + @test all(sol[temp][sol.t .> 1.0] .<= 0.79) && all(sol[temp][sol.t .> 1.0] .>= 0.49) + + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i + @set! x.furnace_on = false + end; initialize = ModelingToolkit.ImperativeAffect(modified = (; + temp)) do x, o, c, i + @set! x.temp = 0.2 + end) + furnace_enable = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_on_threshold], + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i + @set! x.furnace_on = true + end) + @named sys = ODESystem( + eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) + ss = structural_simplify(sys) + prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + sol = solve(prob, Tsit5(); dtmax = 0.01) + @test all(sol[temp][sol.t .> 1.0] .<= 0.79) && all(sol[temp][sol.t .> 1.0] .>= 0.49) + @test all(sol[temp][sol.t .!= 0.0] .<= 0.79) && all(sol[temp][sol.t .!= 0.0] .>= 0.2) +end + +@testset "ImperativeAffect errors and warnings" begin + @variables temp(t) + params = @parameters furnace_on_threshold=0.5 furnace_off_threshold=0.7 furnace_power=1.0 leakage=0.1 furnace_on::Bool=false + eqs = [ + D(temp) ~ furnace_on * furnace_power - temp^2 * leakage + ] + + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.ImperativeAffect( + modified = (; furnace_on), observed = (; furnace_on)) do x, o, c, i + @set! x.furnace_on = false + end) + @named sys = ODESystem(eqs, t, [temp], params; continuous_events = [furnace_off]) + ss = structural_simplify(sys) + @test_logs (:warn, + "The symbols Any[:furnace_on] are declared as both observed and modified; this is a code smell because it becomes easy to confuse them and assign/not assign a value.") prob=ODEProblem( + ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + + @variables tempsq(t) # trivially eliminated + eqs = [tempsq ~ temp^2 + D(temp) ~ furnace_on * furnace_power - temp^2 * leakage] + + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.ImperativeAffect( + modified = (; furnace_on, tempsq), observed = (; furnace_on)) do x, o, c, i + @set! x.furnace_on = false + end) + @named sys = ODESystem( + eqs, t, [temp, tempsq], params; continuous_events = [furnace_off]) + ss = structural_simplify(sys) + @test_throws "refers to missing variable(s)" prob=ODEProblem( + ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + + @parameters not_actually_here + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.ImperativeAffect(modified = (; furnace_on), + observed = (; furnace_on, not_actually_here)) do x, o, c, i + @set! x.furnace_on = false + end) + @named sys = ODESystem( + eqs, t, [temp, tempsq], params; continuous_events = [furnace_off]) + ss = structural_simplify(sys) + @test_throws "refers to missing variable(s)" prob=ODEProblem( + ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.ImperativeAffect(modified = (; furnace_on), + observed = (; furnace_on)) do x, o, c, i + return (; fictional2 = false) + end) + @named sys = ODESystem( + eqs, t, [temp, tempsq], params; continuous_events = [furnace_off]) + ss = structural_simplify(sys) + prob = ODEProblem( + ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + @test_throws "Tried to write back to" solve(prob, Tsit5()) +end + +@testset "Quadrature" begin + @variables theta(t) omega(t) + params = @parameters qA=0 qB=0 hA=0 hB=0 cnt::Int=0 + eqs = [D(theta) ~ omega + omega ~ 1.0] + function decoder(oldA, oldB, newA, newB) + state = (oldA, oldB, newA, newB) + if state == (0, 0, 1, 0) || state == (1, 0, 1, 1) || state == (1, 1, 0, 1) || + state == (0, 1, 0, 0) + return 1 + elseif state == (0, 0, 0, 1) || state == (0, 1, 1, 1) || state == (1, 1, 1, 0) || + state == (1, 0, 0, 0) + return -1 + elseif state == (0, 0, 0, 0) || state == (0, 1, 0, 1) || state == (1, 0, 1, 0) || + state == (1, 1, 1, 1) + return 0 + else + return 0 # err is interpreted as no movement + end + end + qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0], + ModelingToolkit.ImperativeAffect((; qA, hA, hB, cnt), (; qB)) do x, o, c, i + @set! x.hA = x.qA + @set! x.hB = o.qB + @set! x.qA = 1 + @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) + x + end, + affect_neg = ModelingToolkit.ImperativeAffect( + (; qA, hA, hB, cnt), (; qB)) do x, o, c, i + @set! x.hA = x.qA + @set! x.hB = o.qB + @set! x.qA = 0 + @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) + x + end; rootfind = SciMLBase.RightRootFind) + qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π / 2) ~ 0], + ModelingToolkit.ImperativeAffect((; qB, hA, hB, cnt), (; qA)) do x, o, c, i + @set! x.hA = o.qA + @set! x.hB = x.qB + @set! x.qB = 1 + @set! x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) + x + end, + affect_neg = ModelingToolkit.ImperativeAffect( + (; qB, hA, hB, cnt), (; qA)) do x, o, c, i + @set! x.hA = o.qA + @set! x.hB = x.qB + @set! x.qB = 0 + @set! x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) + x + end; rootfind = SciMLBase.RightRootFind) + @named sys = ODESystem( + eqs, t, [theta, omega], params; continuous_events = [qAevt, qBevt]) + ss = structural_simplify(sys) + prob = ODEProblem(ss, [theta => 1e-5], (0.0, pi)) + sol = solve(prob, Tsit5(); dtmax = 0.01) + @test getp(sol, cnt)(sol) == 198 # we get 2 pulses per phase cycle (cos 0 crossing) and we go to 100 cycles; we miss a few due to the initial state +end @testset "Initialization" begin @variables x(t)