fix hessian calculations

Author: ChrisRackauckas

src/systems/optimization/optimizationsystem.jl

@@ -63,11 +63,12 @@ end
 
 function generate_hessian(sys::OptimizationSystem, vs = states(sys), ps = parameters(sys);
                           sparse = false, kwargs...)
-    hes = calculate_hessian(sys)
     if sparse
-        hes = sparse(hes)
+        hess = sparsehessian(equations(sys),[dv() for dv in states(sys)])
+    else
+        hess = calculate_hessian(sys)
     end
-    return build_function(hes, convert.(Variable,vs), convert.(Variable,ps);
+    return build_function(hess, convert.(Variable,vs), convert.(Variable,ps);
                           conv = AbstractSysToExpr(sys),kwargs...)
 end
 
@@ -89,6 +90,10 @@ struct AutoModelingToolkit <: DiffEqBase.AbstractADType end
 DiffEqBase.OptimizationProblem(sys::OptimizationSystem,args...;kwargs...) =
     DiffEqBase.OptimizationProblem{true}(sys::OptimizationSystem,args...;kwargs...)
 
+OptimizationProblemExpr(sys::OptimizationSystem,args...;kwargs...) =
+    OptimizationProblemExpr{true}(sys::OptimizationSystem,args...;kwargs...)
+
+
 """
 ```julia
 function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem,
@@ -166,8 +171,8 @@ struct OptimizationProblemExpr{iip} end
 function OptimizationProblemExpr{iip}(sys::OptimizationSystem, u0,
                                           parammap=DiffEqBase.NullParameters();
                                           lb=nothing, ub=nothing,
-                                          grad = true,
-                                          hes = false, sparse = false,
+                                          grad = false,
+                                          hess = false, sparse = false,
                                           checkbounds = false,
                                           linenumbers = false, parallel=SerialForm(),
                                           kwargs...) where iip

test/modelingtoolkitize.jl

@@ -48,13 +48,14 @@ x0 = zeros(2)
 p  = [1.0,100.0]
 
 prob = OptimizationProblem(rosenbrock,x0,p)
-sys = modelingtoolkitize(prob)
-x0map = states(sys) .=> x0
-parammap = parameters(sys) .=> p
+sys = modelingtoolkitize(prob) # symbolicitize me captain!
 
-prob = OptimizationProblem(sys,x0map,parammap,grad=true)
+prob = OptimizationProblem(sys,x0,p,grad=true,hess=true)
 sol = solve(prob,NelderMead())
 @test sol.minimum < 1e-8
 
 sol = solve(prob,BFGS())
 @test sol.minimum < 1e-8
+
+sol = solve(prob,Newton())
+@test sol.minimum < 1e-8