test: test parameters as initialization unknowns

Author: AayushSabharwal

test/initializationsystem.jl

@@ -1,4 +1,5 @@
 using ModelingToolkit, OrdinaryDiffEq, NonlinearSolve, Test
+using SymbolicIndexingInterface
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 @parameters g
@@ -576,3 +577,132 @@ sol = solve(oprob_2nd_order_2, Rosenbrock23()) # retcode: Success
     sol = solve(prob, Tsit5(), reltol = 1e-4)
     @test all(sol(1.0, idxs = sys.x) .≈ +exp(1)) && all(sol(1.0, idxs = sys.y) .≈ -exp(1))
 end
+
+@testset "Initialization of parameters" begin
+    function test_parameter(prob, sym, val, initialval = zero(val))
+        @test prob.ps[sym] ≈ initialval
+        @test init(prob, Tsit5()).ps[sym] ≈ val
+        @test solve(prob, Tsit5()).ps[sym] ≈ val
+    end
+    function test_initializesystem(sys, u0map, pmap, p, equation)
+        isys = ModelingToolkit.generate_initializesystem(
+            sys; u0map, pmap, guesses = ModelingToolkit.guesses(sys))
+        @test is_variable(isys, p)
+        @test equation in equations(isys) || (0 ~ -equation.rhs) in equations(isys)
+    end
+    @variables x(t) y(t)
+    @parameters p q
+    u0map = Dict(x => 1.0, y => 1.0)
+    pmap = Dict()
+    pmap[q] = 1.0
+    # `missing` default, equation from ODEProblem
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ x * q, D(y) ~ y * p], t; defaults = [p => missing], guesses = [p => 1.0])
+    pmap[p] = 2q
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), pmap)
+    test_parameter(prob, p, 2.0)
+    # `missing` default, provided guess
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ x, p ~ x + y], t; defaults = [p => missing], guesses = [p => 0.0])
+    prob = ODEProblem(sys, u0map, (0.0, 1.0))
+    test_parameter(prob, p, 2.0)
+    test_initializesystem(sys, u0map, pmap, p, 0 ~ p - x - y)
+
+    # `missing` to ODEProblem, equation from default
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ x * q, D(y) ~ y * p], t; defaults = [p => 2q], guesses = [p => 1.0])
+    pmap[p] = missing
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), pmap)
+    test_parameter(prob, p, 2.0)
+    test_initializesystem(sys, u0map, pmap, p, 0 ~ 2q - p)
+    test_parameter(prob2, p, 2.0)
+    # `missing` to ODEProblem, provided guess
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ x, p ~ x + y], t; guesses = [p => 0.0])
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), pmap)
+    test_parameter(prob, p, 2.0)
+    test_initializesystem(sys, u0map, pmap, p, 0 ~ x + y - p)
+
+    # No `missing`, default and guess
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ x * q, D(y) ~ y * p], t; defaults = [p => 2q], guesses = [p => 0.0])
+    delete!(pmap, p)
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), pmap)
+    test_parameter(prob, p, 2.0)
+    test_initializesystem(sys, u0map, pmap, p, 0 ~ 2q - p)
+
+    # Should not be solved for:
+
+    # ODEProblem value with guess, no `missing`
+    @mtkbuild sys = ODESystem([D(x) ~ x * q, D(y) ~ y * p], t; guesses = [p => 0.0])
+    _pmap = merge(pmap, Dict(p => 3q))
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), _pmap)
+    @test prob.ps[p] ≈ 3.0
+    @test prob.f.initializeprob === nothing
+    # Default overridden by ODEProblem, guess provided
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ q * x, D(y) ~ y * p], t; defaults = [p => 2q], guesses = [p => 1.0])
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), _pmap)
+    @test prob.ps[p] ≈ 3.0
+    @test prob.f.initializeprob === nothing
+
+    @mtkbuild sys = ODESystem([D(x) ~ x, p ~ x + y], t; guesses = [p => 0.0])
+    @test_throws ModelingToolkit.MissingParametersError ODEProblem(
+        sys, [x => 1.0, y => 1.0], (0.0, 1.0))
+
+    @testset "Null system" begin
+        @variables x(t) y(t) s(t)
+        @parameters x0 y0
+        @mtkbuild sys = ODESystem([x ~ x0, y ~ y0, s ~ x + y], t; guesses = [y0 => 0.0])
+        prob = ODEProblem(sys, [s => 1.0], (0.0, 1.0), [x0 => 0.3, y0 => missing])
+        test_parameter(prob, y0, 0.7)
+    end
+
+    using ModelingToolkitStandardLibrary.Mechanical.TranslationalModelica: Fixed, Mass,
+                                                                           Spring, Force,
+                                                                           Damper
+    using ModelingToolkitStandardLibrary.Mechanical: TranslationalModelica as TM
+    using ModelingToolkitStandardLibrary.Blocks: Constant
+
+    @named mass = TM.Mass(; m = 1.0, s = 1.0, v = 0.0, a = 0.0)
+    @named fixed = Fixed(; s0 = 0.0)
+    @named spring = Spring(; c = 2.0, s_rel0 = nothing)
+    @named gravity = Force()
+    @named constant = Constant(; k = 9.81)
+    @named damper = TM.Damper(; d = 0.1)
+    @mtkbuild sys = ODESystem(
+        [connect(fixed.flange, spring.flange_a), connect(spring.flange_b, mass.flange_a),
+            connect(mass.flange_a, gravity.flange), connect(constant.output, gravity.f),
+            connect(fixed.flange, damper.flange_a), connect(damper.flange_b, mass.flange_a)],
+        t;
+        systems = [fixed, spring, mass, gravity, constant, damper],
+        guesses = [spring.s_rel0 => 1.0])
+    prob = ODEProblem(sys, [], (0.0, 1.0), [spring.s_rel0 => missing])
+    test_parameter(prob, spring.s_rel0, -3.905)
+end
+
+@testset "Update initializeprob parameters" begin
+    @variables x(t) y(t)
+    @parameters p q
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ x, p ~ x + y], t; guesses = [x => 0.0, p => 0.0])
+    prob = ODEProblem(sys, [y => 1.0], (0.0, 1.0), [p => 3.0])
+    @test prob.f.initializeprob.ps[p] ≈ 3.0
+    @test init(prob, Tsit5())[x] ≈ 2.0
+    prob.ps[p] = 2.0
+    @test prob.f.initializeprob.ps[p] ≈ 3.0
+    @test init(prob, Tsit5())[x] ≈ 1.0
+    ModelingToolkit.defaults(prob.f.sys)[p] = missing
+end
+
+@testset "Equations for dependent parameters" begin
+    @variables x(t)
+    @parameters p q=5 r
+    @mtkbuild sys = ODESystem(
+        D(x) ~ 2x + r, t; parameter_dependencies = [r ~ p + 2q, q ~ p + 3],
+        guesses = [p => 1.0])
+    prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [p => missing])
+    @test length(equations(ModelingToolkit.get_parent(prob.f.initializeprob.f.sys))) == 4
+    integ = init(prob, Tsit5())
+    @test integ.ps[p] ≈ 2
+end

test/parameter_dependencies.jl

@@ -10,7 +10,7 @@ using SymbolicIndexingInterface
 using NonlinearSolve
 
 @testset "ODESystem with callbacks" begin
-    @parameters p1=1.0 p2=1.0
+    @parameters p1=1.0 p2
     @variables x(t)
     cb1 = [x ~ 2.0] => [p1 ~ 2.0] # triggers at t=-2+√6
     function affect1!(integ, u, p, ctx)
@@ -31,7 +31,7 @@ using NonlinearSolve
     prob = ODEProblem(sys, [x => 1.0], (0.0, 1.5), jac = true)
     @test prob.ps[p1] == 1.0
     @test prob.ps[p2] == 2.0
-    @test_nowarn solve(prob, Tsit5())
+    @test SciMLBase.successful_retcode(solve(prob, Tsit5()))
     prob = ODEProblem(sys, [x => 1.0], (0.0, 1.5), [p1 => 1.0], jac = true)
     @test prob.ps[p1] == 1.0
     @test prob.ps[p2] == 2.0