Experiments with Hostile Brothers nonconvex NLP model

In this post, let's do some experiments with the "hostile brothers problem." We have a plot of land. In my test models \([L,U]\times [L,U]\) with \(L=0\), \(U=100\)). We want to place \(n\) brothers on this plot. As they don't get along, we want to spread them out by maximizing the distance between neighbors. In modeling terms, we create a maximin model that maximizes the minimum distance between any two brothers. This yields some non-convex problems, so I'll make the problem small: we just have \(n=10\) brothers.

Model 1: the simplest quadratic model

We start with the most obvious model: 

Model 1
\[\begin{aligned}\max\>&{\color{darkred}z}\\ & {\color{darkred}z} \le \sum_c \left({\color{darkred}x}_{i,c}-{\color{darkred}x}_{j,c}\right)^2 && \forall i\lt j \\ & {\color{darkred}x}_{i,c} \in [{\color{darkblue}L},{\color{darkblue}U}] \\ & c \in \{x,y\} \\ & i,j \in \{1,\dots,n\} \end{aligned}\]

Instead of working with the distance, we use the squared distance. That saves us a square root. The square root may cause problems: it is not differentiable at 0. So it is better to use a quadratic model.

The model is very small:

MODEL STATISTICS

BLOCKS OF EQUATIONS           1     SINGLE EQUATIONS           45

BLOCKS OF VARIABLES           2     SINGLE VARIABLES           21

NON ZERO ELEMENTS           225     NON LINEAR N-Z            180

CODE LENGTH                 451     CONSTANT POOL              16

However, this is not an easy model at all! Proving global optimality is actually very difficult. I gave the solver 1,000 seconds, and it came back with:

    Elapsed time (sec): 999, estimated tree completion: 0.00015

        Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time

     2483289  1774.764542  4812.031822     11 2420538     24   63.12%      25    999

 

Many solvers now produce a progress measure indicating how much of the tree has been explored [1]. I usually am watching the gap. A disadvantage of the gap is that it moves quickly at the beginning and slowly at the end. (Not really a disadvantage: this is much better than the other way around.) The tree completion should progress more linearly with time. Here, we did not really put a dent in this problem in 1,000 seconds: tree completion is about zero, and the gap is very large. This is not an issue with this particular solver (Xpress): other solvers are struggling in the same way.

In all of the following experiments, I'll keep the time limit of 1,000 seconds. Let's see if we can get closer to proving optimality.

Model 2: less nonlinear, but not better

This is a small experiment, that is unlikely to yield great results. Let's try to make the problem less nonlinear as follows:

 

Model 2
\[\begin{aligned}\max\>&{\color{darkred}z}\\ & {\color{darkred}d}_{i,j,c} ={\color{darkred}x}_{i,c}-{\color{darkred}x}_{j,c} && \forall i\lt j, c \\ & {\color{darkred}z} \le \sum_c {\color{darkred}d}_{i,j,c}^2 && \forall i\lt j \\ & {\color{darkred}x}_{i,c} \in [{\color{darkblue}L},{\color{darkblue}U}] \\ & c \in \{x,y\} \\ & i,j \in \{1,\dots,n\} \end{aligned}\]

The nonlinear nonzero count (a measure of how nonlinear the model is) goes down:

MODEL STATISTICS

BLOCKS OF EQUATIONS           2     SINGLE EQUATIONS          135

BLOCKS OF VARIABLES           3     SINGLE VARIABLES          111

NON ZERO ELEMENTS           405     NON LINEAR N-Z             90

CODE LENGTH                 361     CONSTANT POOL              16

I hoped this would help somewhat. It really didn't. We see:

    Elapsed time (sec): 962, estimated tree completion: 0.01635

        Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time

     4902062  1774.764542  5413.400780     11 3749096     48   67.22%       9    998

More of the tree has been explored, but the gap actually deteriorated. Usually, these measures go in the same direction, but not always, as we see. We found the same best solution, as will be the case for all quadratic experiments. That is actually an indication (but not proof) that this solution is optimal or near-optimal.

Model 3: adding a symmetry breaker

Here, we add an ordering constraint to model 1. In model 1, we could just renumber (or rename) the points without changing the situation. Requiring that solutions are ordered by their \(x\)-coordinate can help. 

Model 3
\[\begin{aligned}\max\>&{\color{darkred}z}\\ & {\color{darkred}z} \le \sum_c \left({\color{darkred}x}_{i,c}-{\color{darkred}x}_{j,c}\right)^2 && \forall i\lt j \\ &{\color{darkred}x}_{i,x} \le{\color{darkred}x}_{i+1,x} && \forall i\lt n\\ & {\color{darkred}x}_{i,c} \in [{\color{darkblue}L},{\color{darkblue}U}] \\ & c \in \{x,y\} \\ & i,j \in \{1,\dots,n\} \end{aligned}\]

 Results:

    Elapsed time (sec): 984, estimated tree completion: 0.25586

        Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time

     3113456  1774.764542  2782.885513     10 1379914     34   36.23%      13    999

This is an improvement. And offers some hope. 

In the solution, points with smaller indices \(i\) should be further to the left. This ordering constraint is also used in [2].

Model 4: refining the symmetry breaker

In the previous model, we considered only the \(x\)-coordinate in the ordering constraint. I try to refine that a bit, so that the \(y\)-coordinate also plays a role.

Model 4

\[\begin{aligned}\max\>&{\color{darkred}z}\\ & {\color{darkred}z} \le \sum_c \left({\color{darkred}x}_{i,c}-{\color{darkred}x}_{j,c}\right)^2 && \forall i\lt j \\ &{\color{darkred}x}_{i,x}+0.123{\color{darkred}x}_{i,y} \le{\color{darkred}x}_{i+1,x}+0.123 {\color{darkred}x}_{i+1,y} && \forall i\lt n\\ & {\color{darkred}x}_{i,c} \in [{\color{darkblue}L},{\color{darkblue}U}] \\ & c \in \{x,y\} \\ & i,j \in \{1,\dots,n\} \end{aligned}\]

Why 0.123? No good reason. I suspect that the actual value is not very important. I am just trying to make things unique. Let's see:

    Elapsed time (sec): 988, estimated tree completion: 0.55634

        Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time

     2462808  1774.764542  2484.977151     10 886081     32   28.58%      12    999

That looks definitely better. We have explored more than half of the tree. And the gap has been reduced further. The best found solution is still the same. If we want to prove optimality, this would be a good candidate.

Model 5: Manhattan distance

If we want to use a different metric for the distance between neighbors, e.g., the distance based on the L1 norm, we get a model like:

Model 5
\[\begin{aligned}\max\>&{\color{darkred}z}\\ & {\color{darkred}z} \le \sum_c \left|{\color{darkred}x}_{i,c}-{\color{darkred}x}_{j,c}\right| && \forall i\lt j \\ & {\color{darkred}x}_{i,c} \in [{\color{darkblue}L},{\color{darkblue}U}] \\ & c \in \{x,y\} \\ & i,j \in \{1,\dots,n\} \end{aligned}\]

We can actually use the \(\mathbf{abs}(.)\) function directly in Xpress. That is very nice, as maximizing an L1 norm requires some attention. When I ran this as stated, I saw: 

    ------ unbounded -------

    Concurrent statistics:

               Dual: 90 simplex iterations, 0.00s

    Problem is unbounded if C212 enters the basis

     *** Search completed ***

    Problem is unbounded

 

I think Xpress should be able to handle this a bit better. If an LP relaxation is unbounded, the underlying MIP is not necessarily unbounded. Luckily, we can easily fix this by adding an upper on \({\color{darkred}z}\): \[{\color{darkred}z} \le 1\text{e}6\] Now, it actually reformulates the problem as a linear MIP model. I had high hopes for this, but unfortunately:

    Elapsed time (sec): 999, estimated tree completion: 0.35231    

        Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time

     26588858    57.142857    75.000000     17 9164865     31   23.81%      18    999

Obviously, the best solution is different. After all, we used a different objective. But still no proof of global optimality. Note that Gurobi does not accept this formulation; it requires something simpler, such as \(y=|x|\). This is really somewhat primitive and prevents me from writing things as I would like. I am not completely sure if I should blame Gurobi or GAMS for this.

Model 6: Manhattan distance + symmetry breaker

Here we use:

Model 6

\[\begin{aligned}\max\>&{\color{darkred}z}\\ & {\color{darkred}z} \le \sum_c \left|{\color{darkred}x}_{i,c}-{\color{darkred}x}_{j,c}\right| && \forall i\lt j\\ &{\color{darkred}x}_{i,x}+0.123{\color{darkred}x}_{i,y} \le{\color{darkred}x}_{i+1,x}+0.123 {\color{darkred}x}_{i+1,y} && \forall i\lt n\\ & {\color{darkred}x}_{i,c} \in [{\color{darkblue}L},{\color{darkblue}U}] \\ & c \in \{x,y\} \\ & i,j \in \{1,\dots,n\} \end{aligned}\]

This shows:

    Elapsed time (sec): 1, estimated tree completion: 0.98114

        Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time

       59659    57.142857    58.713838     18     42     27    2.68%       9      2

  *** Search completed ***

Wow, that is quite something. Proven optimality in 2 seconds. Here we see the impact of this symmetry breaker. From hitting a time limit of 1,000 seconds to optimality in 2 seconds. Note that the best solution is the same as in the previous model.

Model 7: Manhattan distance + binary variables

Let's try to implement this ourselves using binary variables. Will I be able to beat what Xpress is doing automatically?

Model 7

\[\begin{aligned}\max\>&{\color{darkred}z}\\ & {\color{darkred}z} \le \sum_c \left( {\color{darkred}d}_{i,j,c}^+ +{\color{darkred}d}_{i,j,c}^- \right)  && \forall i\lt j \\ &{\color{darkred}d}_{i,j,c}^+ -{\color{darkred}d}_{i,j,c}^- = {\color{darkred}x}_{i,c}-{\color{darkred}x}_{j,c} && \forall i \lt j \\& {\color{darkred}d}_{i,j,c}^+ \le {\color{darkblue}M} {\color{darkred}\delta}_{i,j,c} && \forall i \lt j, c \\ & {\color{darkred}d}_{i,j,c}^- \le {\color{darkblue}M} (1-{\color{darkred}\delta}_{i,j,c}) && \forall i \lt j, c \\ &{\color{darkred}x}_{i,x}+0.123{\color{darkred}x}_{i,y} \le{\color{darkred}x}_{i+1,x}+0.123 {\color{darkred}x}_{i+1,y} && \forall i\lt n\\ & {\color{darkred}x}_{i,c} \in [{\color{darkblue}L},{\color{darkblue}U}] \\ & {\color{darkred}d}_{i,j,c}^+,{\color{darkred}d}_{i,j,c}^- \ge 0 \\ & {\color{darkred}\delta}_{i,j,c} \in\{0,1\} \\& {\color{darkblue}M} = 2({\color{darkblue}U}-{\color{darkblue}L})  \\ & c \in \{x,y\} \\ & i,j \in \{1,\dots,n\} \end{aligned}\]

This is more work than we needed when we used \(\mathbf{abs}(.)\) directly. We see:

    Elapsed time (sec): 2, estimated tree completion: 0.89402

        Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time

       28397    57.142857    72.001382     16     57     20   20.64%      32      3

     *** Search completed ***

This needed a bit more time. The number of nodes required seems significantly lower than that of the previous model. So it is not clear who the winner is.

Conclusion

A simple problem where we need to place 10 points on a square is very difficult. We find good solutions (probably optimal or near optimal) very quickly. But proving global optimality is another kettle of fish.

Instead of adding a symmetry breaker only to the \(x\)-coordinate, I suggest also using the \(y\)-coordinate. That seems to improve progress significantly, both on the quadratic and the linear model. This ordering approach is new for me. I think I will use this more often. It really seems to make a difference compared to using just \(x\).

Model

GAMS Models

$ontext      Hostile Brothers Problem     Model 1, 10 points 1000 seconds:     Elapsed time (sec): 999, estimated tree completion: 0.00015         Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time      2483289  1774.764542  4812.031822     11 2420538     24   63.12%      25    999   Model 2, 10 points 1000 seconds:     Elapsed time (sec): 962, estimated tree completion: 0.01635         Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time      4902062  1774.764542  5413.400780     11 3749096     48   67.22%       9    998   Model 3, 10 points 1000 seconds:     Elapsed time (sec): 984, estimated tree completion: 0.25586         Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time      3113456  1774.764542  2782.885513     10 1379914     34   36.23%      13    999   Model 4, 10 points 1000 seconds:     Elapsed time (sec): 988, estimated tree completion: 0.55634         Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time      2462808  1774.764542  2484.977151     10 886081     32   28.58%      12    999   Model 5, 10 points 1000 seconds:     Elapsed time (sec): 999, estimated tree completion: 0.35231            Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time      26588858    57.142857    75.000000     17 9164865     31   23.81%      18    999   Model 6, 10 points 1000 seconds:     Elapsed time (sec): 1, estimated tree completion: 0.98114         Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time        59659    57.142857    58.713838     18     42     27    2.68%       9      2   *** Search completed ***   Model 7, 10 points 1000 seconds:     Elapsed time (sec): 2, estimated tree completion: 0.89402         Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time        28397    57.142857    72.001382     16     57     20   20.64%      32      3      *** Search completed ***     $offtext   * solver option qcp=xpress, dnlp=xpress, mip=xpress;   * time limit option reslim=1000;   sets     i /point1*point10/     c /x,y/     pm /'+','-'/ ; alias (i,j);   set lt(i,j) 'upper triangular part'; lt(i,j) = ord(i) < ord(j);   scalars     L /0/     U /100/ ;   variables     z            'objective value'     x(i,c)       'location of point i in dimension c'     d(i,j,c)     'distance between points i and j in dimension c'     d2(i,j,c,pm) 'split variables for linearizing manhattan distance' ;   x.lo(i,c) = L; x.up(i,c) = U;   positive variable d2;   binary variable delta(i,j,c);   equation    sqdistance   'squared distance'    edelta       'definition of d'    sqdistance2  'squared distance (alternative, using d)'    order        'ordering of points by x-coordinate'    order2       'ordering of points by x- and y-coordinate (alternative)'    l1distance   'manhattan distance'    l1distance2  'manhattan distance (linearized)'    split1       'variable splitting constraint 1'    split2       'variable splitting constraint 2'    split3       'variable splitting constraint 3' ;     sqdistance(lt(i,j)).. z =l= sum(c, sqr(x(i,c) - x(j,c)));   edelta(lt(i,j),c).. d(i,j,c) =e= x(i,c) - x(j,c); sqdistance2(lt(i,j)).. z =l= sum(c, sqr(d(i,j,c)));   order(i+1).. x(i,'x') =l= x(i+1,'x');   order2(i+1).. x(i,'x')+0.123*x(i,'y') =l= x(i+1,'x')+0.123*x(i+1,'y');   * manhattan distance l1distance(lt(i,j)).. z =l= sum(c, abs(x(i,c) - x(j,c)));   * needed for xpress solver z.up = 1e6;     scalar M 'big-M'; M = 2*(U-L); * manhattan distance manually linearized split1(lt(i,j),c).. d2(i,j,c,'+') - d2(i,j,c,'-') =e= x(i,c) - x(j,c); split2(lt(i,j),c).. d2(i,j,c,'+') =l= M*delta(i,j,c); split3(lt(i,j),c).. d2(i,j,c,'-') =l= M*(1-delta(i,j,c)); l1distance2(lt(i,j)).. z =l= sum((c,pm), d2(i,j,c,pm));     model brothers1 /sqdistance/; model brothers2 /edelta, sqdistance2/; model brothers3 /sqdistance, order/; model brothers4 /sqdistance, order2/; model brothers5 /l1distance/; model brothers6 /l1distance, order2/; model brothers7 /l1distance2, split1,split2,split3, order2/;   solve brothers1 maximizing z using qcp; solve brothers2 maximizing z using qcp; solve brothers3 maximizing z using qcp; solve brothers4 maximizing z using qcp; solve brothers5 maximizing z using dnlp; solve brothers6 maximizing z using dnlp; solve brothers7 maximizing z using mip;

 

References

1. Cornuejols, Gerard; Karamanov, Miroslav; Li, Yanjun (2006). Early Estimates of the Size of Branch-and-Bound Trees. Carnegie Mellon University. Journal contribution. https://doi.org/10.1184/R1/6705221.v1
2. [2605.04850] Out-of-the-Box Global Optimization for Packing Problems: New Models and Improved Solutions