Merge pull request #1821 from lamorton/constants2

Adding @Constants

Author: YingboMa

docs/src/basics/ContextualVariables.md

@@ -20,6 +20,21 @@ All modeling projects have some form of parameters. `@parameters` marks a variab
 as being the parameter of some system, which allows automatic detection algorithms
 to ignore such variables when attempting to find the states of a system.
 
+## Constants
+
+Constants are like parameters that:
+- always have a default value, which must be assigned when the constants are
+    declared
+- do not show up in the list of parameters of a system.
+
+The intended use-cases for constants are:
+- representing literals (eg, π) symbolically, which results in cleaner
+    Latexification of equations (avoids turning `d ~ 2π*r` into `d = 6.283185307179586 r`)
+- allowing auto-generated unit conversion factors to live outside the list of
+    parameters
+- representing fundamental constants (eg, speed of light `c`) that should never
+     be adjusted inadvertently.
+
 ## Wildcard Variable Arguments
 
 ```julia
@@ -28,7 +43,7 @@ to ignore such variables when attempting to find the states of a system.
 
 It is possible to define a dependent variable which is an open function as above,
 for which its arguments must be specified each time it is used. This is useful with
-PDEs for example, where one may need to use `u(t, x)` in the equations, but will 
+PDEs for example, where one may need to use `u(t, x)` in the equations, but will
 need to be able to write `u(t, 0.0)` to define a boundary condition at `x = 0`.
 
 ## Variable metadata [Experimental/TODO]

docs/src/tutorials/ode_modeling.md

@@ -14,16 +14,18 @@ But if you want to just see some code and run, here's an example:
 using ModelingToolkit
 
 @variables t x(t)   # independent and dependent variables
-@parameters τ       # parameters
+@parameters τ       # parameters 
+@constants h = 1    # constants have an assigned value
 D = Differential(t) # define an operator for the differentiation w.r.t. time
 
 # your first ODE, consisting of a single equation, the equality indicated by ~
-@named fol = ODESystem([ D(x)  ~ (1 - x)/τ])
+@named fol = ODESystem([ D(x)  ~ (h - x)/τ])
 
 using DifferentialEquations: solve
 using Plots: plot
 
 prob = ODEProblem(fol, [x => 0.0], (0.0,10.0), [τ => 3.0])
+# parameter `τ` can be assigned a value, but constant `h` cannot
 sol = solve(prob)
 plot(sol)
 ```
@@ -42,19 +44,20 @@ first-order lag element:
 ```
 
 Here, ``t`` is the independent variable (time), ``x(t)`` is the (scalar) state
-variable, ``f(t)`` is an external forcing function, and ``\tau`` is a constant
+variable, ``f(t)`` is an external forcing function, and ``\tau`` is a 
 parameter. In MTK, this system can be modelled as follows. For simplicity, we
-first set the forcing function to a constant value.
+first set the forcing function to a time-independent value.
 
 ```julia
 using ModelingToolkit
 
 @variables t x(t)  # independent and dependent variables
 @parameters τ       # parameters
+@constants h = 1    # constants
 D = Differential(t) # define an operator for the differentiation w.r.t. time
 
 # your first ODE, consisting of a single equation, indicated by ~
-@named fol_model = ODESystem(D(x) ~ (1 - x)/τ)
+@named fol_model = ODESystem(D(x) ~ (h - x)/τ)
       # Model fol_model with 1 equations
       # States (1):
       #   x(t)
@@ -89,7 +92,7 @@ intermediate variable `RHS`:
 
 ```julia
 @variables RHS(t)
-@named fol_separate = ODESystem([ RHS  ~ (1 - x)/τ,
+@named fol_separate = ODESystem([ RHS  ~ (h - x)/τ,
                                   D(x) ~ RHS ])
       # Model fol_separate with 2 equations
       # States (2):
@@ -110,7 +113,7 @@ fol_simplified = structural_simplify(fol_separate)
 
 equations(fol_simplified)
       # 1-element Array{Equation,1}:
-      #  Differential(t)(x(t)) ~ (τ^-1)*(1 - x(t))
+      #  Differential(t)(x(t)) ~ (τ^-1)*(h - x(t))
 
 equations(fol_simplified) == equations(fol_model)
       # true
@@ -133,6 +136,12 @@ sol = solve(prob)
 plot(sol, vars=[x, RHS])
 ```
 
+By default, `structural_simplify` also replaces symbolic `constants` with
+their default values. This allows additional simplifications not possible 
+if using `parameters` (eg, solution of linear equations by dividing out
+the constant's value, which cannot be done for parameters since they may
+be zero). 
+
 ![Simulation result of first-order lag element, with right-hand side](https://user-images.githubusercontent.com/13935112/111958403-7e8d3e00-8aed-11eb-9d18-08b5180a59f9.png)
 
 Note that similarly the indexing of the solution works via the names, and so

src/ModelingToolkit.jl

@@ -118,6 +118,7 @@ using .BipartiteGraphs
 
 include("variables.jl")
 include("parameters.jl")
+include("constants.jl")
 
 include("utils.jl")
 include("domains.jl")
@@ -214,7 +215,8 @@ export toexpr, get_variables
 export simplify, substitute
 export build_function
 export modelingtoolkitize
-export @variables, @parameters
+
+export @variables, @parameters, @constants
 export @named, @nonamespace, @namespace, extend, compose, complete
 export debug_system
 

src/constants.jl

@@ -0,0 +1,34 @@
+import SymbolicUtils: symtype, term, hasmetadata, issym
+struct MTKConstantCtx end
+
+isconstant(x::Num) = isconstant(unwrap(x))
+""" Test whether `x` is a constant-type Sym. """
+function isconstant(x)
+    x = unwrap(x)
+    x isa Symbolic && getmetadata(x, MTKConstantCtx, false)
+end
+
+"""
+    toconstant(s::Sym)
+
+Maps the parameter to a constant. The parameter must have a default.
+"""
+function toconstant(s::Sym)
+    hasmetadata(s, Symbolics.VariableDefaultValue) ||
+        throw(ArgumentError("Constant `$(s)` must be assigned a default value."))
+    setmetadata(s, MTKConstantCtx, true)
+end
+
+toconstant(s::Num) = wrap(toconstant(value(s)))
+
+"""
+$(SIGNATURES)
+
+Define one or more constants.
+"""
+macro constants(xs...)
+    Symbolics._parse_vars(:constants,
+                          Real,
+                          xs,
+                          toconstant) |> esc
+end

src/structural_transformation/codegen.jl

@@ -1,6 +1,6 @@
 using LinearAlgebra
 
-using ModelingToolkit: isdifferenceeq, process_events
+using ModelingToolkit: isdifferenceeq, process_events, get_preprocess_constants
 
 const MAX_INLINE_NLSOLVE_SIZE = 8
 
@@ -187,12 +187,14 @@ function gen_nlsolve!(is_not_prepended_assignment, eqs, vars, u0map::AbstractDic
 
     fname = gensym("fun")
     # f is the function to find roots on
+    funex = isscalar ? rhss[1] : MakeArray(rhss, SVector)
+    pre = get_preprocess_constants(funex)
     f = Func([DestructuredArgs(vars, inbounds = !checkbounds)
               DestructuredArgs(params, inbounds = !checkbounds)],
              [],
-             Let(needed_assignments[inner_idxs],
-                 isscalar ? rhss[1] : MakeArray(rhss, SVector),
-                 false)) |> SymbolicUtils.Code.toexpr
+             pre(Let(needed_assignments[inner_idxs],
+                     funex,
+                     false))) |> SymbolicUtils.Code.toexpr
 
     # solver call contains code to call the root-finding solver on the function f
     solver_call = LiteralExpr(quote
@@ -294,6 +296,8 @@ function build_torn_function(sys;
     syms = map(Symbol, states)
 
     pre = get_postprocess_fbody(sys)
+    cpre = get_preprocess_constants(rhss)
+    pre2 = x -> pre(cpre(x))
 
     expr = SymbolicUtils.Code.toexpr(Func([out
                                            DestructuredArgs(states,
@@ -302,10 +306,10 @@ function build_torn_function(sys;
                                                             inbounds = !checkbounds)
                                            independent_variables(sys)],
                                           [],
-                                          pre(Let([torn_expr;
-                                                   assignments[is_not_prepended_assignment]],
-                                                  funbody,
-                                                  false))),
+                                          pre2(Let([torn_expr;
+                                                    assignments[is_not_prepended_assignment]],
+                                                   funbody,
+                                                   false))),
                                      sol_states)
     if expression
         expr, states
@@ -477,17 +481,19 @@ function build_observed_function(state, ts, var_eq_matching, var_sccs,
         push!(subs, sym ← obs[eqidx].rhs)
     end
     pre = get_postprocess_fbody(sys)
-
+    cpre = get_preprocess_constants([obs[1:maxidx];
+                                     isscalar ? ts[1] : MakeArray(ts, output_type)])
+    pre2 = x -> pre(cpre(x))
     ex = Code.toexpr(Func([DestructuredArgs(solver_states, inbounds = !checkbounds)
                            DestructuredArgs(parameters(sys), inbounds = !checkbounds)
                            independent_variables(sys)],
                           [],
-                          pre(Let([collect(Iterators.flatten(solves))
-                                   assignments[is_not_prepended_assignment]
-                                   map(eq -> eq.lhs ← eq.rhs, obs[1:maxidx])
-                                   subs],
-                                  isscalar ? ts[1] : MakeArray(ts, output_type),
-                                  false))), sol_states)
+                          pre2(Let([collect(Iterators.flatten(solves))
+                                    assignments[is_not_prepended_assignment]
+                                    map(eq -> eq.lhs ← eq.rhs, obs[1:maxidx])
+                                    subs],
+                                   isscalar ? ts[1] : MakeArray(ts, output_type),
+                                   false))), sol_states)
 
     expression ? ex : @RuntimeGeneratedFunction(ex)
 end

src/structural_transformation/utils.jl

@@ -179,6 +179,8 @@ function find_eq_solvables!(state::TearingState, ieq, to_rm = Int[], coeffs = no
         a, b, islinear = linear_expansion(term, var)
         a, b = unwrap(a), unwrap(b)
         islinear || (all_int_vars = false; continue)
+        a = ModelingToolkit.fold_constants(a)
+        b = ModelingToolkit.fold_constants(b)
         if a isa Symbolic
             all_int_vars = false
             if !allow_symbolic

src/systems/abstractsystem.jl

@@ -1030,7 +1030,8 @@ The optional argument `io` may take a tuple `(inputs, outputs)`.
 This will convert all `inputs` to parameters and allow them to be unconnected, i.e.,
 simplification will allow models where `n_states = n_equations - n_inputs`.
 """
-function structural_simplify(sys::AbstractSystem, io = nothing; simplify = false, kwargs...)
+function structural_simplify(sys::AbstractSystem, io = nothing; simplify = false,
+                             simplify_constants = true, kwargs...)
     sys = expand_connections(sys)
     sys isa DiscreteSystem && return sys
     state = TearingState(sys)
@@ -1046,6 +1047,18 @@ function structural_simplify(sys::AbstractSystem, io = nothing; simplify = false
     return has_io ? (sys, input_idxs) : sys
 end
 
+function eliminate_constants(sys::AbstractSystem)
+    if has_eqs(sys)
+        eqs = get_eqs(sys)
+        eq_cs = collect_constants(eqs)
+        if !isempty(eq_cs)
+            new_eqs = eliminate_constants(eqs, eq_cs)
+            @set! sys.eqs = new_eqs
+        end
+    end
+    return sys
+end
+
 function io_preprocessing(sys::AbstractSystem, inputs,
                           outputs; simplify = false, kwargs...)
     sys, input_idxs = structural_simplify(sys, (; inputs, outputs); simplify, kwargs...)

src/systems/callbacks.jl

@@ -269,7 +269,13 @@ function compile_condition(cb::SymbolicDiscreteCallback, sys, dvs, ps;
     p = map(x -> time_varying_as_func(value(x), sys), ps)
     t = get_iv(sys)
     condit = condition(cb)
-    build_function(condit, u, t, p; expression, wrap_code = condition_header(), kwargs...)
+    cs = collect_constants(condit)
+    if !isempty(cs)
+        cmap = map(x -> x => getdefault(x), cs)
+        condit = substitute(condit, cmap)
+    end
+    build_function(condit, u, t, p; expression, wrap_code = condition_header(),
+                   kwargs...)
 end
 
 function compile_affect(cb::SymbolicContinuousCallback, args...; kwargs...)
@@ -337,9 +343,11 @@ function compile_affect(eqs::Vector{Equation}, sys, dvs, ps; outputidxs = nothin
         t = get_iv(sys)
         integ = gensym(:MTKIntegrator)
         getexpr = (postprocess_affect_expr! === nothing) ? expression : Val{true}
+        pre = get_preprocess_constants(rhss)
         rf_oop, rf_ip = build_function(rhss, u, p, t; expression = getexpr,
                                        wrap_code = add_integrator_header(integ, outvar),
                                        outputidxs = update_inds,
+                                       postprocess_fbody = pre,
                                        kwargs...)
         # applied user-provided function to the generated expression
         if postprocess_affect_expr! !== nothing
@@ -376,7 +384,9 @@ function generate_rootfinding_callback(cbs, sys::AbstractODESystem, dvs = states
     u = map(x -> time_varying_as_func(value(x), sys), dvs)
     p = map(x -> time_varying_as_func(value(x), sys), ps)
     t = get_iv(sys)
-    rf_oop, rf_ip = build_function(rhss, u, p, t; expression = Val{false}, kwargs...)
+    pre = get_preprocess_constants(rhss)
+    rf_oop, rf_ip = build_function(rhss, u, p, t; expression = Val{false},
+                                   postprocess_fbody = pre, kwargs...)
 
     affect_functions = map(cbs) do cb # Keep affect function separate
         eq_aff = affects(cb)

src/systems/diffeqs/abstractodesystem.jl

@@ -83,13 +83,15 @@ end
 function generate_tgrad(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys);
                         simplify = false, kwargs...)
     tgrad = calculate_tgrad(sys, simplify = simplify)
-    return build_function(tgrad, dvs, ps, get_iv(sys); kwargs...)
+    pre = get_preprocess_constants(tgrad)
+    return build_function(tgrad, dvs, ps, get_iv(sys); postprocess_fbody = pre, kwargs...)
 end
 
 function generate_jacobian(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys);
                            simplify = false, sparse = false, kwargs...)
     jac = calculate_jacobian(sys; simplify = simplify, sparse = sparse)
-    return build_function(jac, dvs, ps, get_iv(sys); kwargs...)
+    pre = get_preprocess_constants(jac)
+    return build_function(jac, dvs, ps, get_iv(sys); postprocess_fbody = pre, kwargs...)
 end
 
 function generate_control_jacobian(sys::AbstractODESystem, dvs = states(sys),
@@ -109,7 +111,9 @@ function generate_dae_jacobian(sys::AbstractODESystem, dvs = states(sys),
     dvs = states(sys)
     @variables ˍ₋gamma
     jac = ˍ₋gamma * jac_du + jac_u
-    return build_function(jac, derivatives, dvs, ps, ˍ₋gamma, get_iv(sys); kwargs...)
+    pre = get_preprocess_constants(jac)
+    return build_function(jac, derivatives, dvs, ps, ˍ₋gamma, get_iv(sys);
+                          postprocess_fbody = pre, kwargs...)
 end
 
 function generate_function(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys);
@@ -163,8 +167,10 @@ function generate_difference_cb(sys::ODESystem, dvs = states(sys), ps = paramete
     end
 
     pre = get_postprocess_fbody(sys)
+    cpre = get_preprocess_constants(body)
+    pre2 = x -> pre(cpre(x))
     f_oop, f_iip = build_function(body, u, p, t; expression = Val{false},
-                                  postprocess_fbody = pre, kwargs...)
+                                  postprocess_fbody = pre2, kwargs...)
 
     cb_affect! = let f_oop = f_oop, f_iip = f_iip
         function cb_affect!(integ)

src/systems/diffeqs/odesystem.jl

@@ -206,7 +206,7 @@ end
 
 function ODESystem(eqs, iv = nothing; kwargs...)
     eqs = scalarize(eqs)
-    # NOTE: this assumes that the order of algebric equations doesn't matter
+    # NOTE: this assumes that the order of algebraic equations doesn't matter
     diffvars = OrderedSet()
     allstates = OrderedSet()
     ps = OrderedSet()
@@ -301,6 +301,13 @@ function build_explicit_observed_function(sys, ts;
     dep_vars = scalarize(setdiff(vars, ivs))
 
     obs = observed(sys)
+
+    cs = collect_constants(obs)
+    if !isempty(cs) > 0
+        cmap = map(x -> x => getdefault(x), cs)
+        obs = map(x -> x.lhs ~ substitute(x.rhs, cmap), obs)
+    end
+
     sts = Set(states(sys))
     observed_idx = Dict(x.lhs => i for (i, x) in enumerate(obs))
     namespaced_to_obs = Dict(states(sys, x.lhs) => x.lhs for x in obs)