Merge branch 'master' into extend_metadata

Author: hersle

docs/src/tutorials/nonlinear.md

@@ -9,31 +9,27 @@ We use (unknown) variables for our nonlinear system.
 ```@example nonlinear
 using ModelingToolkit, NonlinearSolve
 
+# Define a nonlinear system
 @variables x y z
 @parameters σ ρ β
+eqs = [0 ~ σ * (y - x)
+       0 ~ x * (ρ - z) - y
+       0 ~ x * y - β * z]
+guesses = [x => 1.0, y => 0.0, z => 0.0]
+ps = [σ => 10.0, ρ => 26.0, β => 8 / 3]
+@mtkbuild ns = NonlinearSystem(eqs)
 
-# Define a nonlinear system
-eqs = [0 ~ σ * (y - x),
-    0 ~ x * (ρ - z) - y,
-    0 ~ x * y - β * z]
-@mtkbuild ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])
-
-guess = [x => 1.0,
-    y => 0.0,
-    z => 0.0]
-
-ps = [σ => 10.0
-      ρ => 26.0
-      β => 8 / 3]
+guesses = [x => 1.0, y => 0.0, z => 0.0]
+ps = [σ => 10.0, ρ => 26.0, β => 8 / 3]
 
-prob = NonlinearProblem(ns, guess, ps)
+prob = NonlinearProblem(ns, guesses, ps)
 sol = solve(prob, NewtonRaphson())
 ```
 
 We can similarly ask to generate the `NonlinearProblem` with the analytical
 Jacobian function:
 
 ```@example nonlinear
-prob = NonlinearProblem(ns, guess, ps, jac = true)
+prob = NonlinearProblem(ns, guesses, ps, jac = true)
 sol = solve(prob, NewtonRaphson())
 ```

src/systems/diffeqs/abstractodesystem.jl

@@ -725,10 +725,16 @@ function get_u0(
         defs = mergedefaults(defs, parammap, ps)
     end
 
-    obs = filter!(x -> !(x[1] isa Number),
-        map(x -> isparameter(x.rhs) ? x.lhs => x.rhs : x.rhs => x.lhs, observed(sys)))
-    observedmap = isempty(obs) ? Dict() : todict(obs)
-    defs = mergedefaults(defs, observedmap, u0map, dvs)
+    # Convert observed equations "lhs ~ rhs" into defaults.
+    # Use the order "lhs => rhs" by default, but flip it to "rhs => lhs"
+    # if "lhs" is known by other means (parameter, another default, ...)
+    # TODO: Is there a better way to determine which equations to flip?
+    obs = map(x -> x.lhs => x.rhs, observed(sys))
+    obs = map(x -> x[1] in keys(defs) ? reverse(x) : x, obs)
+    obs = filter!(x -> !(x[1] isa Number), obs) # exclude e.g. "0 => x^2 + y^2 - 25"
+    obsmap = isempty(obs) ? Dict() : todict(obs)
+
+    defs = mergedefaults(defs, obsmap, u0map, dvs)
     if symbolic_u0
         u0 = varmap_to_vars(
             u0map, dvs; defaults = defs, tofloat = false, use_union = false, toterm)

src/systems/nonlinear/nonlinearsystem.jl

@@ -193,8 +193,12 @@ function calculate_jacobian(sys::NonlinearSystem; sparse = false, simplify = fal
         return cache[1]
     end
 
-    rhs = [eq.rhs for eq in equations(sys)]
+    # observed equations may depend on unknowns, so substitute them in first
+    # TODO: rather keep observed derivatives unexpanded, like "Differential(obs)(expr)"?
+    obs = Dict(eq.lhs => eq.rhs for eq in observed(sys))
+    rhs = map(eq -> fixpoint_sub(eq.rhs, obs), equations(sys))
     vals = [dv for dv in unknowns(sys)]
+
     if sparse
         jac = sparsejacobian(rhs, vals, simplify = simplify)
     else
@@ -215,7 +219,8 @@ function generate_jacobian(
 end
 
 function calculate_hessian(sys::NonlinearSystem; sparse = false, simplify = false)
-    rhs = [eq.rhs for eq in equations(sys)]
+    obs = Dict(eq.lhs => eq.rhs for eq in observed(sys))
+    rhs = map(eq -> fixpoint_sub(eq.rhs, obs), equations(sys))
     vals = [dv for dv in unknowns(sys)]
     if sparse
         hess = [sparsehessian(rhs[i], vals, simplify = simplify) for i in 1:length(rhs)]

src/utils.jl

@@ -228,7 +228,7 @@ function iv_from_nested_derivative(x, op = Differential)
 end
 
 hasdefault(v) = hasmetadata(v, Symbolics.VariableDefaultValue)
-getdefault(v) = value(getmetadata(v, Symbolics.VariableDefaultValue))
+getdefault(v) = value(Symbolics.getdefaultval(v))
 function getdefaulttype(v)
     def = value(getmetadata(unwrap(v), Symbolics.VariableDefaultValue, nothing))
     def === nothing ? Float64 : typeof(def)

test/guess_propagation.jl

@@ -90,21 +90,21 @@ prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0])
 @variables y(t) = x0
 @mtkbuild sys = ODESystem([x ~ x0, D(y) ~ x], t)
 prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0])
-prob[x] == 1.0
-prob[y] == 1.0
+@test prob[x] == 1.0
+@test prob[y] == 1.0
 
 @parameters x0
 @variables x(t)
 @variables y(t) = x0
 @mtkbuild sys = ODESystem([x ~ y, D(y) ~ x], t)
 prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0])
-prob[x] == 1.0
-prob[y] == 1.0
+@test prob[x] == 1.0
+@test prob[y] == 1.0
 
 @parameters x0
 @variables x(t) = x0
 @variables y(t) = x
 @mtkbuild sys = ODESystem([x ~ y, D(y) ~ x], t)
 prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0])
-prob[x] == 1.0
-prob[y] == 1.0
+@test prob[x] == 1.0
+@test prob[y] == 1.0

test/model_parsing.jl

@@ -263,7 +263,7 @@ end
     @test getdefault(model.cval) == 1
     @test isequal(getdefault(model.c), model.cval + model.jval)
     @test getdefault(model.d) == 2
-    @test_throws KeyError getdefault(model.e)
+    @test_throws ErrorException getdefault(model.e)
     @test getdefault(model.f) == 3
     @test getdefault(model.i) == 4
     @test all(getdefault.(scalarize(model.b2)) .== [1, 3])

test/nonlinearsystem.jl

@@ -283,3 +283,37 @@ sys = structural_simplify(ns)
 @test length(equations(sys)) == 0
 sys = structural_simplify(ns; conservative = true)
 @test length(equations(sys)) == 1
+
+# https://github.com/SciML/ModelingToolkit.jl/issues/2858
+@testset "Jacobian/Hessian with observed equations that depend on unknowns" begin
+    @variables x y z
+    @parameters σ ρ β
+    eqs = [0 ~ σ * (y - x)
+           0 ~ x * (ρ - z) - y
+           0 ~ x * y - β * z]
+    guesses = [x => 1.0, y => 0.0, z => 0.0]
+    ps = [σ => 10.0, ρ => 26.0, β => 8 / 3]
+    @mtkbuild ns = NonlinearSystem(eqs)
+
+    @test isequal(calculate_jacobian(ns), [(-1 - z + ρ)*σ -x*σ
+                                           2x*(-z + ρ) -β-(x^2)])
+    # solve without analytical jacobian
+    prob = NonlinearProblem(ns, guesses, ps)
+    sol = solve(prob, NewtonRaphson())
+    @test sol.retcode == ReturnCode.Success
+
+    # solve with analytical jacobian
+    prob = NonlinearProblem(ns, guesses, ps, jac = true)
+    sol = solve(prob, NewtonRaphson())
+    @test sol.retcode == ReturnCode.Success
+
+    # system that contains a chain of observed variables when simplified
+    @variables x y z
+    eqs = [0 ~ x^2 + 2z + y, z ~ y, y ~ x] # analytical solution x = y = z = 0 or -3
+    @mtkbuild ns = NonlinearSystem(eqs) # solve for y with observed chain z -> x -> y
+    @test isequal(expand.(calculate_jacobian(ns)), [3 // 2 + y;;])
+    @test isequal(calculate_hessian(ns), [[1;;]])
+    prob = NonlinearProblem(ns, unknowns(ns) .=> -4.0) # give guess < -3 to reach -3
+    sol = solve(prob, NewtonRaphson())
+    @test sol[x] ≈ sol[y] ≈ sol[z] ≈ -3
+end

test/odesystem.jl

@@ -1208,3 +1208,23 @@ end
     @test     ModelingToolkit.get_metadata(extend(A2, B1))  == A
     @test Set(ModelingToolkit.get_metadata(extend(A2, B2))) == Set(A ∪ B)
 end
+
+# https://github.com/SciML/ModelingToolkit.jl/issues/2859
+@testset "Initialization with defaults from observed equations (edge case)" begin
+    @variables x(t) y(t) z(t)
+    eqs = [D(x) ~ 0, y ~ x, D(z) ~ 0]
+    defaults = [x => 1, z => y]
+    @named sys = ODESystem(eqs, t; defaults)
+    ssys = structural_simplify(sys)
+    prob = ODEProblem(ssys, [], (0.0, 1.0), [])
+    @test prob[x] == prob[y] == prob[z] == 1.0
+
+    @parameters y0
+    @variables x(t) y(t) z(t)
+    eqs = [D(x) ~ 0, y ~ y0 / x, D(z) ~ y]
+    defaults = [y0 => 1, x => 1, z => y]
+    @named sys = ODESystem(eqs, t; defaults)
+    ssys = structural_simplify(sys)
+    prob = ODEProblem(ssys, [], (0.0, 1.0), [])
+    @test prob[x] == prob[y] == prob[z] == 1.0
+end

test/variable_parsing.jl

@@ -104,16 +104,18 @@ end
     y = 2, [connect = Flow]
 end
 
+@test_throws ErrorException ModelingToolkit.getdefault(x)
 @test !hasmetadata(x, VariableDefaultValue)
 @test getmetadata(x, VariableConnectType) == Flow
 @test getmetadata(x, VariableUnit) == u
 @test getmetadata(y, VariableDefaultValue) === 2
 @test getmetadata(y, VariableConnectType) == Flow
 
 a = rename(value(x), :a)
-@test !hasmetadata(x, VariableDefaultValue)
-@test getmetadata(x, VariableConnectType) == Flow
-@test getmetadata(x, VariableUnit) == u
+@test_throws ErrorException ModelingToolkit.getdefault(a)
+@test !hasmetadata(a, VariableDefaultValue)
+@test getmetadata(a, VariableConnectType) == Flow
+@test getmetadata(a, VariableUnit) == u
 
 @variables t x(t)=1 [connect = Flow, unit = u]