Suboptimal solution for multilinear relations

[ThibautCoudroy]

Hello,

I'm using SCIP to solve a continuous multilinear problem. I have two equal formulations for the objective function but I found a case for which one formulation gives a far suboptimal solution. Yet, SCIP claims an optimal result for both formulations.

Here is the PySCIPOpt code. The main function is called twice, where the first execution/formulation gives the suboptimal result while the second execution/formulation gives the optimal one.

from pyscipopt import Model

RATIO, FLOW = 0, 1

def solve_problem(obj_func=RATIO):
    """Solve reconciliation problem for a given objective function
    The expected result is to have both formulations leading to the same solution,
    but here the RATIO one is significantly worse
    """

    model = Model()

    flows = [model.addVar(f"f{i}") for i in range(5)]
    model.addCons(sum(flows) == 300)

    weights = [model.addVar(f"s{i}_w", ub=1) for i in range(5)]
    for i in range(5):
        model.addCons(flows[i] * (weights[i]+1) == 100)

    # Two equivalent definitions for the objective function (since the constraints above are equalities)
    obj = model.addVar("obj")
    if obj_func == RATIO:
        model.addCons(obj >= sum([abs(0.9 - weights[i]) for i in range(5)]))
    else:
        model.addCons(obj >= sum([abs(0.9 - (100/flows[i] - 1)) for i in range(5)]))

    model.setObjective(obj)
    model.optimize()
    print(model.getBestSol())


if __name__ == "__main__":
    # Suboptimal minimization (obj=0.84)
    solve_problem(obj_func=RATIO)
    # Optimal minimization (obj=0.78)
    solve_problem(obj_func=FLOW)

First call result:
{'f0': 53.234780356212, 'f1': 89.7560785882046, 'f2': 52.258562416300194, 'f3': 52.49201622298299, 'f4': 52.258562416300194, 's0_w': 0.8784711673621274, 's1_w': 0.11413067027880895, 's2_w': 0.9135620131182284, 's3_w': 0.9050516106620469, 's4_w': 0.913562016031791, 'obj': 0.83957379917113}

Second call result:
{'f0': 52.63158306654943, 'f1': 52.631578975054545, 'f2': 89.4736555021569, 'f3': 52.63160348118459, 'f4': 52.631578975054545, 's0_w': 0.8999998512981086, 's1_w': 0.899999995642411, 's2_w': 0.11764741742991985, 's3_w': 0.8999991143293624, 's4_w': 0.899999999000531, 'obj': 0.782352681977184}

Expectation: Both objective values should be equal.

Note that there is a symmetry between the variables, so they should not be compared one to one.

I ran the same code on various situations (changing the number of flows, the value of the sum of flows and value used in the objective function), but this particular choice of values is the one for which the gap is the highest I have found. The second formulation (FLOW) always had a better result than the other (RATIO).

Some things I have tried:

• Removing the upper bound of 1 on the weights solves the issue by providing the optimal result (that has every weight under 1) way faster than the other formulation. Changing it to any value greater than 1 also solves the issue.
• Changing the bounds of weights from [0, 1] to [1, 2] and adapting consequently the constraints and objective function seem to fix the issue, but I'm unsure if it is luck or if it is handled differently

Is it expected that the results could be this different by substituting one element of the objective function with another variable linked by an equality?

[DominikKamp]

It is not expected. The first formulation

model.cip.zip

with feasible solution

model.sol.zip

of value +7.82352941163668e-01 results in

SCIP Status        : problem is solved [optimal solution found]
Solving Time (sec) : 0.88
Solving Nodes      : 31
Primal Bound       : +8.39573799171130e-01 (3 solutions)
Dual Bound         : +8.39573799171130e-01
Gap                : 0.00 %

with

SCIP version 9.2.4 [precision: 8 byte] [memory: block] [mode: optimized] [LP solver: SoPlex 7.1.6] [GitHash: 8f58ab14bc]
Copyright (c) 2002-2025 Zuse Institute Berlin (ZIB)

External libraries: 
  Readline EditLine w  GNU library for command line editing (gnu.org/s/readline)
  SoPlex 7.1.6         Linear programming solver developed at Zuse Institute Berlin (soplex.zib.de) [GitHash: 08bb3155]
  CppAD 20180000.0     Algorithmic Differentiation of C++ algorithms developed by B. Bell (github.com/coin-or/CppAD)
  ZLIB 1.2.12          General purpose compression library by J. Gailly and M. Adler (zlib.net)
  GMP 6.3.0            GNU Multiple Precision Arithmetic Library developed by T. Granlund (gmplib.org)
  ZIMPL 3.7.0          Zuse Institute Mathematical Programming Language developed by T. Koch (zimpl.zib.de)
  AMPL/MP 690e9e7      AMPL .nl file reader library (github.com/ampl/mp)
  PaPILO 2.4.4         parallel presolve for integer and linear optimization (github.com/scipopt/papilo) (built with TBB) [GitHash: d50cb8a6]
  Nauty 2.8.8          Computing Graph Automorphism Groups by Brendan D. McKay (users.cecs.anu.edu.au/~bdm/nauty)
  sassy 1.1            Symmetry preprocessor by Markus Anders (github.com/markusa4/sassy)
  Ipopt 3.14.17        Interior Point Optimizer developed by A. Waechter et.al. (github.com/coin-or/Ipopt)

@svigerske Would you like to have a look?

[svigerske]

Thanks for the cip and sol files. I can reproduce if I switch to SoPlex.
A workaround is to set nlhdlr/convex/enabled = FALSE. I could image that for obj_func=FLOW, nlhdlr_convex already did not get active here.

The model may also run better if the abs is avoided, and you instead have constraints obj[i] >= 0.9 - weights[i] and obj[i] >= weights[i] - 0.9 with an objective function that sums the obj[i].

[svigerske]

The problem is not with nlhdlr_convex. I can also produce a wrong optimal value when turning it off and turning all primal heuristics off.

I couldn't reproduce an issue with CPLEX, but running with SoPlex and LPSCHECK=true doesn't point out a mistake by SoPlex.

With a debug solution, I get something like debugging solution was cut off in local node #32 (0x63b4f80943a8) at depth 9, so no boundchange or cut seems to be responsible.

The only way I can get it work reliably is by turning off symmetry handling. In fact, it seems that the "cut off" of a local node is triggered by symmetry_orbitopal.c,

[symmetry_orbitopal.c:2188] debug: Found 0 reductions during orbitopal reduction for orbitope 0
Col hash                               54; Row hash                                0; Hash                               54
  node 2: (1, 0, 0) 
  node 1: (0, 0, 0) 
[symmetry_orbitopal.c:1656] debug: Cannot repair infeasibility for column 0 (original: 0), min
[symmetry_orbitopal.c:2188] debug: Found 0 reductions during orbitopal reduction for orbitope 0
[symmetry_orbitopal.c:2197] debug: Detected infeasibility during orbitopal reduction for orbitope 0

Also just commenting out

scip/src/scip/prop_symmetry.c

Line 7787 in 1639ded

*result = SCIP_CUTOFF;

makes it solve fine again. (I switch to testing with v10-minor here, but seems to be the same with SCIP 9.2.3).

The symmetries that are detected look fine, though:

Display symmetries of component 0.
Permutation 0:
  (<t_f3>,<t_f4>)
  (<t_s3_w>,<t_s4_w>)
Permutation 1:
  (<t_f2>,<t_f3>)
  (<t_s2_w>,<t_s3_w>)
Permutation 2:
  (<t_f1>,<t_f2>)
  (<t_s1_w>,<t_s2_w>)
Permutation 3:
  (<t_f0>,<t_f1>)
  (<t_s0_w>,<t_s1_w>)

I don't want to conclude that there is an issue with propagations of symmetry here. This method is supposed to cut off optimal solutions, but I may need some help from @christopherhojny to check that this cutoff is indeed correct.
I was building with LPS=spx and DEBUGSOL=true (also Ipopt and Papilo, but shouldn't make a difference here) and running https://github.com/user-attachments/files/21711564/model.cip.zip with model.sol.txt via

./bin/scip -c "set misc debugsol model3.sol read model.cip set nlhdlr conv ena F set heur emph off o"

[christopherhojny]

I could reproduce the issue and I have also analyzed all reductions made by orbitopal reduction. To me it seems that all these reductions (including declaring nodes as infeasible) are correct.

Using a debug solution is difficult, because the lexicographic order underlying the symmetry reductions is different at different nodes of the branch-and-bound tree. Making this lexicographic order the same at each node of the branch-and-bound tree by setting columnordering = SCIP_COLUMNORDERING_NONE in

scip/src/scip/prop_symmetry.c

Line 5677 in 80c1fbf

columnordering = SCIPorbitopalReductionGetDefaultColumnOrdering(propdata->orbitopalreddata);

unfortunately does not allow to debug this: the instance is then solved correctly.

In case it helps for debugging, I provide another solution model2.sol.txt. This solution is not immediately cut off when using it as a debug solution; maybe this allows to get some more insights.

[svigerske]

That other solution was the same that I had send before, but with some lines (both variable names and value together) permuted. I
think you meant to use this one:

objective value:                    0.782352941163668
f0                                   89.4736842095014 	(obj:0)
f1                                   52.6315789473684 	(obj:0)
f2                                   52.6315789473684 	(obj:0)
f3                                   52.6315789473684 	(obj:0)
f4                                   52.6315789473684 	(obj:0)
s0_w                                0.117647058836332 	(obj:0)
s1_w                                              0.9 	(obj:0)
s2_w                                              0.9 	(obj:0)
s3_w                                              0.9 	(obj:0)
s4_w                                              0.9 	(obj:0)
obj                                 0.782352941163668 	(obj:1)

But all I'm getting is again

Col hash                                0; Row hash                                0; Hash                                0
  node 5: (1, 4, 19) 
  node 2: (1, 0, 0) 
  node 1: (0, 0, 0) 
[symmetry_orbitopal.c:1547] debug: Start propagating static orbitope: 
>           1            0            4            3            2 < (IDs)
[t_f1 +50.00,+72.86 t_f0 +77.14,+100.00 t_f4 +50.00,+72.86 t_f3 +50.00,+72.86 t_f2 +50.00,+72.86 ] (row 0)
[symmetry_orbitopal.c:1671] debug: Cannot repair infeasibility for column 0 (original: 1), min
[symmetry_orbitopal.c:2212] debug: Found 0 reductions during orbitopal reduction for orbitope 0
[symmetry_orbitopal.c:2221] debug: Detected infeasibility during orbitopal reduction for orbitope 0
[solve.c:558] debug:  -> propagator <symmetry> detected cutoff
[solve.c:642] debug:  --> domain propagation of node 0x55555c714290 finished: cutoff!
[solve.c:646] debug: solinsubtree: 1
[solve.c:4763] debug: node solving iteration 1 finished: cutoff=1, postpone=0, propagateagain=0, solverelaxagain=0, solvelpagain=0, nlperrors=0, restart=0
[solve.c:319] debug: calling primal heuristics in depth 2 (timing: 8)
  0.1s|     4 |     1 |   127 |  38.3 |  1107k |   2 |   - |  32 |   7 |   3 |   0 |   0 | 5.712500e-01 |      --      |    Inf 
[solve.c:5399] debug: Processing of node 4 in depth 2 finished. 1 siblings, 0 children, 0 leaves left
[solve.c:5401] debug: **********************************************************************
[debug.c:1153] ERROR: debugging solution was cut off in local node #5 (0x55555c714290) at depth 2
[debug.c:1046] ERROR: invalid global upper bound: <t_f0>[89.4736842095014] <= 77.1406249999467

I'll put this aside, as I currently don't know how to debug this further with reasonable effort.

If orbitopal symmetry propagation can conclude that a node of the tree is infeasible, then it should know how the debug solution needs to be adapted to be feasible in another node of the tree. It would be great if it were to adjust the debug solution accordingly then.