Tinkering with values of variables in a solved JuMP model

Specific Domains Optimization (Mathematical)
1

I’d like to solve a JuMP model, use primal_feasibility_report to check which constraints are violated, tweak the variable values until the constraints are all strictly respected, and then compute the new objective value with the new variables. Is it possible to do so?

I can pass a Dict to primal_feasiblity_report with the new values. That’s a bit inconvenient as it seems I have to do scalar by scalar, I can’t pass a vector or a matrix at once, but it works. I see no way to recompute the objective value, though.

2

Does something like this work for you?

julia> using JuMP, HiGHS

julia> profit = [5, 4, 3, 7, 6];

julia> weight = [2, 3, 1, 4, 3];

julia> capacity = 8;

julia> N = length(profit);

julia> model = Model(HiGHS.Optimizer);

julia> set_silent(model);

julia> @variable(model, x[1:N], Bin);

julia> @objective(model, Max, profit' * x);

julia> @constraint(model, c, weight' * x <= capacity);

julia> optimize!(model);

julia> x_sol = value.(x)
5-element Vector{Float64}:
 0.0
 0.0
 1.0
 1.0
 1.0

julia> x_new = copy(x_sol);

julia> x_new[1] = 1.0;

julia> x_new[2] = 0.5;

julia> d = Dict(x .=> x_new)
Dict{VariableRef, Float64} with 5 entries:
  x[4] => 1.0
  x[1] => 1.0
  x[3] => 1.0
  x[2] => 0.5
  x[5] => 1.0

julia> f = v -> d[v]
#14 (generic function with 1 method)

julia> value(f, objective_function(model))
23.0

julia> value(f, c)
11.5

julia> primal_feasibility_report(f, model)
Dict{Any, Float64} with 2 entries:
  x[2] binary                                      => 0.5
  c : 2 x[1] + 3 x[2] + x[3] + 4 x[4] + 3 x[5] ≤ 8 => 3.5

If the Dict is annoying to construct in your case, maybe there are better ways to construct the variable to value function (in this case f)?

3

If the solver says it’s OPTIMAL or it has a suboptimal solution, then that solution is already deemed feasible.

It’s fairly easy and you don’t have to resort to primal_feasibility_report. You only need to fix the variable’s lower_bound and upper_bound to your candidate solution and call optimize! again, the solver will help you compute the objective value, and you fetch that value through objective_value(model).

But the bigger issue is that, just as you’ve used the word “tinker”—it may not be helpful to tweak the solution returned from the solver (What’s the purpose? You typically get inferior objective values).

4

That doesn’t work, solvers can’t impose constraints exactly, they only work up to some tolerance. What I need is to compute an objective value for a solution that respects the constraints exactly (up to numerical precision).

5

Amazing, that works! value(f, objective_function(model)) is what I needed. I also didn’t know about the broadcast syntax for building dictionaries, that’s very convenient, and it works even for symmetric matrices.

using JuMP, LinearAlgebra
import Hypatia

function feasibility_test()
    model = Model(Hypatia.Optimizer)
    @variable(model, x[1:3])
    @constraint(model, x in JuMP.SecondOrderCone())
    @variable(model, ρ[1:2,1:2] in JuMP.PSDCone())
    @constraint(model, x[1] ≤ 1)
    @constraint(model, tr(ρ) == 1)
    @objective(model, Max, x[2] + x[3] + ρ[1, 2])
    optimize!(model)
    display(objective_value(model))
    display(primal_feasibility_report(model))
    newρ = value(ρ)
    newx = value(x)
    newx[1] = 1
    newx[2:3] .-= 1e-5
    d = merge(Dict(ρ .=> newρ), Dict(x .=> newx))
    display(primal_feasibility_report(model, d))
    f = x -> d[x]
    display(value(f, objective_function(model)))
end

Unfortunately it doesn’t work with the case I actually need, complex variables, but that should be easy enough to implement: can't use primal_feasibility_report with complex variables · Issue #3113 · jump-dev/MathOptInterface.jl · GitHub

6

For the record, now the issue has been fixed, and there is proper support for complex variables:

using JuMP, LinearAlgebra
import Hypatia

function complex_test()
    model = Model(Hypatia.Optimizer)
    @variable(model, ρ[1:2,1:2] in HermitianPSDCone())
    @constraint(model, tr(ρ) == 1)
    @objective(model, Max, imag(ρ[1, 2]))
    set_silent(model)
    optimize!(model)
    display(primal_feasibility_report(model))
    newρ = value(ρ)
    newρ[2, 2] = 1 - newρ[1, 1]

    shape = JuMP.HermitianMatrixShape(2)
    vars = [convert(VariableRef, x) for x in vectorize(ρ, shape)]
    dict = Dict(vars .=> vectorize(newρ, shape))
    display(primal_feasibility_report(model, dict))
    f(x) = dict[x]
    return value(f, objective_function(model))
end
7

Ignore this if you knew already, but there is also a technique called ‘Modelling to Generate Alternatives’ (MGA) which seems to be close to what you are trying to do? This is commonly used in energy system optimization, because there are often many solutions that differe only a little bit in the final system cost. There is even a package for it: NearOptimalAlternatives.jl

8

I guess I need to understand why MGA isn’t just multi objective where the first objective is cost and the remaining objectives are different measures of the solution. In which case, we could use GitHub - jump-dev/MultiObjectiveAlgorithms.jl: A Julia package for solving multi-objective optimization problems · GitHub