Update optimization model by changing variables to coefficients

Hi there, :slight_smile:
I am trying to extract shadow prices out of a MIQP model. In order to do this I want to “fix” the mixed integer variables after a first solving and then extract the duals.

Here’s a minimal working example:

using JuMP
using Gurobi

model = Model(Gurobi.Optimizer)

@variable(model, x, Bin)
@variable(model, y)

@constraint(model, x * y >= 1)
@constraint(model, c1, x + y <= 6)
@objective(model, Max, y)

optimize!(model)

solution = Dict(v => value(v) for v in all_variables(model))
for v in all_variables(model)
    is_binary(v) || continue
    unset_binary(v)
    fix(v, solution[v]; force=true)
end
new_model = copy(model)
set_optimizer(new_model, Gurobi.Optimizer)
set_optimizer_attribute(new_model, "QCPDual", 1)
optimize!(new_model)

JuMP.shadow_price(new_model[:c1])

However, for larger and more difficult models, QCPDual is not guaranteed to compute duals according to this post.

I was thinking if there is an option to get fully rid of the binary/integer variables and replace those by the result values in order to end up with a full LP. Is there a way to e.g. iterate through my constraints by using e.g. all_constraints(model; include_variable_in_set_constraints=true) replace the variables and afterwards delete the variables from the model?
In the end I would hope to get then a model in this case which looks like this:

max(y)
s.t.
1 * y >= 1
1 + y <= 6

So far I was not able to figure this out…

This is an interesting question. We don’t have a good way to solve this, unfortunately, and I can’t even recommend a simple work-around.

The fixed variable approach should be sufficient. Do you have an example where the QCP dual failed?

Also note that you can fix the variables with:

undo = fix_discrete_variables(model);

See Computing the duals of a mixed-integer program · JuMP

2 Likes

I don’t think this problem is well posed.

@max Dual variable (shadow prices) exists in convex conic programs, not in QP.
Because QP may be convex or nonconvex. Therefore, if you want a dual variable, you should first formulate your problem as a convex quadratic conic program, not only a QP.

If your intention is to get an LP with adjustable parameters, then just natively write a function in which you build this LP, and the function arguments being those parameters.

The problem is well posed, and Gurobi can compute duals for QPs (and QCPs, if you set QCPDual to 1).

Here’s how I would write the code in JuMP:

julia> using JuMP, Gurobi

julia> begin
           model = Model(Gurobi.Optimizer)
           set_silent(model)
           @variable(model, x, Bin)
           @variable(model, y)
           @constraint(model, x * y >= 1)
           @constraint(model, c1, x + y <= 6)
           @objective(model, Max, y)
           optimize!(model)
           undo = fix_discrete_variables(model)
           set_attribute(model, "QCPDual", 1)
           optimize!(model)
           @assert is_solved_and_feasible(model)
           shadow_price(c1)
       end
1.0
1 Like

I don’t think in a realistic operations research problem we will make use of the dual variable of a general QP (or QCQP).

We tend to write dual programs for convex conic programs (including convex SOCP). (Do we also write dual programs for QP? seems exists but I don’t encounter it often).

The Gurobi solver is usable but some times may have problems like Gurobi Error 10005: Unable to retrieve attribute 'Pi' - #7 by WalterMadelim

I have a bit of experience about the Gurobi solver. Sometimes I encounter strange behaviors. Perhaps tuning those parameters could help make it more stable, like Presolve, DualReductions

Do we also write dual programs for QP

See Duality · MathOptInterface

Thanks for your quick answers! Actually, I cannot provide an easy example where I encounter the failed computation of the Duals by QCPDual

However, I thought it could be an “easy” way of getting rid of the QP formulations (at least in my usecase). I only have QP formulations in my model where I multiply binaries with linear continuous variables. That’s why I thought I could remove those formulations.

Maybe I need to consider reformulating my Model

You can reformulate x * y to:

using JuMP
function bilinear(model, x, y)
    @assert is_binary(x)
    @assert has_lower_bound(y)
    @assert has_upper_bound(y)
    y_L, y_U = lower_bound(y), upper_bound(y)
    xy = @variable(model)
    @constraint(model, xy <= y_U * x)
    @constraint(model, xy >= y_U * x)
    @constraint(model, xy <= y - y_L * (1 - x))
    @constraint(model, xy >= y - y_U * (1 - x))
    return xy
end

model = Model()
@variable(model, x, Bin)
@variable(model, -100 <= y <= 100)
xy = bilinear(model, x, y)
@constraint(model, xy >= 1)
1 Like