How to get the gradient of a constraint of an optimization problem at a fixed point

VPBML · October 13, 2024, 3:00pm

Hi,
I wonder how I can get the value of the gradient of a constraint at some fixed point using JuMP. For example, in the following optimization problem, how can I get the gradient of the constraint with name Temp with respect to variable u[1] at the starting point (provided via the start values)?
I would really appreciate it if you could help me with this problem.

using JuMP, Ipopt

model = Model(Ipopt)
@variable(model, Tc_out[1:2])
@variable(model, T)
@variable(model, Th_out[1:2])
@variable(model, 0<= u[1:2] <= 1, start = 0.5)
    
@NLparameter(model, Tc_in == 41.0)
@NLparameter(model, wc == 97.0)
@NLparameter(model, Th_in[1:2] == 105.0)
set_value(Th_in[2], 206.0)
@NLparameter(model, wh[1:2] == 118.0)
set_value(wh[2], 94.0)
@NLparameter(model, UA[1:2] == 240.0)
set_value(UA[2], 502.0)

@NLobjective(model, Max, 1)
@NLconstraint(model, Tc_outL[i = 1:2], Tc_out[i] >= Tc_in)
@NLconstraint(model, Tc_outH[i = 1:2], Tc_out[i] <= Th_in[i])
@NLconstraint(model, Th_outL[i = 1:2], Th_out[i] >= Tc_in)
@NLconstraint(model, Th_outH[i = 1:2], Th_out[i] <= Th_in[i])
@NLconstraint(model, HTrans[i = 1:2], -1*Tc_out[i] + Tc_in + UA[i]*((Th_in[i]-Tc_out[i])*(Th_out[i]-Tc_in)*0.5*(Th_in[i]-Tc_out[i]+Th_out[i]-Tc_in))^(1/3)/(u[i]*wc*4.2) == 0)
@NLconstraint(model, ECon_cold, -1*T + Tc_out[1]*u[1] + Tc_out[2]*u[2] == 0)
@NLconstraint(model, ECon_HX[i = 1:2], (-1*Th_out[i] + Th_in[i])*wh[i] - (Tc_out[i] - Tc_in)*wc*u[i] == 0)
@NLconstraint(model, MCon, u[1]+u[2] -1 == 0)
@NLconstraint(model, Temp, ((Tc_out[2] - Tc_in)^2)/(Th_in[2] - Tc_in) -  ((Tc_out[1] - Tc_in)^2)/(Th_in[1] - Tc_in) == 0)

set_start_value(model[:Tc_out][1], (41 + 105)*0.5)
set_start_value(model[:Tc_out][2], (41 + 206)*0.5)
set_start_value(model[:T], (2*41 + 105 + 206)*0.25)
set_start_value(model[:Th_out][1], (41 + 105)*0.5)
set_start_value(model[:Th_out][2], (41 + 206)*0.5)

dpo · October 13, 2024, 5:01pm

There may be a way with JuMP alone, but you can wrap your model in a MathProgNLPModel from NLPModelsJuMP.jl and evaluate any problem function or its derivatives at any point.

odow · October 13, 2024, 7:17pm

It’s a little but convoluted, but see Computing Hessians · JuMP.

Here’s how I would write and solve your model.

Note that:

I have swapped from the @NL syntax to the new JuMP nonlinear syntax (see ANN: JuMP v1.15 is released)
I used MOI.eval_constraint_gradient instead of _jacobian because you care only about one constraint. You could use the jacobian if you wanted
The gradient of the function in constraint Temp is 0 with respect to u[1] because u[1] does not appear in it. Is that what you expected?

julia> using JuMP

julia> import Ipopt

julia> import SparseArrays

julia> model = Model(Ipopt.Optimizer)
A JuMP Model
├ solver: Ipopt
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none

julia> @variables(model, begin
           Tc_out[1:2], (start = (41 + 105)*0.5)
           T, (start = 0.25 * (2*41 + 105 + 206))
           Th_out[1:2], (start = (41 + 105)*0.5)
           0 <= u[1:2] <= 1, (start = 0.5)
           Tc_in in Parameter(41.0)
           wc in Parameter(97.0)
           Th_in[1:2] in Parameter(105.0)
           wh[1:2] in Parameter(118.0)
           UA[1:2] in Parameter(240.0)
       end);

julia> set_parameter_value(Th_in[2], 206.0)

julia> set_parameter_value(wh[2], 94.0)

julia> set_parameter_value(UA[2], 502.0)

julia> @constraints(model, begin
           Tc_outL[i in 1:2], Tc_out[i] >= Tc_in
           Tc_outH[i in 1:2], Tc_out[i] <= Th_in[i]
           Th_outL[i in 1:2], Th_out[i] >= Tc_in
           Th_outH[i in 1:2], Th_out[i] <= Th_in[i]
           HTrans[i in 1:2],
               Tc_in - Tc_out[i] + UA[i] * (
                   (Th_in[i] - Tc_out[i]) * (Th_out[i] - Tc_in) * 0.5 * (Th_in[i] - Tc_out[i] + Th_out[i] - Tc_in)
               )^(1 / 3) / (u[i] * wc * 4.2) == 0
           ECon_cold, -T + Tc_out[1] * u[1] + Tc_out[2] * u[2] == 0
           ECon_HX[i in 1:2],
               (Th_in[i] - Th_out[i]) * wh[i] - (Tc_out[i] - Tc_in) * wc * u[i] == 0
           MCon, u[1] + u[2] == 1
           Temp,
               (Tc_out[2] - Tc_in)^2 / (Th_in[2] - Tc_in) -  (Tc_out[1] - Tc_in)^2 / (Th_in[1] - Tc_in) == 0
       end);

julia> function df_dx(model::Model, ci::ConstraintRef, x::VariableRef)
           function my_value(x::VariableRef)
               if is_parameter(x)
                   return parameter_value(x)
               else
                   return start_value(x)
               end
           end
           nlp = MOI.Nonlinear.Model()
           # Add ci as a constraint
           object = constraint_object(ci)
           MOI.Nonlinear.add_constraint(nlp, object.func, object.set)
           # Build the evaluator
           vars = all_variables(model)
           backend = MOI.Nonlinear.SparseReverseMode()
           evaluator = MOI.Nonlinear.Evaluator(nlp, backend, index.(vars))
           # Initialize
           MOI.initialize(evaluator, [:Jac])
           # Initialize the output structure
           sparsity = MOI.constraint_gradient_structure(evaluator, 1)
           ∇g = zeros(length(sparsity))
           # Compute gradient at `my_value(x)`
           MOI.eval_constraint_gradient(evaluator, ∇g, my_value.(vars), 1)
           # Convert to sparsevec
           J = SparseArrays.sparsevec(sparsity, ∇g, length(vars))
           # Return x element
           i = findfirst(==(x), vars)
           return J[i]
       end
df_dx (generic function with 1 method)

julia> df_dx(model, Temp, u[1])
0.0

VPBML · October 14, 2024, 8:37am

Hi @odow,
Thanks for the answer. I actually expected a non-zero gradient. I have tried your suggestion on a simpler example shown below:

using JuMP
import Ipopt
import SparseArrays

model = Model(Ipopt.Optimizer)
@variables(model, begin
u, (start = 10.0)
y1, (start = 30.0)
y2, (start = 30.0)
d in Parameter(10.0)
end)

@constraints(model, begin
def_y1, y1 - 2*u - d == 0
def_y2, y2 - u - 2*d == 0
def_SOC, y1 - y2 == 0
end)

function df_dx(model::Model, ci::ConstraintRef, x::VariableRef)
    function my_value(x::VariableRef)
        if is_parameter(x)
            return parameter_value(x)
        else
            return start_value(x)
        end
    end
    nlp = MOI.Nonlinear.Model()
    # Add ci as a constraint
    object = constraint_object(ci)
    MOI.Nonlinear.add_constraint(nlp, object.func, object.set)
    # Build the evaluator
    vars = all_variables(model)
    backend = MOI.Nonlinear.SparseReverseMode()
    evaluator = MOI.Nonlinear.Evaluator(nlp, backend, index.(vars))
    # Initialize
    MOI.initialize(evaluator, [:Jac])
    # Initialize the output structure
    sparsity = MOI.constraint_gradient_structure(evaluator, 1)
    ∇g = zeros(length(sparsity))
    # Compute gradient at `my_value(x)`
    MOI.eval_constraint_gradient(evaluator, ∇g, my_value.(vars), 1)
    # Convert to sparsevec
    J = SparseArrays.sparsevec(sparsity, ∇g, length(vars))
    # Return x element
    i = findfirst(==(x), vars)
    return J[i]
end

df_dx(model, def_SOC, u)

The result is 0 but we can see that, even though u is not in def_SOC:
def_SOC = y1 - y2 = 2*u + d - (u + 2*d) = u - d
So I want something that gives me df_dx(model, def_SOC, u) = 1 here.

odow · October 15, 2024, 12:51am

JuMP might not be the right tool for this. We cannot compute the derivative of a constraint’s function with respect to a variable as it is constrained by other constraints.

Why do you need this? What are you trying to achieve?

odow · October 15, 2024, 1:03am

Note that if you have a feasible primal-dual solution, then you can fix u to a value, set the function as the objective, and compute the reduced cost:

julia> using JuMP

julia> import Ipopt

julia> begin
           model = Model(Ipopt.Optimizer)
           set_silent(model)
           @variables(model, begin
               u == 10.0
               y1, (start = 30.0)
               y2, (start = 30.0)
               d in Parameter(10.0)
           end)
           @constraint(model, y1 == 2 * u - d)
           @constraint(model, y2 == u - 2 * d)
           @objective(model, Min, y1 - y2)
           optimize!(model)
           @assert is_solved_and_feasible(model)
           reduced_cost(u)
       end
1.0

VPBML · October 20, 2024, 8:22pm

Hi,
The reason why I need to evaluate the gradient of a constraint as such is quite bizarre, but basically I am trying to look for a constraint that satisfies a few criteria, one of which is about the gradient at some fixed points. I have a workaround for this problem of evaluating the gradient though.Thanks for your answers and suggestions. I really appreciate it!