diff --git a/lib/OptimizationNLPModels/Project.toml b/lib/OptimizationNLPModels/Project.toml
index f72e628b4..24cf6c9de 100644
--- a/lib/OptimizationNLPModels/Project.toml
+++ b/lib/OptimizationNLPModels/Project.toml
@@ -8,22 +8,24 @@ ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
 ADTypes = "1.7"
 NLPModels = "0.21"
 Optimization = "4"
 Reexport = "1.2"
+SparseArrays = "1.11.0"
 julia = "1.9"
 
 [extras]
+Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 NLPModelsTest = "7998695d-6960-4d3a-85c4-e1bceb8cd856"
+OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
-Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
-OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
 
 [targets]
 test = ["Test", "NLPModelsTest", "OptimizationOptimJL", "ReverseDiff", "Zygote", "Ipopt", "OptimizationMOI"]
diff --git a/lib/OptimizationNLPModels/src/OptimizationNLPModels.jl b/lib/OptimizationNLPModels/src/OptimizationNLPModels.jl
index 65e67be20..cf9d5166d 100644
--- a/lib/OptimizationNLPModels/src/OptimizationNLPModels.jl
+++ b/lib/OptimizationNLPModels/src/OptimizationNLPModels.jl
@@ -1,6 +1,6 @@
 module OptimizationNLPModels
 
-using Reexport
+using Reexport, SparseArrays
 @reexport using NLPModels, Optimization, ADTypes
 
 """
@@ -21,9 +21,20 @@ function SciMLBase.OptimizationFunction(nlpmodel::AbstractNLPModel,
         cons(res, x, p) = NLPModels.cons!(nlpmodel, x, res)
         cons_j(J, x, p) = (J .= NLPModels.jac(nlpmodel, x))
         cons_jvp(Jv, v, x, p) = NLPModels.jprod!(nlpmodel, x, v, Jv)
+        function lag_h(h, θ, σ, λ)
+            H = NLPModels.hess(nlpmodel, θ, λ; obj_weight = σ)
+            k = 0
+            rows, cols, _ = findnz(H)
+            for (i, j) in zip(rows, cols)
+                if i <= j
+                    k += 1
+                    h[k] = H[i, j]
+                end
+            end
+        end
 
         return OptimizationFunction(
-            f, adtype; grad, hess, hv, cons, cons_j, cons_jvp, kwargs...)
+            f, adtype; grad, hess, hv, cons, cons_j, cons_jvp, lag_h, kwargs...)
     end
 
     return OptimizationFunction(f, adtype; grad, hess, hv, kwargs...)