Do the same for the Hessian

Author: hersle

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)]

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