Update optimization model by changing variables to coefficients

max · January 10, 2025, 1:48pm

Hi there,
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…

odow · January 10, 2025, 7:35pm

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

WalterMadelim · January 11, 2025, 1:11am

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.

odow · January 11, 2025, 1:17am

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

WalterMadelim · January 11, 2025, 2:25am

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

odow · January 11, 2025, 3:35am

Do we also write dual programs for QP

See Duality · MathOptInterface

max · January 13, 2025, 11:33am

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

odow · January 13, 2025, 8:25pm

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)

WalterMadelim · March 24, 2025, 3:01am

In your example x * y >= 1 happens to be a convex constraint.
If I change the sign, then it ends up with Error, like this

julia> using JuMP, Gurobi

julia> begin
           model = Model(Gurobi.Optimizer)
               @variable(model, 0 <= x <= 1)
           @variable(model, y)
           @constraint(model, x * y <= 1)
           @constraint(model, c1, x + y <= 6)
           @objective(model, Max, y)
           set_attribute(model, "QCPDual", 1)
       end
Set parameter Username
Set parameter LicenseID to value 2602363
Academic license - for non-commercial use only - expires 2025-12-20
Set parameter QCPDual to value 1

julia> optimize!(model)
Set parameter QCPDual to value 1
Gurobi Optimizer version 12.0.1 build v12.0.1rc0 (win64 - Windows 11.0 (26100.2))

CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Non-default parameters:
QCPDual  1

Optimize a model with 1 rows, 2 columns and 2 nonzeros
Model fingerprint: 0x38be9e11
Model has 1 quadratic constraint
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  QMatrix range    [1e+00, 1e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [6e+00, 6e+00]
  QRHS range       [1e+00, 1e+00]

Error: Continuous model is non-convex but QCP duals are requested.
       Either set QCPDuals to 0, or set NonConvex to 2 to solve model as a MIP.


User-callback calls 31, time in user-callback 0.00 sec
ERROR: Gurobi Error 10020: Constraint Q not PSD (diagonal adjustment of 1.0e+00 would be required). Set NonConvex parameter to -1 or 2 to solve model.
Stacktrace:
 [1] _check_ret
   @ K:\judepot1114\packages\Gurobi\QK0tn\src\MOI_wrapper\MOI_wrapper.jl:417 [inlined]
 [2] _check_ret_GRBoptimize(model::Gurobi.Optimizer)
   @ Gurobi K:\judepot1114\packages\Gurobi\QK0tn\src\MOI_wrapper\MOI_wrapper.jl:434
 [3] optimize!(model::Gurobi.Optimizer)
   @ Gurobi K:\judepot1114\packages\Gurobi\QK0tn\src\MOI_wrapper\MOI_wrapper.jl:2799
 [4] optimize!
   @ K:\judepot1114\packages\MathOptInterface\TYq6d\src\Bridges\bridge_optimizer.jl:367 [inlined]
 [5] optimize!
   @ K:\judepot1114\packages\MathOptInterface\TYq6d\src\MathOptInterface.jl:122 [inlined]
 [6] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{…})
   @ MathOptInterface.Utilities K:\judepot1114\packages\MathOptInterface\TYq6d\src\Utilities\cachingoptimizer.jl:327
 [7] optimize!(model::Model; ignore_optimize_hook::Bool, _differentiation_backend::MathOptInterface.Nonlinear.SparseReverseMode, kwargs::@Kwargs{})
   @ JuMP K:\judepot1114\packages\JuMP\xlp0s\src\optimizer_interface.jl:595
 [8] optimize!(model::Model)
   @ JuMP K:\judepot1114\packages\JuMP\xlp0s\src\optimizer_interface.jl:546
 [9] top-level scope
   @ REPL[3]:1
Some type information was truncated. Use `show(err)` to see complete types.

WalterMadelim · March 24, 2025, 3:07am

If we follow the Error message, then

set_attribute(model, "NonConvex", 2)
optimize!(model)
JuMP.dual_status(model) # NO_SOLUTION::ResultStatusCode = 0
JuMP.shadow_price(c1)
ERROR: The shadow price is not available because no dual result is available.

odow · March 24, 2025, 3:11am

Note that in the original example x is binary. It is then fixed to a value. So your posts and errors are because you have relaxed the bounds on x, not because of the sign of x * y <= 1 (which is always non-convex, regardless of >= or <=).

WalterMadelim · March 24, 2025, 4:52am

I took a closer look. Here is your original code (which works)

using JuMP, Gurobi
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) # 🍅 This is an MIQCP
    undo = fix_discrete_variables(model)
    # 🍅 At this line, `model` is continuous, convex and quadratic 
    # Therefore, we can proceed with
    set_attribute(model, "QCPDual", 1)
    optimize!(model)
    assert_is_solved_and_feasible(model; dual = true)
    shadow_price(c1)
end

Actually I retained the bounds via 0 <= x <= 1. I relaxed the integerality constraint of x to let it become a continuous QCP, which is easier than it was.

And about the bilinear constraint x * y >= 1, it is a bit subtle.

This constraint is itself nonconvex
However, if we restrict x >= 0 additionally (which is then equivalent to x > 0), then Gurobi can identify it as convex.

Here are 2 related examples who are both convex (and Gurobi can identify)
The first is

begin
    model = Model(Gurobi.Optimizer)
    @variable(model, x >= 0)
    @variable(model, y)
    @constraint(model, x * y >= 1)
    @constraint(model, c1, x + y <= 6)
    @objective(model, Max, y)
    set_attribute(model, "QCPDual", 1)
end
optimize!(model) # convex
assert_is_solved_and_feasible(model; dual = true)

The second is

begin
    model = Model(Gurobi.Optimizer)
    @variable(model, x <= 0)
    @variable(model, y)
    @constraint(model, x * y >= 1)
    @constraint(model, c1, x + y <= 6)
    @objective(model, Max, y)
    set_attribute(model, "QCPDual", 1)
end
optimize!(model) # convex
assert_is_solved_and_feasible(model; dual = true)

An exceptional point for me is that Gurobi can Identify this program as convex quadratic

    model = Model(Gurobi.Optimizer)
    @variable(model, y)
    @variable(model, x)
    set_lower_bound(x, 1)
    set_upper_bound(x, 1)
    @constraint(model, x * y <= 1)
    @objective(model, Max, y)
    optimize!(model)
    JuMP.solution_summary(model; verbose = true)

Although it still fails to directly identify it as a simpler linear program

    @variable(model, y)
    @constraint(model, y <= 1)
    @objective(model, Max, y)

Another exceptional point for me is that Gurobi can even provide the dual variable of non-affine constraint like x * y >= 1, although I’m unsure if it has a usage in reality.

Topic		Replies	Views
Dual vs shadow price: Which one to use? Optimization (Mathematical) question , jump	8	1584	May 23, 2023
Relaxation of MIP in order to use getdual Optimization (Mathematical) question , jump	6	1193	December 4, 2021
JuMP Support for solver specific methods to create fixed MIP Optimization (Mathematical)	5	246	February 27, 2023
Creating fixed model associated with a MIP model using Gurobi in JuMP Optimization (Mathematical) jump	1	591	October 19, 2022
[JuMP + Gurobi Bug?] value of variable in solved model changes after deleting another variable Optimization (Mathematical) jump , bug	7	1147	October 8, 2019

Update optimization model by changing variables to coefficients

Related topics