Add option to eliminate constants during structural_simplify.

Author: lamorton

src/systems/abstractsystem.jl

@@ -952,16 +952,20 @@ $(SIGNATURES)
 
 Structurally simplify algebraic equations in a system and compute the
 topological sort of the observed equations. When `simplify=true`, the `simplify`
-function will be applied during the tearing process. It also takes kwargs
+function will be applied during the tearing process. When `simplify_constants=true`
+`eliminate_constants` will be applied prior to tearing. It also takes kwargs
 `allow_symbolic=false` and `allow_parameter=true` which limits the coefficient
 types during tearing.
 
 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)
+    if simplify_constants
+        sys = eliminate_constants(sys)
+    end
     state = TearingState(sys)
     has_io = io !== nothing
     has_io && markio!(state, io...)
@@ -979,6 +983,20 @@ function structural_simplify(sys::AbstractSystem, io = nothing; simplify = false
     return has_io ? (sys, input_idxs) : sys
 end
 
+
+""" Replace constants in dynamical equations with their literal values."""
+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...)

test/constants.jl

@@ -12,16 +12,22 @@ eqs = [D(x) ~ a]
 prob = ODEProblem(sys, [0], [0.0, 1.0], [])
 sol = solve(prob, Tsit5())
 
+newsys = eliminate_constants(sys)
+@test isequal(equations(newsys), [D(x) ~ 1])
+
 # Test structural_simplify substitutions & observed values
 eqs = [D(x) ~ 1,
     w ~ a]
 @named sys = ODESystem(eqs)
-simp = structural_simplify(sys);
+simp = structural_simplify(sys, simplify_constants = false);
 @test isequal(simp.substitutions.subs[1], eqs[2])
 @test isequal(equations(simp)[1], eqs[1])
 prob = ODEProblem(simp, [0], [0.0, 1.0], [])
 sol = solve(prob, Tsit5())
 @test sol[w][1] == 1
+# Now eliminate the constants first
+simp = structural_simplify(sys, simplify_constants = true);
+@test isequal(simp.substitutions.subs[1], w ~ 1)
 
 #Constant with units
 @constants β=1 [unit = u"m/s"]
@@ -31,3 +37,8 @@ MT.get_unit(β)
 D = Differential(t)
 eqs = [D(x) ~ β]
 sys = ODESystem(eqs, name = :sys)
+# Note that the equation won't unit-check now
+#    b/c the literal value doesn't have units on it
+#    Testing that units checking is bypassed in the constructors
+simp = structural_simplify(sys)
+@test_throws MT.ValidationError MT.check_units(simp.eqs...)