From f736a7857a883f4e46eee2ed0b6d79ee2210a63c Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Tue, 16 Jul 2024 16:36:28 +0200 Subject: [PATCH 01/12] Substitute in observed before calculating jacobian --- src/systems/nonlinear/nonlinearsystem.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl index db13bb4541..445f01f10e 100644 --- a/src/systems/nonlinear/nonlinearsystem.jl +++ b/src/systems/nonlinear/nonlinearsystem.jl @@ -193,8 +193,11 @@ 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 + obs = map(eq -> eq.lhs => eq.rhs, observed(sys)) + rhs = map(eq -> substitute(eq.rhs, obs), equations(sys)) vals = [dv for dv in unknowns(sys)] + if sparse jac = sparsejacobian(rhs, vals, simplify = simplify) else From 890d8593ad8def87fc9e202cd98bd8961c607aeb Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Tue, 16 Jul 2024 16:52:32 +0200 Subject: [PATCH 02/12] Added test --- test/nonlinearsystem.jl | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/test/nonlinearsystem.jl b/test/nonlinearsystem.jl index 75c94a6d60..099d1726bc 100644 --- a/test/nonlinearsystem.jl +++ b/test/nonlinearsystem.jl @@ -283,3 +283,32 @@ 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 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 +end From 15c44d5295ab5d18230534c6b022bda97e1550b1 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Tue, 16 Jul 2024 17:06:24 +0200 Subject: [PATCH 03/12] Clean up example a bit --- docs/src/tutorials/nonlinear.md | 25 ++++++++++--------------- 1 file changed, 10 insertions(+), 15 deletions(-) diff --git a/docs/src/tutorials/nonlinear.md b/docs/src/tutorials/nonlinear.md index 2043bd4be1..3ccc1d9861 100644 --- a/docs/src/tutorials/nonlinear.md +++ b/docs/src/tutorials/nonlinear.md @@ -9,24 +9,19 @@ We use (unknown) variables for our nonlinear system. ```@example nonlinear using ModelingToolkit, NonlinearSolve +# Define a nonlinear system @variables x y z @parameters σ ρ β +@mtkbuild ns = NonlinearSystem([ + 0 ~ σ * (y - x) + 0 ~ x * (ρ - z) - y + 0 ~ x * y - β * z +]) -# 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()) ``` @@ -34,6 +29,6 @@ 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()) ``` From b477f736b30bee9e41ca7b9e95490cf21c190805 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Tue, 16 Jul 2024 17:12:13 +0200 Subject: [PATCH 04/12] Add two TODOs --- src/systems/nonlinear/nonlinearsystem.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl index 445f01f10e..7459f3cd8b 100644 --- a/src/systems/nonlinear/nonlinearsystem.jl +++ b/src/systems/nonlinear/nonlinearsystem.jl @@ -194,6 +194,8 @@ function calculate_jacobian(sys::NonlinearSystem; sparse = false, simplify = fal end # observed equations may depend on unknowns, so substitute them in first + # TODO: must do the same fix in e.g. calculate_hessian? + # TODO: rather keep observed derivatives unexpanded, like "Differential(obs)(expr)"? obs = map(eq -> eq.lhs => eq.rhs, observed(sys)) rhs = map(eq -> substitute(eq.rhs, obs), equations(sys)) vals = [dv for dv in unknowns(sys)] From 597a56bca04be92b1002d3b9c8ef03c9209c9a87 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Tue, 16 Jul 2024 17:31:26 +0200 Subject: [PATCH 05/12] Format --- test/nonlinearsystem.jl | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/test/nonlinearsystem.jl b/test/nonlinearsystem.jl index 099d1726bc..fb4acb618a 100644 --- a/test/nonlinearsystem.jl +++ b/test/nonlinearsystem.jl @@ -288,20 +288,15 @@ sys = structural_simplify(ns; conservative = true) @testset "Jacobian 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 - ] + 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) - ]) - + @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()) From dbd2dc5ce96f0a6700c4111b2e179a6e070bc435 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Tue, 16 Jul 2024 17:34:38 +0200 Subject: [PATCH 06/12] Clean up my own mess --- docs/src/tutorials/nonlinear.md | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/docs/src/tutorials/nonlinear.md b/docs/src/tutorials/nonlinear.md index 3ccc1d9861..fc525cc988 100644 --- a/docs/src/tutorials/nonlinear.md +++ b/docs/src/tutorials/nonlinear.md @@ -12,11 +12,12 @@ using ModelingToolkit, NonlinearSolve # Define a nonlinear system @variables x y z @parameters σ ρ β -@mtkbuild ns = NonlinearSystem([ - 0 ~ σ * (y - x) - 0 ~ x * (ρ - z) - y - 0 ~ x * y - β * z -]) +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) guesses = [x => 1.0, y => 0.0, z => 0.0] ps = [σ => 10.0, ρ => 26.0, β => 8 / 3] From 7d9d1fda64c8ac040dbb5892bb2463a63065f18b Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Tue, 16 Jul 2024 17:56:35 +0200 Subject: [PATCH 07/12] Use fixpoint_sub --- src/systems/nonlinear/nonlinearsystem.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl index 7459f3cd8b..e1ad9959af 100644 --- a/src/systems/nonlinear/nonlinearsystem.jl +++ b/src/systems/nonlinear/nonlinearsystem.jl @@ -196,8 +196,8 @@ function calculate_jacobian(sys::NonlinearSystem; sparse = false, simplify = fal # observed equations may depend on unknowns, so substitute them in first # TODO: must do the same fix in e.g. calculate_hessian? # TODO: rather keep observed derivatives unexpanded, like "Differential(obs)(expr)"? - obs = map(eq -> eq.lhs => eq.rhs, observed(sys)) - rhs = map(eq -> substitute(eq.rhs, obs), equations(sys)) + obs = map(eq -> eq.lhs => eq.rhs, observed(sys)) |> todict + rhs = map(eq -> fixpoint_sub(eq.rhs, obs), equations(sys)) vals = [dv for dv in unknowns(sys)] if sparse From d68cc621bacf7cad08f01abaf95964ada55407db Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Tue, 16 Jul 2024 19:23:59 +0200 Subject: [PATCH 08/12] Add test that requires substitution of chained observeds --- test/nonlinearsystem.jl | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/test/nonlinearsystem.jl b/test/nonlinearsystem.jl index fb4acb618a..19ebc12374 100644 --- a/test/nonlinearsystem.jl +++ b/test/nonlinearsystem.jl @@ -306,4 +306,13 @@ sys = structural_simplify(ns; conservative = true) 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 + 2*z + 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;;]) + 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 From d24e326297131d9f4e042038abfed59e90c3543f Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Tue, 16 Jul 2024 19:40:43 +0200 Subject: [PATCH 09/12] Handle case without observed equations --- src/systems/nonlinear/nonlinearsystem.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl index e1ad9959af..ba6ee21ef6 100644 --- a/src/systems/nonlinear/nonlinearsystem.jl +++ b/src/systems/nonlinear/nonlinearsystem.jl @@ -196,7 +196,7 @@ function calculate_jacobian(sys::NonlinearSystem; sparse = false, simplify = fal # observed equations may depend on unknowns, so substitute them in first # TODO: must do the same fix in e.g. calculate_hessian? # TODO: rather keep observed derivatives unexpanded, like "Differential(obs)(expr)"? - obs = map(eq -> eq.lhs => eq.rhs, observed(sys)) |> todict + 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)] From a1a17a6785d46a8a029637101246827b47b36e39 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Tue, 16 Jul 2024 19:45:53 +0200 Subject: [PATCH 10/12] Do the same for the Hessian --- src/systems/nonlinear/nonlinearsystem.jl | 3 ++- test/nonlinearsystem.jl | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl index ba6ee21ef6..f980ad5d5d 100644 --- a/src/systems/nonlinear/nonlinearsystem.jl +++ b/src/systems/nonlinear/nonlinearsystem.jl @@ -220,7 +220,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)] diff --git a/test/nonlinearsystem.jl b/test/nonlinearsystem.jl index 19ebc12374..11ebaaf32b 100644 --- a/test/nonlinearsystem.jl +++ b/test/nonlinearsystem.jl @@ -285,7 +285,7 @@ sys = structural_simplify(ns; conservative = true) @test length(equations(sys)) == 1 # https://github.com/SciML/ModelingToolkit.jl/issues/2858 -@testset "Jacobian with observed equations that depend on unknowns" begin +@testset "Jacobian/Hessian with observed equations that depend on unknowns" begin @variables x y z @parameters σ ρ β eqs = [0 ~ σ * (y - x) @@ -312,6 +312,7 @@ sys = structural_simplify(ns; conservative = true) eqs = [0 ~ x^2 + 2*z + 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 From ce3d5972850abaf907a5557c2fbd45d043ea92ed Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Tue, 16 Jul 2024 19:48:48 +0200 Subject: [PATCH 11/12] Format --- test/nonlinearsystem.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/nonlinearsystem.jl b/test/nonlinearsystem.jl index 11ebaaf32b..7d982ff7da 100644 --- a/test/nonlinearsystem.jl +++ b/test/nonlinearsystem.jl @@ -309,9 +309,9 @@ sys = structural_simplify(ns; conservative = true) # system that contains a chain of observed variables when simplified @variables x y z - eqs = [0 ~ x^2 + 2*z + y, z ~ y, y ~ x] # analytical solution x = y = z = 0 or -3 + 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(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()) From 3088c80794fd9b03c4f41faa17ff1fb720c82d32 Mon Sep 17 00:00:00 2001 From: Herman Sletmoen Date: Tue, 16 Jul 2024 19:50:33 +0200 Subject: [PATCH 12/12] Removed TODO --- src/systems/nonlinear/nonlinearsystem.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl index f980ad5d5d..43d4bc8cc5 100644 --- a/src/systems/nonlinear/nonlinearsystem.jl +++ b/src/systems/nonlinear/nonlinearsystem.jl @@ -194,7 +194,6 @@ function calculate_jacobian(sys::NonlinearSystem; sparse = false, simplify = fal end # observed equations may depend on unknowns, so substitute them in first - # TODO: must do the same fix in e.g. calculate_hessian? # 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))