Juniper behaviour: No incumbent when MIP feasible?

Specific Domains Optimization (Mathematical)
1

I’m getting some unintuitive behaviour with Juniper (HiGHS & Ipopt). HiGHS is finding a solution that satisfies the constraints, but beyond that an incumbent solution won’t be reported and Juniper will time-out at an infeasible point. If I disable the objective function, I do get a valid solution. I feel like that warning from the feasibility pump code is important, but am unsure what to do about it.

Julia output & reproducible code follow.

#Variables: 1752
#IntBinVar: 1752
Obj Sense: Min

Start values are not feasible.
Status of relaxation: LOCALLY_SOLVED
Time for relaxation: 4.153654098510742
Relaxation Obj: 1.756802304325599e-15

       MIPobj              NLPobj       Time 
=============================================

[..............]

Solving report
  Status            Optimal
  Primal bound      102.222208073
  Dual bound        102.211991515
  Gap               0.00999% (tolerance: 0.01%)
  Solution status   feasible
                    102.222208073 (objective)
[..............]

┌ Warning: NLP couldn't be solved to optimality
└ @ Juniper ~/.julia/packages/Juniper/HEO6p/src/fpump.jl:405
      102.2222               -          53.8 

FP: 53.75900387763977 s
FP: 1 round
FP: No integral solution found

 ONodes   CLevel          Incumbent                   BestBound            Gap    Time   Restarts  GainGap  
============================================================================================================
[ Info: Breaking out of strong branching as the time limit of 100.0 seconds got reached.
    2       2                 -                          0.0                -    102.6      0         -     
    3       3                 -                          0.0                -    107.4      -      5867.1%

[..............]

#branches: 27
BnB time: 244.76
% solve child time: 100.0
Solve node time get idx: 102.6
Solve node time branch: 142.15
Branch time: 142.15
Get idx time: 102.6
Upd gains time: 0.0
Obj: NaN
termination_status(model) = MathOptInterface.TIME_LIMIT
primal_status(model) = MathOptInterface.INFEASIBLE_POINT

Code sample (apologies for length).

import HiGHS
import Ipopt
import Juniper

using JuMP

# Convenience
function _to_bool(m)::Matrix{Bool}
    return [x == true for x in m]
end

function _test()

    # Juniper-based approach
    optimizer = Juniper.Optimizer
    nl_solver= optimizer_with_attributes(Ipopt.Optimizer,
                                         "print_level" => 0)
    mip_solver = optimizer_with_attributes(HiGHS.Optimizer,
                                           "output_flag" => true)
    model = Model(optimizer_with_attributes(
        optimizer,
        "nl_solver"=>nl_solver,
        "mip_solver"=>mip_solver,
        "log_levels"=>[:Info, :Table, :Timing, :AllOptions]))
                                            
    set_time_limit_sec(model, 300)

    num_teams = 6
    stints_per_session = 6
    hours_per_stint = 2
    
    # These work
    # stints_per_session = 2
    # hours_per_stint = 1

    # There are 4 overlapping "sessions"
    # over a 30 hour period, each having 6 stints.
    sessions = _to_bool(
        [1 1 1 1 1 1 0 0 0 0 0 0 0 0 0
         0 0 0 0 1 1 1 1 1 1 0 0 0 0 0
         0 0 0 0 0 0 0 1 1 1 1 1 1 0 0
         0 0 0 0 0 0 0 0 0 1 1 1 1 1 1])
    
    (num_sessions, event_length) = size(sessions)

    # Driver availability for each stint
    avail_matrix = _to_bool(
        [1  1  1  1  0  0  0  0  0  1  1  1  1  1  1
         0  1  1  1  1  1  1  1  1  1  0  0  0  0  0
         1  1  1  1  0  0  0  1  1  1  1  1  1  1  1
         0  1  1  1  0  0  0  0  1  1  1  1  1  1  1
         1  1  1  1  0  0  0  0  0  1  1  1  1  1  1
         0  1  1  1  1  0  0  0  1  1  1  1  1  1  1
         0  0  1  1  1  1  1  1  1  1  1  1  1  1  1
         0  0  1  1  1  0  0  1  0  0  0  0  0  0  1
         1  1  1  1  1  0  0  0  1  1  1  1  1  1  1
         0  0  0  0  0  0  0  1  1  1  1  1  0  0  1
         0  0  0  0  0  0  0  1  1  1  1  1  1  1  1
         0  1  1  1  0  0  0  1  1  1  1  1  1  1  1
         0  0  0  1  1  1  1  0  0  0  0  0  0  0  0
         0  0  0  1  1  1  1  1  1  1  1  1  1  0  0
         0  0  0  0  1  1  1  1  1  1  0  0  0  0  0
         1  1  0  0  1  1  0  0  0  1  1  1  1  0  0
         1  1  1  1  1  0  0  0  1  1  1  0  0  0  0
         1  1  1  1  1  1  1  1  1  1  1  1  1  1  1])

    num_drivers = size(avail_matrix)[1]
    println("  num_drivers = " * string(num_drivers))

    # Fix everyone's ideal hours for now.
    ideal_hours = ones(num_drivers) * 4

    # A schedule is a boolean matrix of:
    # drivers x stints x teams
    @variable(model,
              schedule[1:num_drivers, 1:event_length, 1:num_teams], Bin)

    # A driver cannot be scheduled if they are not available
    # and can only be scheduled for one team in a given stint.
    @constraint(model,
                sum(schedule, dims=3) .<= avail_matrix)

    # team_sessions: num_sessions x num_teams of booleans
    #                True for the session the team is running.
    @variable(model,
              team_sessions[1:num_sessions, 1:num_teams], Bin)
    # A team runs in one and only one session.
    @constraint(model,
                team_sessions' * ones(Bool, num_sessions)
                .== ones(Bool, num_teams))

    # A team must have exactly 1 driver for each stint in their session
    @constraint(model,
            [i in 1:num_teams],
            schedule[:,:,i]' * ones(Bool, num_drivers)
            .== sessions' * team_sessions[:,i])

    # A utility variable for team composition.
    @variable(model,
              team_drivers[1:num_drivers, 1:num_teams], Bin)

    # A driver for a team must drive at least 1 stint.
    @constraint(model,
                [i in 1:num_teams],
                sum(schedule[:,:,i], dims=(2)) .>= team_drivers[:,i])
    # And (double-duty) (big-M)
    # 1. A driver on the team should drive less than stints_per_session
    # 2. (more important) A driver (x) _not_ on the team does not
    #    have a scheduled stint.
    #       (team_drivers[x,i] == 0 -> schedule <= 0
    @constraint(model,
                [i in 1:num_teams],
                sum(schedule[:,:,i], dims=(2)) .<=
                    team_drivers[:,i] * stints_per_session)

    # Leverage the team_drivers to constrain team size
    # Limit a team to max_team_size distinct drivers
    min_team_size = 2
    max_team_size = 3
    @constraint(model,
                [i in 1:num_teams],
                sum(team_drivers[:,i]) <= max_team_size)
    @constraint(model,
                [i in 1:num_teams],
                sum(team_drivers[:,i]) >= min_team_size)

    # A driver is on exactly one team.
    #@constraint(model, team_drivers * ones(num_teams) .== ones(num_drivers))
    # A driver is on one or two teams.
    @constraint(model, team_drivers * ones(num_teams) .>= ones(num_drivers))
    @constraint(model, team_drivers * ones(num_teams) .<= ones(num_drivers)*2)

    # Is this an optimization? Pinning the first driver to a team?
    # @constraint(model, team_drivers[1,1] == true)
    fix(team_drivers[1, 1], true; force = true)

    # Minimize sum of squares of scheduled hours - ideal hours
    # for each driver.
    hours_err_pt_1 = @expression(model,
       sum(schedule, dims=(2,3)) * hours_per_stint - ideal_hours)

    # TODO register .^ ?
    if true
        @NLobjective(model, Min,
                     hours_err_pt_1[1] ^ 2
                     + hours_err_pt_1[2] ^ 2
                     + hours_err_pt_1[3] ^ 2
                     + hours_err_pt_1[4] ^ 2
                     + hours_err_pt_1[5] ^ 2
                     + hours_err_pt_1[6] ^ 2
                     + hours_err_pt_1[7] ^ 2
                     + hours_err_pt_1[8] ^ 2
                     + hours_err_pt_1[9] ^ 2
                     + hours_err_pt_1[10] ^ 2
                     + hours_err_pt_1[11] ^ 2
                     + hours_err_pt_1[12] ^ 2
                     + hours_err_pt_1[13] ^ 2
                     + hours_err_pt_1[14] ^ 2
                     + hours_err_pt_1[15] ^ 2
                     + hours_err_pt_1[16] ^ 2
                     + hours_err_pt_1[17] ^ 2
                     + hours_err_pt_1[18] ^ 2)
    end

    optimize!(model)

    @show termination_status(model)
    
    @show primal_status(model)
    
    solution_summary(model)

    return model
end

@time model = _test()
2

The feasibility pump attempts to find a valid assignment of the discrete variables using the MIP solver and an approximation of the nonlinear constraints, then the NLP solver is used to try to find a solution to the nonlinear system. In this case it seems the first phase works but then the second phase fails.

Having a quick look at your model I might try adding a variable for hours_err_pt_1 rather than an expression. Also your model looks like it falls into the category MIQP, if so, this should be better handled by solvers like Guorbi, CPlex, Express or maybe even Pavito.

Juniper is most useful for MINLP models with many non-convex constraints that more complex than quadratic.

3

Thank you Carleton. Changing hours_err_pt_1 to be a variable worked. I think I understand intuitively why that helped, but not why it was necessary. If you want to recommend a reference that might help someone like me get up-to-speed quickly, feel free to do so.

Maybe I’ll look at Pavito again (Guorbi, CPlex, and Express are all commercial — which, is fine, but I’m checking out the free stuff first :wink: )

Thanks again.

4

I unfortunately don’t know a good reference for the art of NLP modeling. I suppose it is still a field that is maturing. The rules of thumb that I use is, try to avoid having equality constraints in NLP (i.e. project out any intermediate variables with expressions, as you have done here); this is mostly to improve runtime performance. If you run into convergence issues: (1) upgrade to the linear system solver, HSL tends to be a big improvement over the open-access alternatives, in my experience; (2) then try carefully lifting the model (by adding back some variables with equality expressions) to make the objective and constraints linear or quadratic. This tends to improve reliability at the cost of slower convergence.

One nice example of how the lifting can be problematic to performance can be seen in this paper, [2005.14087] The Impacts of Convex Piecewise Linear Cost Formulations on AC Optimal Power Flow, where we show that different linear constraint forms can have a big impact on NLP convergence.

2 Likes