Solving bilinear optimization problem using JuMP

Shuvomoy Das Gupta

March 9, 2022

In this blog, we will discuss how to solve a bilinear optimization problem in JuMP+Julia using Gurobi and /orKNITRO as the solver. Here we note that, Gurobi is able to solve these problems to global optimality, whereas KNITRO can provide only solve these problems only locally.

Table of contents

Problem in consideration
JuMP code
Solve from scratch
Warm-starting Gurobi with a known feasible solution
Warm-starting via JuMP
Warm-starting via Gurobi.jl
Finding a locally optimal solution using KNITRO
Some helpful online links on speeding up the solution process

Problem in consideration

We want to solve the following bilinear problem:

\begin{array}{ll} \textrm{maximize} & x_{1}x_d - x_2 x_{d-1}\\ \textrm{subject to} & \sum_{i=1}^{d}x_{i}\leq10\\ & x_{i}x_{i+1}\leq2,\quad i=1,\ldots,d-1\\ & \sum_{i=1}^{d-1}x_{i}x_{i+1}=1\\ & x\geq0, \end{array}

where $x\in\mathbf{R}^d$ is the decision variables.

`JuMP` code

The JuMP code is as follows.

using Gurobi, JuMP, BenchmarkTools

function solve_bilinear(d; start_value_provided = :off, x_start_value = zeros(d), warm_start_type = :none, silent_flag = :off)

  # declare the model
  # -----------------

  nonlinear_model = direct_model(Gurobi.Optimizer())
  # using direct_model has results in smaller memory allocation
  # we could also use
  # nonlinear_model = Model(Gurobi.Optimizer)
  # but this requires more memory allocation

  set_optimizer_attribute(nonlinear_model, "MIPFocus", 3)
  # If you are more interested in good quality feasible solutions, you can select MIPFocus=1.
  # If you believe the solver is having no trouble finding the optimal solution, and wish to focus more
  # attention on proving optimality, select MIPFocus=2.
  # If the best objective bound is moving very slowly (or not at all), you may want to try MIPFocus=3 to focus on the bound.


  set_optimizer_attribute(nonlinear_model, "NonConvex", 2)
  # "NonConvex" => 2 tells Gurobi to use its nonconvex algorithm

  # declare the variable along with the non-negativity constraitn
  @variable(nonlinear_model, x[1:d] >=0 )

  # warm starting using Gurobi.jl
  if  start_value_provided == :on && warm_start_type == :Gurobi_direct
    for i in 1:d
        MOI.set(nonlinear_model, Gurobi.VariableAttribute("Start"), x[i], x_start_value[i])
    end
  end

  # warm starting using JuMP
  if  start_value_provided == :on && warm_start_type == :JuMP
    for i in 1:d
        set_start_value(x[i], x_start_value[i])
    end
  end

  # add objective
  @objective(nonlinear_model, Max, x[1]*x[d] - x[2]*x[d-1])

  # add the constraint Σx[i] <= 10
  @constraint(nonlinear_model, sum(x[i] for i in 1:d) <= 10)

  # add the bilinear inequality constraint ∀i  x[i]x[i+1] <= 2
  for i in 1:d-1
    @constraint(nonlinear_model, x[i]*x[i+1] <= 2)
  end

  # add the bilinear equality constraint ∑_{i} x[i]x[i+1] == 1
  @constraint(nonlinear_model, sum(x[i]*x[i+1] for i in 1:d-1) == 1)

  # optimize the model
  if silent_flag == :on
     set_silent(nonlinear_model)
   end
  optimize!(nonlinear_model)

  # store the important values and return them
  status = JuMP.termination_status(nonlinear_model)
  x_sol = JuMP.value.(x)
  obj_value = JuMP.objective_value(nonlinear_model)

  return status, x_sol, obj_value

end

Solve from scratch

d = 10
status, x_sol, obj_value =  solve_bilinear(d; start_value_provided = :off, x_start_value = zeros(d), warm_start_type = :none, silent_flag = :off)

The output is:

Continuous model is non-convex -- solving as a MIP.

Presolve time: 0.02s
Presolved: 61 rows, 22 columns, 120 nonzeros
Presolved model has 20 bilinear constraint(s)
Variable types: 22 continuous, 0 integer (0 binary)
Presolve removed 39 rows and 2 columns
Presolved: 22 rows, 20 columns, 59 nonzeros


Root relaxation: objective 4.955556e+01, 19 iterations, 0.01 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0   49.55556    0   11          -   49.55556      -     -    0s
     0     0   25.21339    0    2          -   25.21339      -     -    0s
     0     0   25.10018    0   11          -   25.10018      -     -    0s
     0     0   25.09269    0   11          -   25.09269      -     -    0s
     0     0   25.03708    0   10          -   25.03708      -     -    0s
     0     0   25.00988    0   10          -   25.00988      -     -    0s
     0     0   24.99742    0   10          -   24.99742      -     -    0s
     0     0   24.99291    0   10          -   24.99291      -     -    0s
H    0     0                    -100.0000000   24.99291   125%     -    0s
H    0     0                      23.9993796   24.99291  4.14%     -    0s
     0     2   24.99291    0   10   23.99938   24.99291  4.14%     -    0s

Cutting planes:
  RLT: 22

Explored 81 nodes (380 simplex iterations) in 0.18 seconds
Thread count was 8 (of 8 available processors)

Solution count 2: 23.9994 -100

Optimal solution found (tolerance 1.00e-04)
Best objective 2.399937963333e+01, best bound 2.400151700222e+01, gap 0.0089%

User-callback calls 298, time in user-callback 0.01 sec

Let us benchmark the results.

b1 = @benchmark solve_bilinear(d; start_value_provided = :off, x_start_value = zeros(d), warm_start_type = :none,  silent_flag = :on)
println("benchmark for original code without start value")
println("***************************")
io = IOBuffer()
show(io, "text/plain", b1)
s = String(take!(io))
println(s)

Benchmarked output is as follows.

benchmark for original code without start value
***************************
BenchmarkTools.Trial:
  memory estimate:  48.24 KiB
  allocs estimate:  884
  --------------
  minimum time:     45.491 ms (0.00% GC)
  median time:      55.348 ms (0.00% GC)
  mean time:        56.950 ms (0.00% GC)
  maximum time:     70.355 ms (0.00% GC)
  --------------
  samples:          88
  evals/sample:     1

Warm-starting `Gurobi` with a known feasible solution

What if we have a good guess of a solution for this problem? In that case, we can feed that vaue using set_start_value(variable_name, variable_start_value) function or MOI.set(nonlinear_model, Gurobi.VariableAttribute("Start"), variable_name, variable_start_value). To test this, we just feed the optimal solution that we just found to test if the warm-starting is working.

Warm-starting via `JuMP`

status_2, x_sol_2, obj_value_2 = solve_bilinear(d; start_value_provided = :on, x_start_value = x_sol, warm_start_type = :JuMP, silent_flag = :off)

The output is:

Continuous model is non-convex -- solving as a MIP.


Loaded user MIP start with objective 23.9994

Presolve time: 0.00s
Presolved: 61 rows, 22 columns, 120 nonzeros
Presolved model has 20 bilinear constraint(s)
Variable types: 22 continuous, 0 integer (0 binary)
Presolve removed 39 rows and 2 columns
Presolved: 22 rows, 20 columns, 59 nonzeros


Root relaxation: objective 4.955556e+01, 19 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0   49.55556    0   11   23.99938   49.55556   106%     -    0s
     0     0   25.21341    0    2   23.99938   25.21341  5.06%     -    0s
     0     0   25.09295    0    8   23.99938   25.09295  4.56%     -    0s
     0     0   25.03730    0   10   23.99938   25.03730  4.32%     -    0s
     0     0   25.01008    0   11   23.99938   25.01008  4.21%     -    0s
     0     0   24.99757    0    6   23.99938   24.99757  4.16%     -    0s
     0     0   24.99757    0   10   23.99938   24.99757  4.16%     -    0s
     0     0   24.99299    0    7   23.99938   24.99299  4.14%     -    0s
     0     2   24.99299    0    7   23.99938   24.99299  4.14%     -    0s

Cutting planes:
  RLT: 22

Explored 81 nodes (440 simplex iterations) in 0.08 seconds
Thread count was 8 (of 8 available processors)

Solution count 1: 23.9994

Optimal solution found (tolerance 1.00e-04)
Best objective 2.399937963333e+01, best bound 2.400120366930e+01, gap 0.0076%

User-callback calls 274, time in user-callback 0.00 sec

b2 = @benchmark solve_bilinear(d; start_value_provided = :on, x_start_value = x_sol, warm_start_type = :JuMP, silent_flag = :on)
println("benchmark for original code with start value via JuMP")
println("***************************")
io = IOBuffer()
show(io, "text/plain", b2)
s = String(take!(io))
println(s)

The benchmarking result is as follows.

benchmark for original code with start value via JuMP
***************************
BenchmarkTools.Trial:
  memory estimate:  48.43 KiB
  allocs estimate:  904
  --------------
  minimum time:     43.148 ms (0.00% GC)
  median time:      47.913 ms (0.00% GC)
  mean time:        49.426 ms (0.00% GC)
  maximum time:     84.259 ms (0.00% GC)
  --------------
  samples:          102
  evals/sample:     1

Warm-starting via `Gurobi.jl`

status_3, x_sol_3, obj_value_3 = solve_bilinear(d; start_value_provided = :on, x_start_value = x_sol, warm_start_type =  :Gurobi_direct, silent_flag = :off)

The output is:

Continuous model is non-convex -- solving as a MIP.


Loaded user MIP start with objective 23.9994

Presolve time: 0.00s
Presolved: 61 rows, 22 columns, 120 nonzeros
Presolved model has 20 bilinear constraint(s)
Variable types: 22 continuous, 0 integer (0 binary)
Presolve removed 39 rows and 2 columns
Presolved: 22 rows, 20 columns, 59 nonzeros


Root relaxation: objective 4.955556e+01, 19 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0   49.55556    0   11   23.99938   49.55556   106%     -    0s
     0     0   25.21341    0    2   23.99938   25.21341  5.06%     -    0s
     0     0   25.09295    0    8   23.99938   25.09295  4.56%     -    0s
     0     0   25.03730    0   10   23.99938   25.03730  4.32%     -    0s
     0     0   25.01008    0   11   23.99938   25.01008  4.21%     -    0s
     0     0   24.99757    0    6   23.99938   24.99757  4.16%     -    0s
     0     0   24.99757    0   10   23.99938   24.99757  4.16%     -    0s
     0     0   24.99299    0    7   23.99938   24.99299  4.14%     -    0s
     0     2   24.99299    0    7   23.99938   24.99299  4.14%     -    0s

Cutting planes:
  RLT: 22

Explored 81 nodes (440 simplex iterations) in 0.09 seconds
Thread count was 8 (of 8 available processors)

Solution count 1: 23.9994

Optimal solution found (tolerance 1.00e-04)
Best objective 2.399937963333e+01, best bound 2.400120366930e+01, gap 0.0076%

User-callback calls 267, time in user-callback 0.00 sec

Let's take a look at the benchmarking output for this case.

b3 = @benchmark solve_bilinear(d; start_value_provided = :on, x_start_value = x_sol, warm_start_type = :Gurobi_direct, silent_flag = :on)
println("benchmark for original code with start value via Gurobi.jl")
println("***************************")
io = IOBuffer()
show(io, "text/plain", b3)
s = String(take!(io))
println(s)

The output is as follows this time.

benchmark for original code with start value via Gurobi.jl
***************************
BenchmarkTools.Trial:
  memory estimate:  48.46 KiB
  allocs estimate:  910
  --------------
  minimum time:     31.733 ms (0.00% GC)
  median time:      36.254 ms (0.00% GC)
  mean time:        37.185 ms (0.00% GC)
  maximum time:     73.067 ms (0.00% GC)
  --------------
  samples:          135
  evals/sample:     1

We see that in both warm-starting procedures, Gurobi accepts the provided solution, which makes the branch and bound method work somewhat better (number of explored nodes goes down from 397 to 225 after warm-starting). Of course, because the problem size is small, it does not impact the runtime significantly.

Finding a locally optimal solution using `KNITRO`

using JuMP, KNITRO
function solve_bilinear_locally_KNITRO(d; start_value_provided = :off, x_start_value = zeros(d), warm_start_type = :none, silent_flag = :off)

  # declare the model
  # -----------------

  nonlinear_model = Model(optimizer_with_attributes(KNITRO.Optimizer, "convex" => 0))

  # declare the variable along with the non-negativity constraitn
  @variable(nonlinear_model, x[1:d] >=0 )


  # warm starting using JuMP
  if  start_value_provided == :on
    for i in 1:d
        set_start_value(x[i], x_start_value[i])
    end
  end

  # add objective
  @objective(nonlinear_model, Max, x[1]*x[d] - x[2]*x[d-1])

  # add the constraint Σx[i] <= 10
  @constraint(nonlinear_model, sum(x[i] for i in 1:d) <= 10)

  # add the bilinear inequality constraint ∀i  x[i]x[i+1] <= 2
  for i in 1:d-1
    @constraint(nonlinear_model, x[i]*x[i+1] <= 2)
  end

  # add the bilinear equality constraint ∑_{i} x[i]x[i+1] == 1
  @constraint(nonlinear_model, sum(x[i]*x[i+1] for i in 1:d-1) == 1)

  # optimize the model
  if silent_flag == :on
     set_silent(nonlinear_model)
   end
  optimize!(nonlinear_model)

  # store the important values and return them
  status = JuMP.termination_status(nonlinear_model)
  x_sol = JuMP.value.(x)
  obj_value = JuMP.objective_value(nonlinear_model)

  return status, x_sol, obj_value

end

Solve the problem to local optimality using KNITRO.

d = 10
status_KNITRO, x_sol_KNITRO, obj_value_KNITRO = solve_bilinear_locally_KNITRO(d; start_value_provided = :off, x_start_value = zeros(d), warm_start_type = :none, silent_flag = :off)

The output is:

=======================================
             Trial License
       (NOT FOR COMMERCIAL USE)
         Artelys Knitro 12.3.0
=======================================

Knitro presolve eliminated 0 variables and 0 constraints.

convex:                  0
datacheck:               0
hessian_no_f:            1
The problem is identified as a QCQP.

Problem Characteristics                                 (   Presolved)
-----------------------
Objective goal:  Maximize
Objective type:  quadratic
Number of variables:                                 10 (          10)
    bounded below only:                              10 (          10)
    bounded above only:                               0 (           0)
    bounded below and above:                          0 (           0)
    fixed:                                            0 (           0)
    free:                                             0 (           0)
Number of constraints:                               11 (          11)
    linear equalities:                                0 (           0)
    quadratic equalities:                             1 (           1)
    gen. nonlinear equalities:                        0 (           0)
    linear one-sided inequalities:                    1 (           1)
    quadratic one-sided inequalities:                 9 (           9)
    gen. nonlinear one-sided inequalities:            0 (           0)
    linear two-sided inequalities:                    0 (           0)
    quadratic two-sided inequalities:                 0 (           0)
    gen. nonlinear two-sided inequalities:            0 (           0)
Number of nonzeros in Jacobian:                      38 (          38)
Number of nonzeros in Hessian:                       11 (          11)

Knitro using the Interior-Point/Barrier Direct algorithm.

  Iter      Objective      FeasError   OptError    ||Step||    CGits
--------  --------------  ----------  ----------  ----------  -------
       0    0.000000e+00   1.000e+00
      10    9.413565e-01   8.790e-02   8.971e-01   1.222e-01        0
      20    2.397843e+01   1.159e-05   9.642e-03   1.307e-02        0
      25    2.400000e+01   6.485e-07   9.957e-08   1.186e-03        0

EXIT: Locally optimal solution found.

Final Statistics
----------------
Final objective value               =   2.39999992955648e+01
Final feasibility error (abs / rel) =   6.48e-07 / 6.48e-07
Final optimality error  (abs / rel) =   9.96e-08 / 1.99e-08
# of iterations                     =         25
# of CG iterations                  =          9
# of function evaluations           =          0
# of gradient evaluations           =          0
# of Hessian evaluations            =          0
Total program time (secs)           =       0.029 (     0.000 CPU time)
Time spent in evaluations (secs)    =       0.000

Note that KNITRO is providing the same objective value as Gurobi, but is not able certify the global optimality.