Using CPLEX (solver dependent) callbacks for custom branching strategy in Julia

valdemarsoegaard · May 18, 2023, 12:08pm

Greetings!

I am working on my master thesis on a large-scale flow problem, where I have implemented a surrogate model derived from a neural network in the form of a set of mixed integer linear constraints which replaces a complicating non-convex, non-linear constraint. Luckily, that part is working fine, but now I want to implement my own custom branching strategy to test if I can speed up the solution time.

I have looked at my options, and it seems that CPLEX is one of (if not the only) commercial solver that will allow me to make my own custom branching strategy. To my understanding, what I am looking to use is the CPX_CALLBACKCONTEXT_BRANCHING callback, but I do not really understand how it works especially from a Julia/JuMP perspective, as they only have examples/pseudocode for the C api.
I am looking to do some simple calculations on some continuous variables, and based on that outcome, choose an associated fractional binary to branch on.

I have previously used Gurobi callbacks with solver-independent MOI callbacks (I fed a heuristic solution to the MIPnode) and that was very simple to implement. Looking through source code and examples online for the CPLEX solver-dependent callbacks though, it seems a bit more complicated. At least to my understanding, it does not fit into the categories of either ‘LazyConstraints’ ‘HeuristicSolution’ or ‘UserCuts’, so I cannot use the syntax for any of these e.g. MOI.submit(model, MOI.LazyConstraint(cb_data), con).

Anyway, while I have decent experience with coding and programming, mostly through Python and a bit of Julia, I have a hard time trying figure out exactly what kind of syntax etc. I am looking for.

Hope someone can help! If you need more information please let me know and I’ll be happy to provide you.

TL;DR I hope someone can help, maybe with a dummy example, on how to use the solver-dependent callbacks in CPLEX in a Julia context.

UPDATE:

I tried to look a bit more into the CPLEX documentation and the following example which seems to do something more like what i am looking for, but not exactly. I think I need to call the following function status = CPXcallbackmakebranch(context, varcnt, varind, varlu, varbd, rcnt, nzcnt, rhs, sense, rmatbeg, rmatind, rmatval, nodeest, seqnum_p) to create the desired branch after i have defined the desired fractional binary to branch on. Exactly how to provide the arguments correctly, i am still quite confused about though… I found this IBM Documentation that describes the parameters needed, but i don’t understand why i need to provide constraints, I only want to branch on a single variable, so in my mind, only the new upper and lower bound for the to child nodes are relevant, i guess there is obviously something about the function i just don’t understand

For context, here is my callback function so far:

 function my_callback_function(cb_data::CPLEX.CallbackContext, context_id::Clong)
    if context_id == CPX_CALLBACKCONTEXT_BRANCHING
        @info "callback called"

        # Before querying `callback_value`, you must call:
        CPLEX.load_callback_variable_primal(cb_data, context_id)

        # get the relevant values from the callback
        y_val = callback_value.(cb_data, m[:y])
        # y_hat_val = callback_value(cb_data, y_hat)
        xp_val = callback_value.(cb_data, m[:xp])
        xn_val =  callback_value.(cb_data, m[:xn])
        b_val = callback_value.(cb_data, m[:b])

        #Compute the maximum distance to the ReLU hyperplane
        d = (.-xn_val.+xp_val.+1)./(y_val.+1)
        #get the index of the associated fractional binary (if the ReLU is not enforced it will be > 1 and otherwise = 1)
        b_chs = argmax(d)
        @info "b_chs: $(b_chs)"
        #now create the child nodes (somehow)
  
        #Somehow add the appropriate arguments to the function call... 
        # status = CPXcallbackmakebranch(context, varcnt, varind, varlu, varbd, rcnt, nzcnt, rhs, sense, rmatbeg, rmatind, rmatval, nodeest, seqnum_p)
        #I only know that the context needs to be cb_data, and i think that the varind should be b_chs, but I don't know the rest of the arguments...
        status = CPXcallbackmakebranch(cb_data, ???, b_chs, ???, ???, ???, ???, ???, ???, ???, ???, ???, ???, ???)
    else
        return
    end
end

Best,
Valdemar

odow · May 18, 2023, 8:53pm

it seems that CPLEX is one of (if not the only) commercial solver that will allow me to make my own custom branching strategy

Correct.

as they only have examples/pseudocode for the C api.

“They” being JuMP or the IBM docs? The Julia-C API is essentially identical. The only difference is that you need to convert from JuMP variables to the C integer column using CPLEX.column(cb_data, x). And something like int* seqnum_p is seqnum_p = Ref{Cint}()

but i don’t understand why i need to provide constraints

You don’t. These are optional.

I have a hard time trying figure out exactly what kind of syntax etc. I am looking for.

Yeah. Using the solver-dependent callbacks requires some familiarity with C.

maybe with a dummy example, on how to use the solver-dependent callbacks in CPLEX in a Julia context.

Here you go. Does this help? I haven’t tested, so there might be bugs and typos. But it should point you in the right direction.

function my_callback_function(cb_data::CPLEX.CallbackContext, context_id::Clong)
    if context_id != CPX_CALLBACKCONTEXT_BRANCHING
      return
    end
    # Let's assume I cant to set the upper bound of m[:x] to 1.0
    column = CPLEX.column(cb_data, index(m[:x]))
    seqnum_p = Ref{Cint}(0)
    status = CPXcallbackmakebranch(
        cb_data, 
        1,              # varcnt
        Cint[column],   # varind
        Cchar['U'],     # varlu
        Cdouble[1.0],   # varbd
        0,              # rcnt
        0,              # nzcnt
        Cdouble[],      # rhs
        Cchar[],        # sense
        Cint[],         # rmatbeg
        Cint[],         # rmatind
        Cdouble[],      # rmatval
        0.0,            # nodeest  TODO: what is an objective estimate for the branch?
        seqnum_p,       # seqnum_p
    )
end

valdemarsoegaard · May 19, 2023, 11:23am

Thank you so much for the help! I have not really dealt with anything related to C before, so I would never have guessed how to use this syntax on my own.

So, it runs now, but I am a bit in doubt about the output I get.

My function now looks the following way:

 function my_callback_function(cb_data::CPLEX.CallbackContext, context_id::Clong)
    if context_id != CPX_CALLBACKCONTEXT_BRANCHING

        @info "callback called"

        # Before querying `callback_value`, you must call:
        CPLEX.load_callback_variable_primal(cb_data, context_id)

        # get the relevant values from the callback
        y_val = callback_value.(cb_data, m[:y])
        xp_val = callback_value.(cb_data, m[:xp])
        xn_val =  callback_value.(cb_data, m[:xn])
        b_val = callback_value.(cb_data, m[:b])

        #Compute the maximum distance to the ReLU hyperplane
        d = (.-xn_val.+xp_val.+1)./(y_val.+1) #add the '+1' to avoid division by zero
        #get the index of the associated fractional binary (if the ReLU is not enforced it will be > 1 and otherwise = 1)
        b_chs = argmax(d)

        #load the variable i want to branch on with proper CPLEX syntax:
        column = CPLEX.column(cb_data, index(m[:b][b_chs]))
  
        status = CPXcallbackmakebranch(
            cb_data, 
            1,              # varcnt
            Cint[column],   # varind
            Cchar['U'],     # varlu
            Cdouble[1.0],   # varbd #need to make two child nodes... 
            0,              # rcnt
            0,              # nzcnt
            Cdouble[],      # rhs
            Cchar[],        # sense
            Cint[],         # rmatbeg
            Cint[],         # rmatind
            Cdouble[],      # rmatval
            503,            # nodeest  TODO: what is an objective estimate for the branch? -> i set the objective to 503 as that is slightly above the objective of the root node (and i know the solution is around there)
            seqnum_p,       # seqnum_p
        )
    if status != 0
        @warn "CPXcallbackmakebranch failed with status $(status)"
    else
        @info "I submitted a new branching strategy that branches on $(b_chs), and the status was: $(status) (successful)"
    end
end

MOI.set(m, MOI.NumberOfThreads(), 1) #It crashed when running multithread... 
MOI.set(m, CPLEX.CallbackFunction(), my_callback_function) #setting the callback function through MOI

The info logs that are outputted keep repeating until CPLEX does a restart (which CPLEX just does every once in a while I guess?) Do you know if that is expected behavior, or does that mean it is not branching the way I expect it to?

Version identifier: 20.1.0.0 | 2020-11-10 | 9bedb6d68
CPXPARAM_Threads                                 1
Generic callback                                 0xfe
Tried aggregator 3 times.
MIP Presolve eliminated 1930 rows and 202 columns.
Aggregator did 13110 substitutions.
Reduced MIP has 16013 rows, 16304 columns, and 55651 nonzeros.
Reduced MIP has 4987 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.05 sec. (50.47 ticks)
Probing fixed 673 vars, tightened 8265 bounds.
Probing time = 0.08 sec. (41.28 ticks)
Tried aggregator 1 time.
MIP Presolve eliminated 1346 rows and 1346 columns.
MIP Presolve modified 8060 coefficients.
Reduced MIP has 14667 rows, 14958 columns, and 51782 nonzeros.
Reduced MIP has 4314 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.02 sec. (15.07 ticks)
Probing fixed 12 vars, tightened 2476 bounds.
Probing time = 0.08 sec. (37.24 ticks)
Clique table members: 5464.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: none, using 1 thread.
Root relaxation solution time = 4.61 sec. (4736.57 ticks)

        Nodes                                         Cuts/
   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap

      0     0      502.5709  1126                    502.5709    16619
      0     0      502.5709   504                  Cuts: 2182    18094
      0     0      502.5709   590                  Cuts: 2567    20193
      0     0      502.5709   364                   Cuts: 771    20862
      0     0      502.5709   392                  Cuts: 1245    21631
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
      0     1      502.5709   274                    502.5709    21631
Elapsed time = 44.69 sec. (44797.98 ticks, tree = 0.01 MB, solutions = 0)
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
... # it repeats a lot in between 
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
     20     1      502.5709   274                    502.5709    21631
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
...
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
     60     1      502.5709   274                    502.5709    21631
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
...
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
    200     1      502.5709   274                    502.5709    21631
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
...
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
    427     1      502.5709   274                    502.5709    21631
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
...
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
    650     1      502.5709   274                    502.5709    21631
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
...
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
    937     1      502.5709   274                    502.5709    21631
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
   2110     1      502.5709   274                    502.5709    21631
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
[ Info: I submitted a new branching strategy that branches on CartesianIndex(10, 20, 2, 5), and the status was: 0 (successful)
# Then it performs a restart...
Performing restart 1

Repeating presolve.
Tried aggregator 2 times.
MIP Presolve eliminated 24 rows and 24 columns.
MIP Presolve modified 3548 coefficients.
Aggregator did 3 substitutions.
Reduced MIP has 14640 rows, 14931 columns, and 51282 nonzeros.
Reduced MIP has 4302 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.02 sec. (21.43 ticks)
Tried aggregator 1 time.
MIP Presolve eliminated 92 rows and 92 columns.
MIP Presolve modified 593 coefficients.
Reduced MIP has 14548 rows, 14839 columns, and 51052 nonzeros.
Reduced MIP has 4256 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.02 sec. (14.91 ticks)
Represolve time = 0.17 sec. (82.43 ticks)
   2116     0      502.5709  1770                  Cuts: 3737    50313
   2116     0      502.5709   701                  Cuts: 3737    52719
   2116     0      502.5709   771                  Cuts: 3737    54692
   2116     0      502.5709  1282                  Cuts: 3737    67158
   2116     0      502.5709   820                  Cuts: 3737    69958
   2116     0      502.5709   772                  Cuts: 3737    72132

[ Info: I submitted a new branching strategy that branches on CartesianIndex(12, 20, 2, 5), and the status was: 0 (successful)
... #after restart it finds a new variable

Also, since the variables that I am branching on are binary, I want to create two child nodes. One where it is set to 0 and the other to 1. Should I explicitly tell CPLEX to make two child nodes, by doing something like this:

seqnum_p = Ref{Cint}(0)
    status = CPXcallbackmakebranch(
        cb_data, 
        2,              # varcnt
        Cint[column, column],   # varind UPDATE: added the same binary (column) twice, one call for upper and one for lower bound. 
        Cchar['L', 'U'],     # varlu UPDATE: added 'L'
        Cdouble[0.0, 1.0],   # varbd    UPDATE: ADDED 0  to the lower bound
        0,              # rcnt
        0,              # nzcnt
        Cdouble[],      # rhs
        Cchar[],        # sense
        Cint[],         # rmatbeg
        Cint[],         # rmatind
        Cdouble[],      # rmatval
        503,            # nodeest 
        seqnum_p,       # seqnum_p
    )

or does it already implicitly know that it should create the other child node, or is there something completely different I need to do?

Last thing, which is a bit off-topic. Do you know if I can supply CPLEX with a partial solution to some binaries, together with the branching strategy? due to the way the mixed integer constraints are made, I know that a few (maybe 10-15 out of all the fractional binary variables) could be set to integer values, as they already satisfy the associated constraints (but the solver does not explicitly know this when they reach the MIPnode). I did that previously in Gurobi using a Heuristic callback, which seemed to work, but I don’t know if can both provide a solver-independent heuristic callback to CPLEX on top of the solver-dependent branching callback that I am implementing now, but maybe CPLEX has some other way to do it?

Either way, thank you SO MUCH for your help so far, it’s highly appreciated!

odow · May 19, 2023, 10:01pm

I’ve made a mistake. It should be column = CPLEX.column(cb_data, index(m[:x])) - 1. See

github.com

jump-dev/CPLEX.jl/blob/ccde44778f93d65af3c6492e44571af03bc015cc/src/MOI/MOI_wrapper.jl#L684-L691


      
          """
              column(model::Optimizer, x::MOI.VariableIndex)
          
          Return the 1-indexed column associated with `x`.
          
          The C API requires 0-indexed columns.
          """
          function column(model::Optimizer, x::MOI.VariableIndex)

But also:

Cchar[‘U’], # varlu
Cdouble[1.0], # varbd #need to make two child nodes…

If you want to force the variable to 1, don’t you need to set the lower bound to 1?

or does it already implicitly know that it should create the other child node, or is there something completely different I need to do?

Not sure. You’d have to carefully read the IBM documentation. You probably need to make two calls to CPXcallbackmakebranch?

but I don’t know if can both provide a solver-independent heuristic callback to CPLEX on top of the solver-dependent branching callback that I am implementing now, but maybe CPLEX has some other way to do it?

Modify the if statement to include

if context_id == CPX_CALLBACKCONTEXT_RELAXATION
    # ... heuristic stuff ...
elseif context_id == CPX_CALLBACKCONTEXT_BRANCHING
    # ... branching stuff ...
end

you can see how we post the heuristic solution here https://github.com/jump-dev/CPLEX.jl/blob/ccde44778f93d65af3c6492e44571af03bc015cc/src/MOI/MOI_callbacks.jl#L350-L357

valdemarsoegaard · May 23, 2023, 2:01pm

Thanks! I also found another mistake own my own part, but now i am branching on the correct variable, and i get the expected output from the callback!

and yes, you’re right, thanks!

I found a pseudo-code example, and here they do indeed make two separate calls, so I’m doing that too now and everything seems to work!

Thanks! i’ll try this too.

valdemarsoegaard · May 23, 2023, 6:40pm

I tried using this, but the default solving strategy ‘CPXCALLBACKSOLUTION_SOLVE’ does not seem to work very well in my case, where I am only providing a partial solution (it is stuck here for a very long time, probably trying to solve the partial solution to a feasible one, which is tough for my MIP).

I tried to replicate the heuristic call that you referred me to ret = CPXcallbackpostheursoln( cb.callback_data, Cint(length(variables)), Cint[_info(model, var).column - 1 for var in variables], values, NaN, CPXCALLBACKSOLUTION_SOLVE, ) and changed the solving strategy to CPXCALLBACKSOLUTION_PROPEGATE. I do not completely understand how this works from reading their description of the strategy, but at least it runs now… Do you know if makes sense that changing to this strategy works better in my case, where finding feasible solutions is tough, and I only try to fix a few fractional variables before starting my custom branching strategy starts?

odow · May 23, 2023, 7:01pm

Your analysis makes sense. In general heuristic callback are only useful if you can provide a feasible point for all integer variables.