Objective Sensitivity

src/NonLinearProgram/NonLinearProgram.jl

@@ -27,10 +27,12 @@ end
 Base.@kwdef struct ForwCache
     primal_Δs::Dict{MOI.VariableIndex,Float64}  # Sensitivity for primal variables (indexed by VariableIndex)
     dual_Δs::Vector{Float64}           # Sensitivity for constraints and bounds (indexed by ConstraintIndex)
+    dual_p::Dict{MOI.ConstraintIndex,Float64}  # Dual wrt parameters
 end
 
 Base.@kwdef struct ReverseCache
     Δp::Vector{Float64}  # Sensitivity for parameters
+    dual_p::Dict{MOI.ConstraintIndex,Float64}  # Dual wrt parameters
 end
 
 # Define the form of the NLP
@@ -513,15 +515,21 @@ function DiffOpt.forward_differentiate!(model::Model; tol = 1e-6)
         end
 
         # Compute Jacobian
-        Δs = _compute_sensitivity(model; tol = tol)
+        Δs, df_dp = _compute_sensitivity(model; tol = tol)
 
         # Extract primal and dual sensitivities
         primal_Δs = Δs[1:length(model.cache.primal_vars), :] * Δp # Exclude slacks
         dual_Δs = Δs[cache.index_duals, :] * Δp # Includes constraints and bounds
 
+        # Dual wrt parameters
+        varorder =
+            sort(collect(keys(form.var2ci)); by = x -> form.var2ci[x].value)
+        dual_p = [df_dp[form.var2param[var_idx].value] for var_idx in varorder]
+
         model.forw_grad_cache = ForwCache(;
             primal_Δs = Dict(model.cache.primal_vars .=> primal_Δs),
             dual_Δs = dual_Δs,
+            dual_p = Dict([form.var2ci[var_idx] for var_idx in varorder] .=> dual_p),
         )
     end
     return nothing
@@ -533,7 +541,7 @@ function DiffOpt.reverse_differentiate!(model::Model; tol = 1e-6)
         form = model.model
 
         # Compute Jacobian
-        Δs = _compute_sensitivity(model; tol = tol)
+        Δs, df_dp = _compute_sensitivity(model; tol = tol)
         num_primal = length(cache.primal_vars)
         # Fetch primal sensitivities
         Δx = zeros(num_primal)
@@ -576,7 +584,12 @@ function DiffOpt.reverse_differentiate!(model::Model; tol = 1e-6)
             sort(collect(keys(form.var2ci)); by = x -> form.var2ci[x].value)
         Δp = [Δp[form.var2param[var_idx].value] for var_idx in varorder]
 
-        model.back_grad_cache = ReverseCache(; Δp = Δp)
+        # Dual wrt parameters
+        dual_p = [df_dp[form.var2param[var_idx].value] for var_idx in varorder]
+
+        model.back_grad_cache = ReverseCache(; Δp = Δp, 
+            dual_p = Dict([form.var2ci[var_idx] for var_idx in varorder] .=> dual_p),
+        )
     end
     return nothing
 end
@@ -610,4 +623,16 @@ function MOI.get(
     return MOI.Parameter{T}(model.back_grad_cache.Δp[ci.value])
 end
 
+function MOI.get(
+    model::Model,
+    ::MOI.ConstraintDual,
+    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Parameter{T}},
+) where {T}
+    if !isnothing(model.forw_grad_cache)
+        return model.forw_grad_cache.dual_p[ci]
+    else
+        return model.back_grad_cache.dual_p[ci]
+    end
+end
+
 end # module NonLinearProgram

src/NonLinearProgram/nlp_utilities.jl

@@ -27,6 +27,13 @@ function _fill_off_diagonal(H::SparseMatrixCSC)
     return ret
 end
 
+function _compute_gradient(model::Model)
+    evaluator = model.cache.evaluator
+    grad = zeros(length(model.x))
+    MOI.eval_objective_gradient(evaluator, grad, model.x)
+    return grad
+end
+
 """
     _compute_optimal_hessian(evaluator::MOI.Nonlinear.Evaluator, rows::Vector{JuMP.ConstraintRef}, x::Vector{JuMP.VariableRef})
 
@@ -104,7 +111,7 @@ function _create_evaluator(form::Form)
         backend,
         MOI.VariableIndex.(1:form.num_variables),
     )
-    MOI.initialize(evaluator, [:Hess, :Jac])
+    MOI.initialize(evaluator, [:Hess, :Jac, :Grad])
     return evaluator
 end
 
@@ -496,5 +503,10 @@ function _compute_sensitivity(model::Model; tol = 1e-6)
     ∂s[num_w+num_cons+1:num_w+num_cons+num_lower, :] *= _sense_multiplier
     # Dual bounds upper
     ∂s[num_w+num_cons+num_lower+1:end, :] *= -_sense_multiplier
-    return ∂s
+
+    # dual wrt parameter
+    primal_idx = [i.value for i in model.cache.primal_vars]
+    df_dx = _compute_gradient(model)[primal_idx]
+    df_dp = df_dx'∂s[1:num_vars, :]
+    return ∂s, df_dp
 end

src/moi_wrapper.jl

@@ -732,6 +732,26 @@ function MOI.get(
     )
 end
 
+function MOI.get(
+    model::Optimizer,
+    attr::MOI.ConstraintDual,
+    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Parameter{T}},
+) where {T}
+    try
+        return MOI.get(
+            _checked_diff(model, attr, :forward_differentiate!),
+            attr,
+            model.index_map[ci],
+        )
+    catch
+        return MOI.get(
+            _checked_diff(model, attr, :reverse_differentiate!),
+            attr,
+            model.index_map[ci],
+        )
+    end
+end
+
 function MOI.supports(
     ::Optimizer,
     ::ReverseVariablePrimal,

test/nlp_program.jl

@@ -161,6 +161,8 @@ function test_analytical_simple(; P = 2) # Number of parameters
         # Compute derivatives
         DiffOpt.forward_differentiate!(m)
 
+        @test all(isapprox.(dual.(ParameterRef.(p)), dual.(con); atol = 1e-8))
+
         # Test sensitivities 
         @test_throws ErrorException MOI.get(
             m.moi_backend.optimizer.model.diff.model,
@@ -736,6 +738,8 @@ function test_ReverseConstraintDual()
     # Compute derivatives
     DiffOpt.reverse_differentiate!(m)
 
+    @test all(isapprox.(dual.(ParameterRef.(p)), dual.(con); atol = 1e-8))
+
     # Test sensitivities ReverseConstraintSet
     @test all(
         isapprox(