Optimizing a quadratic binary function on CPLEX/JuMP

I am using CPLEX 12.8 in JuMP to minimize a quadratic binary nonconvex function, according to quadratic function by CPLEX.
In particular, my function is the following:

\sum_{i=1}^{m-1} \sum_{f=1}^{F} \sum_{\underset{\bar{f} \neq f}{\bar{f}=1}}^F \sum_{h \in H} \left( D_{f \bar{f}} \cdot \gamma_{i\bar{f}h} \cdot \gamma_{i+1,f,h} \right)

over a set X \neq \varnothing, where \gamma_{ifh} \in \{0,1\} and D_{f\bar{f}} > 0 (distance symetric matrix).

However, when I apply CPLEX in this function (even not knowing if it is convex), the error appears:

Matrix Q must be either symmetric or triangular
error(::String) at error.jl:33
add_qpterms!(::CPLEX.Model, ::Array{Int32,1}, ::Array{Int32,1}, ::Array{Float64,1}) at cpx_quad.jl:17
set(::CPLEX.Optimizer, ::MathOptInterface.ObjectiveFunction{MathOptInterface.ScalarQuadraticFunction{Float64}}, ::MathOptInterface.ScalarQuadraticFunction{Float64}) at MOI_wrapper.jl:727
set(::MathOptInterface.Bridges.LazyBridgeOptimizer{CPLEX.Optimizer}, ::MathOptInterface.ObjectiveFunction{MathOptInterface.ScalarQuadraticFunction{Float64}}, ::MathOptInterface.ScalarQuadraticFunction{Float64}) at bridge_optimizer.jl:718
_pass_attributes(::MathOptInterface.Bridges.LazyBridgeOptimizer{CPLEX.Optimizer}, ::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}, ::Bool, ::MathOptInterface.Utilities.IndexMap, ::Array{MathOptInterface.AbstractModelAttribute,1}, ::Tuple{}, ::Tuple{}, ::Tuple{}, ::typeof(MathOptInterface.set)) at copy.jl:148
pass_attributes at copy.jl:112 [inlined]
pass_attributes at copy.jl:111 [inlined]
default_copy_to(::MathOptInterface.Bridges.LazyBridgeOptimizer{CPLEX.Optimizer}, ::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}, ::Bool) at copy.jl:337
#automatic_copy_to#109 at copy.jl:15 [inlined]
automatic_copy_to at copy.jl:14 [inlined]
#copy_to#3 at bridge_optimizer.jl:268 [inlined]
(::MathOptInterface.var"#copy_to##kw")(::NamedTuple{(:copy_names,),Tuple{Bool}}, ::typeof(MathOptInterface.copy_to), ::MathOptInterface.Bridges.LazyBridgeOptimizer{CPLEX.Optimizer}, ::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}) at bridge_optimizer.jl:268
attach_optimizer(::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}) at cachingoptimizer.jl:149
optimize!(::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}) at cachingoptimizer.jl:185
optimize!(::Model, ::Nothing; bridge_constraints::Bool, ignore_optimize_hook::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at optimizer_interface.jl:131
optimize! at optimizer_interface.jl:107 [inlined]
optimize!(::Model) at optimizer_interface.jl:107
top-level scope at Modelo_avancado_graus_dias.jl.jl:139
include_string(::Function, ::Module, ::String, ::String) at loading.jl:1088

I have linerized this function by adding the constraints:

\gamma_{i\bar{f}h} + \gamma_{i+1,f,h} - u_{if\bar{fh}} \leq 1, \quad \quad i=1,...,m,\, \quad f=1,...,F, \quad \bar{f} \in \{1,...,F\} \setminus f, \quad h \in H
u_{if\bar{fh}} \leq \gamma_{i\bar{f}h}, \quad \quad i=1,...,m,\, \quad f=1,...,F, \quad \bar{f} \in \{1,...,F\} \setminus f, \quad h \in H
u_{if\bar{fh}} \leq \gamma_{i+1,f,h}, \quad \quad i=1,...,m,\, \quad f=1,...,F, \quad \bar{f} \in \{1,...,F\} \setminus f, \quad h \in H

where 0 \leq u_{if\bar{fh}} = \gamma_{i\bar{f}h} \cdot \gamma_{i+1,f,h}. But it was terrible to optimize the linear version:

\sum_{i=1}^{m-1} \sum_{f=1}^{F} \sum_{\underset{\bar{f} \neq f}{\bar{f}=1}}^F \sum_{h \in H} D_{f \bar{f}} \cdot u_{if\bar{f}h}.

Questions:

  1. How can I fix this?
  2. How can I acess the matrix Q?
  3. There is another way to minimize this function by using CPLEX, even I am not knowing if is convex?

Please provide a link when cross-posting: mixed integer programming - Minimizing a quadratic binary nonconvex function by CPLEX - Operations Research Stack Exchange

I’ve answered your question there.

Thank you. It’s working now. The problem was solved.
However, I noted a weird fact: when I am optimizing the quadratic objective and I use initial feasible solution, like this:

set_start_value(variable, value)

the algorithm ignores my solution That’s true? How can I insert my initial feasible solution?

Look my CPLEX report (I provided a complete integer feasible solution):

CPXPARAM_OptimalityTarget                        3
CPXPARAM_Emphasis_MIP                            1
CPXPARAM_TimeLimit                               120
CPXPARAM_MIP_Tolerances_MIPGap                   0.01
9 of 10 MIP starts provided solutions.
MIP start 'm9' defined initial solution with objective 9402.0000.
Tried aggregator 1 time.
MIQP Presolve eliminated 65292 rows and 3399 columns.
MIQP Presolve modified 5557 coefficients.
Reduced MIQP has 4800 rows, 17018 columns, and 71282 nonzeros.
Reduced MIQP has 11039 binaries, 49 generals, 0 SOSs, and 0 indicators.
Reduced MIQP objective Q matrix has 5760 nonzeros.
Presolve time = 0.17 sec. (184.33 ticks)
Probing fixed 1988 vars, tightened 1853 bounds.
Probing changed sense of 1 constraints.
Probing time = 0.27 sec. (154.62 ticks)
Cover probing fixed 46 vars, tightened 321 bounds.
Clique table members: 4870.
MIP emphasis: integer feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 12 threads.
Root relaxation solution time = 0.28 sec. (329.15 ticks)

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

      0     0      840.4286   828                    840.4286       39
      0     0      876.3884   667                   Cuts: 242    13737
      0     0      890.9483   675                  Cuts: 3022    29421
      0     0      902.6731   474                  Cuts: 1627    62699
      0     2      905.5138   475                    902.6731    80836
Elapsed time = 79.89 sec. (105305.63 ticks, tree = 0.02 MB, solutions = 9)
      1     3      911.2117   894                    905.5138   103486

GUB cover cuts applied:  13
Clique cuts applied:  192
Cover cuts applied:  121
Implied bound cuts applied:  1570
Flow cuts applied:  26
Mixed integer rounding cuts applied:  44
Zero-half cuts applied:  53
Gomory fractional cuts applied:  2

Root node processing (before b&c):
  Real time             =   65.50 sec. (86743.34 ticks)
Parallel b&c, 12 threads:
  Real time             =   54.63 sec. (57737.96 ticks)
  Sync time (average)   =   44.66 sec.
  Wait time (average)   =    0.00 sec.
                          ------------
Total (root+branch&cut) =  120.13 sec. (144481.30 ticks)

9 of 10 MIP starts provided solutions.
MIP start ‘m9’ defined initial solution with objective 9402.0000.

This looks like it received an initial solution. What makes you think it didn’t?

set_start_value is the correct way to provide an initial solution.

To say more, you will need to provide a minimal working example.

1 Like

I don’t think I considered my initial solution because the “Best Integer” column is blank.
At the end of the round, when I enter
has_values(model)
the result is false, i.e. no integer feasible solution has been determined.

When I change this objective (quadratic) to a linear one, and I provide the same initial solution, the CPLEX report starts from this solution. Note the CPLEX report below:

Version identifier: 20.1.0.0 | 2020-11-10 | 9bedb6d68
CPXPARAM_OptimalityTarget                        3   
CPXPARAM_Preprocessing_QToLin                    1   
CPXPARAM_TimeLimit                               600 
CPXPARAM_MIP_Tolerances_MIPGap                   0.01
1 of 1 MIP starts provided solutions.
MIP start 'm1' defined initial solution with objective 24.0000.
Tried aggregator 1 time.
MIQP Presolve eliminated 14766 rows and 6407 columns.
MIQP Presolve modified 10208 coefficients.
Reduced MIQP has 2757 rows, 11130 columns, and 56016 nonzeros.
Reduced MIQP has 11031 binaries, 49 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.09 sec. (129.10 ticks)
Probing fixed 938 vars, tightened 88 bounds.
Probing time = 0.05 sec. (66.39 ticks)
Cover probing fixed 1 vars, tightened 100 bounds.
Tried aggregator 3 times.
MIP Presolve eliminated 315 rows and 3029 columns.
MIP Presolve modified 1411 coefficients.
Aggregator did 9 substitutions.
Reduced MIP has 2432 rows, 8092 columns, and 40846 nonzeros.
Reduced MIP has 8001 binaries, 56 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.05 sec. (59.63 ticks)
Probing time = 0.00 sec. (4.42 ticks)
Tried aggregator 1 time.
Detecting symmetries...
Reduced MIP has 2432 rows, 8092 columns, and 40846 nonzeros.
Reduced MIP has 8001 binaries, 56 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.02 sec. (29.08 ticks)
Probing time = 0.02 sec. (6.80 ticks)
Clique table members: 2772.
MIP emphasis: balance optimality and feasibility.    
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 12 threads.
Root relaxation solution time = 0.11 sec. (172.81 ticks)

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

*     0+    0                           24.0000        1.0000            95.83%
      0     0        7.0000   622       24.0000        7.0000     1859   70.83%
*     0+    0                           18.0000        7.0000            61.11%
      0     0        7.0000   620       18.0000     Cuts: 138     2883   61.11%
      0     0        7.0000   646       18.0000     Cuts: 229     4389   61.11%
      0     0        7.0000   519       18.0000      Cuts: 35     4840   61.11%
*     0+    0                           15.0000        7.0000            53.33%
*     0+    0                           14.0000        7.0000            50.00%
*     0+    0                           12.0000        7.0000            41.67%
      0     0  -1.00000e+75     0       12.0000        7.0000     4840   41.67%
      0     0        7.0000   655       12.0000     Cuts: 148     5814   41.67%
      0     2        7.0000   463       12.0000        7.1891     5814   40.09%
Elapsed time = 3.01 sec. (3202.91 ticks, tree = 0.02 MB, solutions = 5)
      6     8        7.9908   473       12.0000        7.1891     9390   40.09%
    139    90        9.5716   421       12.0000        7.1891    29500   40.09%
    348   202       10.8446   274       12.0000        7.1891    53508   40.09%
    567   327        cutoff             12.0000        7.4073    79388   38.27%
    661   437        cutoff             12.0000        7.4073   104480   38.27%
    804   483       10.0126   275       12.0000        7.4408   117518   37.99%
    939   591       10.5889   319       12.0000        7.6040   144016   36.63%
   1081   606        9.4408   487       12.0000        7.6040   149197   36.63%
   1155   694        8.4930   409       12.0000        7.6040   180246   36.63%
   1592  1008        9.6312   320       12.0000        8.0617   290219   32.82%
Elapsed time = 8.50 sec. (6334.62 ticks, tree = 13.39 MB, solutions = 5)
   2169  1296       10.1874   371       12.0000        8.3531   379319   30.39%
   2956  1772       10.9632   297       12.0000        8.4911   471490   29.24%
   3583  2106        9.6416   321       12.0000        8.6525   564238   27.90%

Performing restart 1

Repeating presolve.
Tried aggregator 1 time.
MIP Presolve eliminated 0 rows and 27 columns.
MIP Presolve modified 764 coefficients.
Reduced MIP has 2432 rows, 8065 columns, and 40840 nonzeros.
Reduced MIP has 7974 binaries, 56 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.03 sec. (35.82 ticks)
Tried aggregator 1 time.
Reduced MIP has 2432 rows, 8065 columns, and 40840 nonzeros.
Reduced MIP has 7974 binaries, 56 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.01 sec. (29.74 ticks)
Represolve time = 0.19 sec. (145.37 ticks)
   3919     0        7.0000   504       12.0000     Cuts: 117   628916   27.57%
   3919     0        7.0000   635       12.0000     Cuts: 213   630133   27.57%
   3919     0        7.0000   681       12.0000      Cuts: 40   631631   27.57%
   3919     0        7.0180   745       12.0000     Cuts: 112   634617   27.57%
   3919     0        7.0556   615       12.0000     Cuts: 122   634989   27.57%
   3919     0        7.0569   671       12.0000     Cuts: 210   635539   27.57%
   3919     2        7.0569   603       12.0000        8.6912   635539   27.57%
   3923     5        7.3976   467       12.0000        8.6912   636755   27.57%
   3927     4        7.4936   593       12.0000        8.6912   637974   27.57%
   3931     8        7.7211   483       12.0000        8.6912   643032   27.57%
   3939    10        8.1176   526       12.0000        8.6912   646686   27.57%
   3951    20        9.2834   336       12.0000        8.6912   658075   27.57%
   3969    26        9.1465   404       12.0000        8.6912   670469   27.57%
Elapsed time = 30.64 sec. (22273.34 ticks, tree = 0.04 MB, solutions = 5)
   4059    49       10.9057   331       12.0000        8.6912   704932   27.57%
   4170   105       10.7076   216       12.0000        8.6912   728846   27.57%
   4247   109        8.7460   426       12.0000        8.6912   732961   27.57%
   4383   173       10.2643   369       12.0000        8.6912   771117   27.57%
   4430   198       10.9686   396       12.0000        8.6912   781149   27.57%
   4474   257        8.8931   454       12.0000        8.6912   816886   27.57%
   4576   296       10.8352   255       12.0000        8.6912   845946   27.57%
   4744   310       10.5107   322       12.0000        8.6912   858149   27.57%
   4942   395        9.4332   337       12.0000        8.6912   887999   27.57%
   5138   517        9.6968   281       12.0000        8.6912   950382   27.57%
Elapsed time = 51.98 sec. (31970.56 ticks, tree = 1.17 MB, solutions = 5)
   5431   659        9.9065   328       12.0000        8.6912   996467   27.57%
   5864   913       10.8874   253       12.0000        8.6912  1113875   27.57%
   6128   922        cutoff             12.0000        8.6912  1120209   27.57%
   6379  1168       10.2230   361       12.0000        8.6950  1237950   27.54%
   6516  1265       10.6848   175       12.0000        8.7476  1297612   27.10%
   6816  1459        cutoff             12.0000        8.9255  1374906   25.62%
   7044  1583        9.7794   295       12.0000        8.9255  1457025   25.62%
   7408  1770       10.7136   281       12.0000        9.0320  1493785   24.73%
   7838  2009       10.5825   317       12.0000        9.0535  1566250   24.55%
   8256  2132        9.9651   359       12.0000        9.1431  1638258   23.81%
Elapsed time = 73.23 sec. (41529.41 ticks, tree = 18.42 MB, solutions = 5)
   8883  2404        cutoff             12.0000        9.1961  1744040   23.37%
   9524  2661       10.2027   391       12.0000        9.2902  1817003   22.58%
  10162  3240       10.2478   389       12.0000        9.3311  1948103   22.24%
  10565  3450       10.9238   270       12.0000        9.3590  1989334   22.01%
  11006  3697       10.4513   202       12.0000        9.4249  2078468   21.46%
  11520  3809       10.1083   408       12.0000        9.4473  2125970   21.27%
  11949  4149        cutoff             12.0000        9.4786  2229034   21.01%
  12115  4440       10.2114   296       12.0000        9.5036  2337859   20.80%
  12319  4515       10.7618   285       12.0000        9.5036  2380972   20.80%
  12719  4558       10.7662   397       12.0000        9.5206  2409816   20.66%
Elapsed time = 96.72 sec. (51095.21 ticks, tree = 72.08 MB, solutions = 5)
  13216  4699       10.8668   303       12.0000        9.5534  2452498   20.39%
  13670  5150        cutoff             12.0000        9.5676  2590353   20.27%
  13897  5178       10.4353   347       12.0000        9.5875  2604432   20.10%
  14057  5392        9.8034   273       12.0000        9.6077  2718810   19.94%
  14448  5430        cutoff             12.0000        9.6286  2748605   19.76%
  15003  5645       10.8860   391       12.0000        9.6397  2821621   19.67%
  15446  5751       10.1852   222       12.0000        9.6480  2860475   19.60%
  15918  5959       11.0000   282       12.0000        9.6614  2943024   19.49%
  16407  6273       10.9406   173       12.0000        9.6680  3044624   19.43%
  17042  6524       10.3652   311       12.0000        9.6988  3110733   19.18%
Elapsed time = 121.78 sec. (60662.49 ticks, tree = 120.26 MB, solutions = 5)
.......

Could you, please, check if CPLEX consider my initial solution is implemented algorithmically when objective is nonconvex and quadratic? Maybe a bug exists.

Here’s where we pass the initial solutions: CPLEX.jl/MOI_wrapper.jl at 71e64ea78fea5dc83c556ee76df2b5ea1df124e4 · jump-dev/CPLEX.jl · GitHub

This may be a bug in CPLEX. Can you reproduce it in a simple model?

Note that CPLEX.jl is maintained by the community; it is not officially supported by IBM. ILOG CPLEX Optimization Studio - ILOG CPLEX Optimization Studio | IBM

You could also try Gurobi. They can find optimal solutions to non-convex MIQP problems.

Thank you. I will provide an example to you see tomorrow, in tthe same post

The licence of Gurobi is free?

Where can I find it?

I looked into this in more detail. Because you set a time limit, CPLEX times out before it gets to a point in the solve where it wants to consider your MIP start. I guess it decides that it doesn’t need to return you the feasible point you provided, because you already knew it was an integer feasible point.

Play around with this code, and you will see it either terminate immediately if the myopic start value is optimal, or try quite hard to find an improvement. I suggest you also run your code for a much longer time limit and see what happens.

using JuMP, CPLEX
function main(N = 100)
    model = Model(CPLEX.Optimizer)
    set_optimizer_attribute(model, "CPXPARAM_OptimalityTarget", 3)

    a = rand(N)
    b = 0.4 * N
    @variable(model, 0 <= x[1:N] <= 2, Int)
    @constraint(model, a' * x <= b)
    @objective(model, Max, sum(x' * x))

    x0 = zeros(N)
    lhs_collected = 0.0
    for i in sortperm(a)
        lhs_collected += 2 * a[i]
        if lhs_collected > b
            break
        end
        x0[i] = 2
    end
    set_start_value.(x, x0)
    @show sum(x0' * x0)
    optimize!(model)
    return model
end
model = main(10)

Gurobi requires a paid license, but it is free for academics: GitHub - jump-dev/Gurobi.jl: Julia interface for Gurobi Optimizer

1 Like

Look my example to you better underestand. I suggest you execute these commands to see my question.

Consider my original model below (including set data):

using Random,JuMP, CPLEX
rng = MersenneTwister(2020);
###############################################################################
#Data set
m=9;
k=50;
F=4;
L=rand(10:40,k,1);
P=rand(180:200,k,1);
r=[3 3 2 2];
ℓ = L./sum(L)
Di = round.(50 .+ 20*randn(rng,F,F))
for i in 1:F
    Di[i,i] = 0
end
Dist=0.5*(Di+Di')
Dist0 = round.(80 .+ 20*randn(rng,F,1))
Lmax = 150 
Φ = 30 
H = 1:Φ 
C = round.(20 .+ 5*randn(rng,Φ)) 
###############################################################################

modelo = Model(CPLEX.Optimizer)
set_optimizer_attribute(modelo,"CPX_PARAM_TILIM",20)
set_optimizer_attribute(modelo,"CPX_PARAM_EPGAP",0.01)
set_optimizer_attribute(modelo,"CPX_PARAM_SCRIND",1)
set_optimizer_attribute(modelo, "CPXPARAM_OptimalityTarget",3)
set_optimizer_attribute(modelo,"CPX_PARAM_QTOLININD",1)

@variable(modelo, x[i in 1:m,j in 1:k], Bin)
@variable(modelo, y[i in 1:m,j in 1:F], Bin)
@variable(modelo, z[i in 1:m,j in 1:k,h in H], Bin)
@variable(modelo, w[i in 1:m,f in 1:F,h in H] >=0)
@variable(modelo, v1[i in 1:m-1,f in 1:F,h in H]>=0)
@variable(modelo, v2[i in 1:m-1,f in 1:F,h in H]>=0)
@variable(modelo, 0 <= e[j in 1:k] <= minimum(C)*Lmax - 100)

#A linear objective
@objective(modelo,Min, sum(y[i,f] for i in 1:m for f in 1:F) )

@constraint(modelo,[i in 1:m,f in 1:F,h in H],
sum(z[i,j,h] for j in sum(r[1:f-1])+1:sum(r[1:f])) == w[i,f,h]
)

@constraint(modelo,[i in 1:m-1, f in 1:F, h in H],
v1[i,f,h] - v2[i,f,h] == w[i+1,f,h] - w[i,f,h]
)

@constraint(modelo,[i in 1:m,h in H],
sum(z[i,j,h] for j in 1:k) <= 1
)

@constraint(modelo,[i in 1:m,h in H,f in 1:F,j in sum(r[1:f-1])+1:sum(r[1:f])],
z[i,j,h] <= y[i,f]
)
@constraint(modelo,[i in 1:m,j in 1:k],
x[i,j] <= sum(z[i,j,h] for h in H)
)
@constraint(modelo,[i in 1:m,j in 1:k],
sum(z[i,j,h] for h in H) <= Φ*x[i,j]
)

#As máquinas devem fazer a colheita de toda a cana em cada fazenda
@constraint(modelo,[j in 1:k],
sum(Lmax*C[h]*z[i,j,h] for i in 1:m for h in H) - e[j] == P[j]*L[j]
)

status=optimize!(modelo)

which has a linear objective. Then, I will save the final solution:

X = value.(x)
Y = value.(y)
Z = value.(z)
W = value.(w)
E = value.(e)
V1 = value.(v1)
V2 = value.(v2)

My intention is use the final solution of the previous problem as initial solution of the same problem except the objetive function (quadratic). Then, I set:

@objective(modelo,Min,sum(Dist0[f]*(v1[i,f,h]+v2[i,f,h]) for i in 1:m-1 for f in 1:F for h in H) +
sum(Dist0[f]*(w[1,f,h] + w[m,f,h]) for f in 1:F for h in H) +
sum(Dist[f,fbar]*w[i+1,f,h]*w[i,fbar,h] for i in 1:m-1 for f in 1:F for fbar in setdiff(1:F,f) for h in H) 
)

#Without mipstart
status=optimize!(modelo)

and resolve the quadratic problem (without mipstart). The CPLEX report is the following:

Version identifier: 20.1.0.0 | 2020-11-10 | 9bedb6d68
CPXPARAM_OptimalityTarget                        3   
CPXPARAM_Preprocessing_QToLin                    1   
CPXPARAM_TimeLimit                               20
CPXPARAM_MIP_Tolerances_MIPGap                   0.01
4 of 4 MIP starts provided solutions.
MIP start 'm1' defined initial solution with objective 1076.0000.
Tried aggregator 1 time.
MIQP Presolve eliminated 2751 rows and 552 columns.
MIQP Presolve modified 3483 coefficients.
Reduced MIQP has 3209 rows, 19220 columns, and 60494 nonzeros.
Reduced MIQP has 13509 binaries, 0 generals, 0 SOSs, and 0 indicators.
Reduced MIQP objective Q matrix has 5472 nonzeros.
Presolve time = 0.05 sec. (73.00 ticks)
Probing fixed 126 vars, tightened 1017 bounds.
Probing time = 0.06 sec. (71.20 ticks)
Clique table members: 1542.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 12 threads.
Root relaxation solution time = 0.14 sec. (148.94 ticks)

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

      0     0      182.7330   323                    182.7330       15
      0     0      274.2959   635                   Cuts: 417    15282
      0     0      274.3084   678                   Cuts: 289    18452
      0     0      274.3084   677                   Cuts: 265    23063
*     0+    0                         1914.0000      274.3084            85.67%
*     0+    0                         1914.0000      274.3084            85.67%
*     0+    0                         1076.0000      274.3084            74.51%
      0     0      274.3084   697     1076.0000     Cuts: 231    26875   74.51%
      0     2      274.3084   697     1076.0000      274.3084    26875   74.51%
Elapsed time = 13.89 sec. (14904.47 ticks, tree = 0.02 MB, solutions = 3)
      1     3      275.1115   742     1076.0000      274.3591    27123   74.50%
      2     4      275.2371   656     1076.0000      274.3591    28348   74.50%
      3     5      275.2371   633     1076.0000      274.3591    28376   74.50%

GUB cover cuts applied:  77
Clique cuts applied:  1
Cover cuts applied:  153
Implied bound cuts applied:  110
Flow cuts applied:  17
Mixed integer rounding cuts applied:  117
Zero-half cuts applied:  61
Gomory fractional cuts applied:  17

Root node processing (before b&c):
  Real time             =   13.50 sec. (14513.16 ticks)
Parallel b&c, 12 threads:
  Real time             =    6.58 sec. (5491.88 ticks)
  Sync time (average)   =    4.87 sec.
  Wait time (average)   =    0.00 sec.
                          ------------
Total (root+branch&cut) =   20.08 sec. (20005.04 ticks)

As you can observe, in the line 1 of the report, column Best Integeris blanck, as expected.

Now, I will provide (try) a feasible integer solution to my problem. Then, I run:

#With mipstart
for i in 1:m,j in 1:k
    set_start_value(x[i,j],X[i,j])
end
for j in 1:k
    set_start_value(e[j],E[j])
end
for i in 1:m,f in 1:F
    set_start_value(y[i,f],Y[i,f])
end
for i in 1:m,j in 1:k,h in H
    set_start_value(z[i,j,h],Z[i,j,h])
end
for i in 1:m,f in 1:F,h in H
    set_start_value(w[i,f,h],W[i,f,h])
end
for i in 1:m-1,f in 1:F,h in H
    set_start_value(v1[i,f,h],V1[i,f,h])
    set_start_value(v2[i,f,h],V2[i,f,h])
end

and run the objective again and in the sequence, optimize:

#A Quadratic objective
@objective(modelo,Min,sum(Dist0[f]*(v1[i,f,h]+v2[i,f,h]) for i in 1:m-1 for f in 1:F for h in H) +
sum(Dist0[f]*(w[1,f,h] + w[m,f,h]) for f in 1:F for h in H) +
sum(Dist[f,fbar]*w[i+1,f,h]*w[i,fbar,h] for i in 1:m-1 for f in 1:F for fbar in setdiff(1:F,f) for h in H) 
)

status=optimize!(modelo)

I will expeted a CPLEX report different of previous one. But are the same:

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

      0     0      182.7330   323                    182.7330       15
      0     0      274.2959   635                   Cuts: 417    15282
      0     0      274.3084   678                   Cuts: 289    18452
      0     0      274.3084   677                   Cuts: 265    23063
*     0+    0                         1914.0000      274.3084            85.67%
*     0+    0                         1914.0000      274.3084            85.67%
*     0+    0                         1076.0000      274.3084            74.51%
      0     0      274.3084   697     1076.0000     Cuts: 231    26875   74.51%
      0     2      274.3084   697     1076.0000      274.3084    26875   74.51%
Elapsed time = 13.98 sec. (14904.47 ticks, tree = 0.02 MB, solutions = 3)
      1     3      275.1115   742     1076.0000      274.3591    27123   74.50%
      2     4      275.2371   656     1076.0000      274.3591    28348   74.50%

e.g, my feasible solution was ignored (see the column Best integer, line 1, is blanck).

However, when I execute he code:


#With mipstart
for i in 1:m,j in 1:k
    set_start_value(x[i,j],X[i,j])
end
for j in 1:k
    set_start_value(e[j],E[j])
end
for i in 1:m,f in 1:F
    set_start_value(y[i,f],Y[i,f])
end
for i in 1:m,j in 1:k,h in H
    set_start_value(z[i,j,h],Z[i,j,h])
end
for i in 1:m,f in 1:F,h in H
    set_start_value(w[i,f,h],W[i,f,h])
end
for i in 1:m-1,f in 1:F,h in H
    set_start_value(v1[i,f,h],V1[i,f,h])
    set_start_value(v2[i,f,h],V2[i,f,h])
end

status=optimize!(modelo)

a new CPLEX report is generated, where my feasible solution is considered. Look:

Version identifier: 20.1.0.0 | 2020-11-10 | 9bedb6d68
CPXPARAM_OptimalityTarget                        3
CPXPARAM_Preprocessing_QToLin                    1
CPXPARAM_TimeLimit                               20
CPXPARAM_MIP_Tolerances_MIPGap                   0.01
1 of 4 MIP starts provided solutions.
MIP start 'm4' defined initial solution with objective 2816.0000.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 12 threads.

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

      4     6      384.4832   610     1076.0000      274.3591    46194   74.50%
Elapsed time = 4.34 sec. (1639.71 ticks, tree = 0.07 MB, solutions = 3)
      8     7      401.9055   673     1076.0000      274.3591    46468   74.50%
     12     9      275.2371   680     1076.0000      275.1115    46735   74.43%

GUB cover cuts applied:  77
Clique cuts applied:  1
Cover cuts applied:  153
Implied bound cuts applied:  154
Flow cuts applied:  17
Mixed integer rounding cuts applied:  117
Zero-half cuts applied:  61
Gomory fractional cuts applied:  17

Root node processing (before b&c):
  Real time             =    0.06 sec. (43.33 ticks)
Parallel b&c, 12 threads:
  Real time             =   20.08 sec. (4300.27 ticks)
  Sync time (average)   =    5.53 sec.
  Wait time (average)   =    0.00 sec.
                          ------------
Total (root+branch&cut) =   20.14 sec. (4343.61 ticks)

My question is: Why, when I run the objective quadratic function again, CPLEX looses my feasible solution??

Please, help me to underestand this.

Note: when the objetive is linear, nothing wrong happens.

Maybe you comment is correct. Because, when I use mipstart and the column of Best integer is totaly blank, the algorithm provide the initial solution (it seems that the bad quality solution does not enter into the branch-and-bound, differently when the objetive is linear and convex).

There might be a couple of things going on here:

  1. We are retaining the previous MIP starts, rather than clearing them on subsequent solves: Delete previous MIP starts if they exist · Issue #352 · jump-dev/CPLEX.jl · GitHub. I don’t know what, if anything, fixing this will do.
  2. Warm-starts are a hint to the solver. They may choose to use them if they think it will help. It certainly looks like in this case the solver receives them at the start of the log, but ignores them during the branch-and-bound. I don’t know why. You would have to follow up with CPLEX support for more details.

p.s. As a quick way, you can set the start value from a previous solution as follows:

vars = all_variables(model)
set_start_value.(vars, value.(vars))
1 Like

Perfect.
Thank you so much for the hints.
I have installed Gurobi and the results (either linear and quadratic model) were to better.
Moreover, the algorithm implemented by GUROBI without mipstart determined feasible solutions.

1 Like