# 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.

### 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

@objective(nonlinear_model, Max, x*x[d] - x*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.

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.

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

@objective(nonlinear_model, Max, x*x[d] - x*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:

=======================================
(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
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.