Output of MPEC methods wrong under certain conditions

[reinderien]

I'm fairly sure this is an APMonitor bug and not a Gekko bug. This script represents the minimal case where I have been able to reproduce the problem:

from gekko import GEKKO

MIDPOINT = 0.0103

m = GEKKO(remote=False)
try:
    m.options.solver = 1
    m.solver_options = [
        'minlp_as_nlp 0',
        'minlp_gap_tol 1e-3',
        'minlp_integer_tol 1e-3',
    ]

    zero = m.Param(0, name='zero')
    shards_cont = m.Param(-1.867, name=f'shards_cont')
    building = m.Var(name=f'buildings', integer=True)
    shards_each = m.Var(name=f'shards_each', integer=True)
    power = m.Intermediate(building**-0.6 * 103e6, name=f'power')
    shards_pos = m.max2(zero, shards_cont)

    m.Equations((
        shards_each > shards_pos,
        shards_each < shards_pos + (1 - MIDPOINT),
        building <= 50,
    ))
    m.Minimize(power)
    m.solve()

    (cont,), (pos,), (each,) = shards_cont, shards_pos, shards_each
    if round(each) != 0:
        print(f'The shard solution is wrong:\n'
              f'max(0, {cont}) = {pos} == {MIDPOINT} != 0\n'
              f'which produces an incorrect shard count of {each}')
finally:
    m.cleanup()

The corresponding APM script is

Model
Parameters
    zero = 0
    shards_cont = -1.867
End Parameters
Variables
    int_buildings = 0
    int_shards_each = 0
    v3 = 0
End Variables
Intermediates
    power=((((int_buildings)^(-0.6)))*(103000000.0))
End Intermediates
Equations
    int_shards_each>v3
    int_shards_each<(v3+0.9897)
    int_buildings<=50
    minimize power
End Equations
Connections
    zero = max2_1.x[1]
    shards_cont = max2_1.x[2]
    v3 = max2_1.y
End Connections
Objects
    max2_1 = max
End Objects

End Model

The output is

The shard solution is wrong:
max(0, -1.867) = 0.0103 == 0.0103 != 0
which produces an incorrect shard count of 1.0

APM version is 0.9.2 with Gekko 0.2.8 via pip.

I don't know why, but under certain numerical and model conditions (including the above, 100% of the time), the output of max2 is just wrong. It prefers to output what, in this case, is MIDPOINT. It seems that the optimizer is fighting to respect the upper bound - int_shards_each<(v3+0.9897) - and in doing so, it fails to respect the output of max.

[reinderien]

The issue seems not to be limited to max2. When I try a workaround that replaces max2(0, x) with (1 + sign2(x))/2, I see similar edge-case numerical errors (for different numerical conditions than the one presented in the repro above).

[APMonitor]

MPEC solutions aren't very reliable and can get stuck at the switching point as you observed. I recommend using the max3() function instead that uses binary switching variables. The problem is solved correctly using max3().

from gekko import GEKKO

MIDPOINT = 0.0103

m = GEKKO(remote=False)
try:
    m.options.solver = 1
    m.solver_options = [
        'minlp_as_nlp 0',
        'minlp_gap_tol 1e-3',
        'minlp_integer_tol 1e-3',
    ]

    zero = m.Param(0, name='zero')
    shards_cont = m.Param(-1.867, name=f'shards_cont')
    building = m.Var(name=f'buildings', integer=True)
    shards_each = m.Var(name=f'shards_each', integer=True)
    power = m.Intermediate(building**-0.6 * 103e6, name=f'power')
    shards_pos = m.max3(zero, shards_cont)

    m.Equations((
        shards_each > shards_pos,
        shards_each < shards_pos + (1 - MIDPOINT),
        building <= 50,
    ))
    m.Minimize(power)
    m.solve()

    (cont,), (pos,), (each,) = shards_cont, shards_pos, shards_each
    if round(each) != 0:
        print(f'The shard solution is wrong:\n'
              f'max(0, {cont}) = {pos} == {MIDPOINT} != 0\n'
              f'which produces an incorrect shard count of {each}')
finally:
    m.cleanup()

 ----------------------------------------------------------------
 APMonitor, Version 0.9.2
 APMonitor Optimization Suite
 ----------------------------------------------------------------

 --------- APM Model Size ------------
 Each time step contains
   Objects      :  0
   Constants    :  0
   Variables    :  11
   Intermediates:  1
   Connections  :  0
   Equations    :  8
   Residuals    :  7

 Number of state variables:    9
 Number of total equations: -  6
 Number of slack variables: -  5
 ---------------------------------------
 Degrees of freedom       :    -2

 * Warning: DOF <= 0
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.00 NLPi:    2 Dpth:    0 Lvs:    2 Obj:  9.85E+06 Gap:       NaN
Iter:     2 I: -1 Tm:      0.00 NLPi:    1 Dpth:    1 Lvs:    1 Obj:  9.85E+06 Gap:       NaN
--Integer Solution:   9.85E+06 Lowest Leaf:   9.85E+06 Gap:   0.00E+00
Iter:     3 I:  0 Tm:     -0.00 NLPi:    2 Dpth:    1 Lvs:    1 Obj:  9.85E+06 Gap:  0.00E+00
 Successful solution

 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  0.031 sec
 Objective      :  9850430.747837381
 Successful solution
 ---------------------------------------------------

[reinderien]

It seems like max3 does not run into the same issues, but encounters different issues: if I combine this model with a secondary optimization objective, max3 executes much more quickly and accurately but actively hinders the optimality of the other objective - in my case, by a large margin of 56%. So it seems there's no silver bullet.

[APMonitor]

Try setting the solver options to give you a better solution with:

m.solver_options = ['minlp_gap_tol 1.0e-4',\
                    'minlp_maximum_iterations 10000',\
                    'minlp_max_iter_with_int_sol 500']
m.options.solver = 1

You can't necessarily compare to the max2() solution because it was violating constraints. You may also want to try solving with m.options.SOLVER=3 initially to get a better initial guess for the MINLP solution when you switch back to m.options.SOLVER=1.

[reinderien]

Interesting stuff. The solver_options shown did not help at all, but doing an interior-point pre-solve "mostly" did the job, and switching to breadth-first greatly helped find a good solution faster.