Monday, January 26, 2015

The Hostile Brothers Problem (3)

Follow up on http://yetanothermathprogrammingconsultant.blogspot.com/2015/01/the-hostile-brothers-problem-1.html and http://yetanothermathprogrammingconsultant.blogspot.com/2015/01/the-hostile-brothers-problem-2.html.

Now some experiments that did not work.

Find global solutions with Baron

This is a very difficult job. One idea that may help is to add some symmetry breaking constraints:

image

This can reduce the search space considerably. We run the equation over all i–1: in GAMS this is notation that means we drop the first possible equation (the equation related to ord(i)=1).

In addition we can pass on a good solution we found with techniques described earlier. Unfortunately this helps but not enough to bring us close to proving global optimality:

  Iteration    Open nodes       Total time    Lower bound      Upper bound
          1             1        000:00:01     0.100000E-08     2.00000   
          1             1        000:00:07     0.200000E-08     2.00000   
*         2             2        000:00:10     0.163965E-01     2.00000   
*         2             2        000:00:18     0.173057E-01     2.00000   
          6             4        000:00:51     0.173057E-01     1.67193   
         12             8        000:01:29     0.173057E-01     1.56665   
         17            12        000:02:31     0.173057E-01     1.53265   
         21            15        000:03:09     0.173057E-01     1.50781   
         25            18        000:03:58     0.173057E-01     1.49840   
         29            20        000:04:29     0.173057E-01     1.49070   
         33            22        000:05:23     0.173057E-01     1.49070
   

Default. Note that Baron actually finds good solutions quickly. But after a little while the lower bound will not move again and the upper bound will move slower and slower.

  Iteration    Open nodes       Total time    Lower bound      Upper bound
          1             1        000:00:02     0.363949E-03     2.00000   
*         1             1        000:00:24     0.111980E-02     1.01351   
*         1             1        000:00:29     0.139747E-01     1.01351   
*         1             1        000:00:43     0.151056E-01     1.01351   
          1             1        000:00:49     0.151056E-01     1.01351   
          3             2        000:01:23     0.151056E-01     1.01351   
          4             3        000:02:00     0.151056E-01     1.01351   
          6             4        000:02:42     0.151056E-01     1.00880   
          8             5        000:03:29     0.151056E-01     1.00802   
         11             7        000:04:17     0.151056E-01     1.00802   
         12             8        000:04:50     0.151056E-01     1.00802   
         15             8        000:05:35     0.151056E-01     1.00191
   

Here we added symmetry breaking constraints. This helps with the upper bound.

Iteration    Open nodes       Total time    Lower bound      Upper bound
        1             1        000:00:04     0.175061E-01     2.00000   
        1             1        000:00:19     0.175061E-01     1.01351   
        4             3        000:00:52     0.175061E-01     1.01351   
        6             4        000:01:30     0.175061E-01     1.00930   
        8             6        000:02:13     0.175061E-01     1.00930   
       10             7        000:02:55     0.175061E-01     1.00783   
       11             7        000:03:26     0.175061E-01     1.00783   
       13             7        000:04:01     0.175061E-01     1.00783   
       15             8        000:04:39     0.175061E-01     1.00244   
       17            10        000:05:13     0.175061E-01     1.00000
   

Here we use a starting point and symmetry breakers. This gives best performance but still not sufficient to get close to a global solution.

Understandably global solvers can not handle just anything we throw at them. Heaving said this I think Baron does a good job finding good solutions if we don’t provide our own initial points.

Other Global solvers

None of the global solvers can solve this problem quickly. It is interesting to see which objective can be used with different global solvers.

image[5]

Solver Discontinuous Formulated allowed Smooth Formulation Accepted
Baron Can not handle function ‘min’ OK
LocalSolver OK OK (n=75 problem ran out of demo license)
OQNLP OK OK
LindoGlobal OK automatic reformulated into bigM constraints OK (n=75 problem: Model size exceeds LindoGlobal limits of (2000,3000)).
LGO OK OK (n=75 problem:Reduced Model Size exceeds LGO limits of (2000,3000)).

 

Can we solve as a MIP?

The quadratic expressions form a problem. Although Cplex can solve non-convex QPs it does not handle non-convex quadratic constraints (I believe SCIP does). We can use some piecewise linear approximation, but it seems to make sense to first try a linear version with a different distance function: the Manhattan distance. Here the distance between two points is the sum of the distance in the x direction and in the y direction. The formulation can look like:   

image

This model turns out not easy to solve. Even after 2 hours compute time the resulting locations are not very good:

image

Note: the plot is made with R and SQLite:

image