Functions like set_binary, etc. does not work in a user-defined function

Specific Domains Optimization (Mathematical)
1

Hello,
I created a simple function that solves a JuMP MIP model with arrays as inputs.

The problem I want to solve is

minimize c’*x
s.t.
Ax = b
l <= x <= u
some x[i] in {0, 1}

[vartype is a vector that defines the corresponding variable type of x (:Cont or :Bin)]

function milprog(c, A, sense, b, l, u, vartype, solver)
    N = length(c)
    model = Model(solver)
    @variable(model, l[i] <= x[i=1:N] <= u[i])
    for i in eachindex(vartype)
        if vartype[i] == :Bin
            set_binary(x[i])
        end
        println(i, is_binary(x[i]))
    end
    @objective(model, Min, c' * x)
    eq_rows, ge_rows, le_rows = sense .== '=', sense .== '>', sense .== '<'
    @constraint(model, A[eq_rows, :] * x .== b[eq_rows])
    @constraint(model, A[ge_rows, :] * x .>= b[ge_rows])
    @constraint(model, A[le_rows, :] * x .<= b[le_rows])
    optimize!(model)
    return (
        status = termination_status(model),
        objval = objective_value(model),
        sol = value.(x)
    )
end

However, calling milprog just solved LP problem instead of MILP, even though I verified that some x[i] was changed to binary by using is_binary. Anything wrong in my code?

Thanks for your help.

2

Can you provide a reproducible example? Which solver?

For example, this works:

julia> using JuMP, HiGHS

julia> function milprog(c, A, sense, b, l, u, vartype, solver)
           m, n = size(A)
           @assert length(sense) == length(b) == m
           @assert length(c) == length(l) == length(u) == length(vartype) == n
           model = Model(solver)
           @variable(model, l[i] <= x[i in 1:n] <= u[i], binary = (vartype[i] == :Bin))
           @objective(model, Min, c' * x)
           for i in 1:m
               if sense[i] == '='
                   @constraint(model, A[i, :]' * x == b[i])
               elseif sense[i] == '>'
                   @constraint(model, A[i, :]' * x >= b[i])
               else
                   @assert sense[i] == '<'
                   @constraint(model, A[i, :]' * x <= b[i])
               end
           end
           optimize!(model)
           if !is_solved_and_feasible(model)
               return (; status = termination_status(model),)
           end
           return (;
               status = termination_status(model),
               objval = objective_value(model),
               sol = value.(x),
           )
       end
milprog (generic function with 1 method)

julia> milprog(
           [12, 20],
           [6 8; 7 12],
           ['>', '>'],
           [100, 120],
           [0, 0],
           [Inf, 3],
           [:Cont, :Bin],
           HiGHS.Optimizer
       )
Running HiGHS 1.8.1 (git hash: 4a7f24ac6): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [6e+00, 1e+01]
  Cost   [1e+01, 2e+01]
  Bound  [1e+00, 1e+00]
  RHS    [1e+02, 1e+02]
Presolving model
2 rows, 2 cols, 4 nonzeros  0s
2 rows, 2 cols, 4 nonzeros  0s

Solving MIP model with:
   2 rows
   2 cols (1 binary, 0 integer, 0 implied int., 1 continuous)
   4 nonzeros
MIP-Timing:      0.0024 - starting analytic centre calculation

Src: B => Branching; C => Central rounding; F => Feasibility pump; H => Heuristic; L => Sub-MIP;
     P => Empty MIP; R => Randomized rounding; S => Solve LP; T => Evaluate node; U => Unbounded;
     z => Trivial zero; l => Trivial lower; u => Trivial upper; p => Trivial point

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
Src  Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   185.1428571     inf                  inf        0      0      0         0     0.0s
 T       0       0         0   0.00%   185.1428571     205.1428571        9.75%        0      0      0         1     0.0s
         1       0         1 100.00%   205.1428571     205.1428571        0.00%        0      0      0         1     0.0s

Solving report
  Status            Optimal
  Primal bound      205.142857143
  Dual bound        205.142857143
  Gap               0% (tolerance: 0.01%)
  P-D integral      6.90514725065e-07
  Solution status   feasible
                    205.142857143 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            0.00 (total)
                    0.00 (presolve)
                    0.00 (solve)
                    0.00 (postsolve)
  Max sub-MIP depth 0
  Nodes             1
  Repair LPs        0 (0 feasible; 0 iterations)
  LP iterations     1 (total)
                    0 (strong br.)
                    0 (separation)
                    0 (heuristics)
(status = MathOptInterface.OPTIMAL, objval = 205.14285714285714, sol = [15.428571428571429, 1.0])
1 Like
3

Thanks for your help.

The problem seems to be due to my bad. I tried to replicate in another machine and everything looks fine. I will look into what happend in yesterday…

1 Like
4

I assume it is probably that you had previously tried a different version of the milprog function in the REPL, and it was more specific (that is, it had a type annotation and your example above does not). Therefore, even though you thought you were calling a new model it was still calling the old one.

For example:

julia> function foo(x::Int)
           return 2 * x
       end
foo (generic function with 1 method)

julia> function foo(x)
           return -x
       end
foo (generic function with 2 methods)

julia> foo(1)
2

The second line doesn’t replace the existing function foo, it adds a new method.

A good trick when encountering unexpected results in Julia is to start from a fresh session. You’ll typically find that the bug was because some earlier method or variable was being used unexpectedly :smile:

1 Like