Multi-objective MINLP problem in JuMP

Running HiGHS 1.5.3 [date: 1970-01-01, git hash: 45a127b78]
Copyright (c) 2023 HiGHS under MIT licence terms
Hessian has dimension 2558 but no nonzeros, so is ignored
Presolving model
2552 rows, 1825 cols, 8804 nonzeros
2088 rows, 1479 cols, 8450 nonzeros
1854 rows, 1417 cols, 7778 nonzeros

Solving MIP model with:
   1854 rows
   1417 cols (525 binary, 0 integer, 0 implied int., 892 continuous)
   7778 nonzeros

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   0.9981586674    -inf                 inf        0      0      0         0     0.1s
         0       0         0   0.00%   0.9981586674    -inf                 inf        0      0      0       704     0.1s
 L       0       0         0   0.00%   0.9981586674    0.9966088199       0.16%     3415    389    174      1361     1.3s
 L       0       0         0   0.00%   0.9981586674    0.9981586674       0.00%     3439    394    190      1741     1.4s

Solving report
  Status            Optimal
  Primal bound      0.998158667416
  Dual bound        0.998158667416
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    0.998158667416 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            1.45 (total)
                    0.09 (presolve)
                    0.00 (postsolve)
  Nodes             1
  LP iterations     1889 (total)
                    0 (strong br.)
                    663 (separation)
                    509 (heuristics)
ERROR:   Cannot solve MIQP problems with HiGHS
ERROR:   Cannot solve MIQP problems with HiGHS
Hessian has dimension 2558 but no nonzeros, so is ignored
Presolving model
2552 rows, 1825 cols, 8804 nonzeros
2088 rows, 1479 cols, 8450 nonzeros
1854 rows, 1417 cols, 7778 nonzeros

Solving MIP model with:
   1854 rows
   1417 cols (525 binary, 0 integer, 0 implied int., 892 continuous)
   7778 nonzeros

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   0.9981586674    -inf                 inf        0      0      0         0     0.1s
         0       0         0   0.00%   0.9981586674    -inf                 inf        0      0      0       704     0.1s
 L       0       0         0   0.00%   0.9981586674    0.9966088199       0.16%     3415    389    174      1361     1.6s
 L       0       0         0   0.00%   0.9981586674    0.9981586674       0.00%     3439    394    190      1741     1.7s

Solving report
  Status            Optimal
  Primal bound      0.998158667416
  Dual bound        0.998158667416
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    0.998158667416 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            1.74 (total)
                    0.08 (presolve)
                    0.00 (postsolve)
  Nodes             1
  LP iterations     1889 (total)
                    0 (strong br.)
                    663 (separation)
                    509 (heuristics)
* Solver : MOA[algorithm=MultiObjectiveAlgorithms.Dichotomy, optimizer=HiGHS]

* Status
  Result count       : 0
  Termination status : OTHER_ERROR
  Message from the solver:
  "Solve complete. Found 0 solution(s)"

* Candidate solution (result #1)
  Primal status      : NO_SOLUTION
  Dual status        : NO_SOLUTION
  Objective bound    : [NaN,9.98159e-01]
  Relative gap       : 0.00000e+00

* Work counters
  Solve time (sec)   : 1.45600e+00
  Simplex iterations : 1889
  Barrier iterations : -1
  Node count         : 1

but this is a version with different inputs as your rewritten code

HiGHS does not support mixed integer quadratic programs. (See the errors in the log.)

Use a solver like Gurobi instead.

Thanks, but for this problem Gurobi is very fast than HiGHS but with very far than optimal solutions (compare to the solutions of HiGHS for the single-objective problem). Can you guess anything about this?

Another point that is a subject of the question for me is if HiGHS does not support MIQP what is the results repented after running the code (that is, what problem it wanted try to solve)?

but for this problem Gurobi is very fast than HiGHS but with very far than optimal solutions (compare to the solutions of HiGHS for the single-objective problem)

To re-iterate, HiGHS does not support MIQP. It will not find a solution.

The Gurobi solution will be optimal. If HiGHS is “better” it is because that solution isn’t feasible.

is if HiGHS does not support MIQP what is the results repented after running the code

Nothing. Per the log

* Status
  Result count       : 0
  Termination status : OTHER_ERROR
  Message from the solver:
  "Solve complete. Found 0 solution(s)"

Thanks again. I misexplained my mean! I added the quadratic term as the second objective to my model. However, fastness of Gurobi and best results of HiGHS (non-optimal results of Gurobi) were related to my single-objective mixed-integer linear programming!

  • I compared the results with those generated by simulations.

Do you have an example? If Gurobi does not find an optimal solution that is a bug.

Based on your rewritten code:

using JuMP 
using HiGHS
A = 5*rand(1,365)
B = rand(1,365)
T = length(A)
M = 10^6
model = Model(); 
set_optimizer(model,HiGHS.Optimizer); # test with Gurobi
@variables(model, begin
    0 <= S
    1 <= R <= 2
    y1[1:T], Bin
    y2[1:T], Bin
    0 <= Y[t = 1:T] <= A[t]
    0 <= O[1:T]
    0 <= V[0:T]
end);
fix(V[0], 0; force = true)
@expressions(model, begin
    Q[t in 1:T], B[t] * R
end)
@constraints(model, begin
    [t in 1:T], Y[t] <= V[t-1]
    [t in 1:T], Y[t] >= A[t] - M * (1 - y1[t])
    [t in 1:T], Y[t] >= V[t-1] - M * y1[t]
    [t in 1:T], V[t] <= V[t-1] + Q[t] - O[t] - Y[t]
    [t in 1:T], V[t] <= S - Y[t]
    [t in 1:T], V[t] >= V[t-1] + Q[t] - O[t] - Y[t] - M * (1 - y2[t])
    [t in 1:T], V[t] >= S - Y[t] - M * y2[t]
end);
@objective(model, Max, sum(Y[t] for t in 1:T)/sum(A[t] for t in 1:T))
optimize!(model)
solution_summary(model)
println(objective_value(model)) 
println("S = ", value.(S), "\n","R = ",value.(R)) 

but just with my real data. It is expected S has a smaller amount and follows the simulation results, while the Gurobi’s result for this variable (with real data) is so far from what is expected!

My problem with Gurobi is because of the value of the S variable, not the value of the objective function

What is the output of the Gurobi log and HiGHS log?

with my own data or random data?
with my own data:
for Gurobi

Set parameter Username
Academic license - for non-commercial use only - expires 2024-08-24
Gurobi Optimizer version 10.0.2 build v10.0.2rc0 (win64)

CPU model: Intel(R) Core(TM) i7-8665U CPU @ 1.90GHz, instruction set [SSE2|AVX|AVX2]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 2555 rows, 1828 columns and 8817 nonzeros
Model fingerprint: 0x16e88fd6
Variable types: 1098 continuous, 730 integer (730 binary)
Coefficient statistics:
  Matrix range     [7e-02, 1e+06]
  Objective range  [6e-04, 6e-04]
  Bounds range     [1e+00, 7e+00]
  RHS range        [1e+06, 1e+06]
Found heuristic solution: objective -0.0000000
Presolve removed 1600 rows and 1032 columns
Presolve time: 0.11s
Presolved: 955 rows, 796 columns, 2798 nonzeros
Found heuristic solution: objective 0.1011473
Variable types: 497 continuous, 299 integer (299 binary)

Root relaxation: objective 9.983067e-01, 514 iterations, 0.00 seconds (0.00 work units)

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

H    0     0                       0.9983067    0.99831  0.00%     -    0s
     0     0    0.99831    0   38    0.99831    0.99831  0.00%     -    0s

Explored 1 nodes (514 simplex iterations) in 0.12 seconds (0.07 work units)
Thread count was 8 (of 8 available processors)

Solution count 3: 0.998307 0.101147 -0 
No other solutions better than 0.998307

Optimal solution found (tolerance 1.00e-04)
Best objective 9.983067403479e-01, best bound 9.983067403479e-01, gap 0.0000%

User-callback calls 4051, time in user-callback 0.00 sec
0.9983067403479037
S = 4888.936547829832
R = 2.0

and for HiGHS

Running HiGHS 1.5.3 [date: 1970-01-01, git hash: 45a127b78]
Copyright (c) 2023 HiGHS under MIT licence terms
Presolving model
2552 rows, 1825 cols, 8804 nonzeros
2140 rows, 1576 cols, 8547 nonzeros
1859 rows, 1419 cols, 7795 nonzeros

Solving MIP model with:
   1859 rows
   1419 cols (526 binary, 0 integer, 0 implied int., 893 continuous)
   7795 nonzeros

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   0.9983067403    -inf                 inf        0      0      0         0     0.1s
         0       0         0   0.00%   0.9983067403    -inf                 inf        0      0      0       662     0.1s

Solving report
  Status            Optimal
  Primal bound      0.998306740348
  Dual bound        0.998306740348
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    0.998306740348 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            0.71 (total)
                    0.10 (presolve)
                    0.00 (postsolve)
  Nodes             1
  LP iterations     1334 (total)
                    0 (strong br.)
                    302 (separation)
                    370 (heuristics)
0.9983067403479303
S = 197.0587824
R = 2.0

and as I said my point is on S variable and not objective function

If they have the same objective function value, then I assume there are multiple possible values that S can take. Why does it need to be small? This just looks like two equivalent solutions to your problem.

1 Like

As I learned in optimization courses, if we do not have an upper bound for a variable, after solving the problem to the optimality, that variable should take the lowest value that satisfies the problem (keeps it feasible). My point in this regard was that the HiGHS solution is consistent with this fact, but Gorubi’s is not

This is not true in general. Solvers may return any feasible solution with an optimal objective value.

To secondarily minimize S, try adding a small negative term to the objective like - 0.0001S.

1 Like

That’s a nice technique. However, it is fine that “any feasible solution” becomes the best one!

2 Likes

@odow is correct, and it’s what I suggested in reply to your email. You model has a non-unique solution and, without introducing a tie-breaking term into the objective as Oscar suggested, different optimization solvers will give you different optimal solutions. Indeed, if you change the order in which constraints are defined, the same solver may give a different solution, since numerics will change subtly, and algorithmic ties will be broken differently

1 Like

HiGHS reports that

Hessian has dimension 2558 but no nonzeros, so is ignored

so you’ve not defined a non-zero Hessian.

If you had defined a non-zero Hessian, HiGHS would have reported that it doesn’t solve MIQPs and returned an error.

2 Likes

Thanks for the response.
Two follow up question (but not related to Julia):

  1. How can one mathematically prove that the proposed model does not have a unique solution for S and, for example, adding a regularized term -0.00001S to the objective function leads to a unique solution for S?
  2. According to the model, do dynamic variables such as Y[t] have unique solutions?

In general it is pretty hard to prove uniqueness of the primal solution for MILP. Why do you need to?

Since my original formula does not have a unique solution for S, I want to justify (using proof) why I added a regularized term to the objective (as a support for my optimization results), considering that S is a continuous variable.