Another very small but very difficult global NLP model

The goal of this exercise is to fill a square area \([0,250]\times[0,100]\) with 25 circles. The model can choose the \(x\) and \(y\) coordinates of the center of each circle and the radius. So we have as variables \(\color{darkred}x_i\), \(\color{darkred}y_i\), and \(\color{darkred}r_i\). The circles placed inside the area should not overlap. The objective is to maximize the total area covered. 

A solution is:

The optimization model can look like: 

Non-convex Quadratic Model

\[\begin{align} \max&\sum_i \color{darkblue}\pi\cdot\color{darkred}r_i^2 && &&\text{maximize covered area}\\ &(\color{darkred}x_i-\color{darkred}x_j)^2 + (\color{darkred}y_i-\color{darkred}y_j)^2 \ge (\color{darkred}r_i+\color{darkred}r_j)^2 && \forall i\lt j && \text{no overlap} \\ & \color{darkred}x_i-\color{darkred}r_i \ge 0 && && \text{stay inside area} \\ & \color{darkred}x_i+\color{darkred}r_i \le \color{darkblue}{\mathit{sizeX}} \\ & \color{darkred}y_i-\color{darkred}r_i \ge 0 \\ & \color{darkred}y_i+\color{darkred}r_i \le \color{darkblue}{\mathit{sizeY}} \\ & \color{darkred}r_i \le \color{darkred}r_{i-1} && \forall i\gt 1&& \text{optional symmetry breaker}\\ & \color{darkred}x_i \in [0,\color{darkblue}{\mathit{sizeX}}] && && \text{bounds} \\ &\color{darkred}y_i \in [0,\color{darkblue}{\mathit{sizeY}}] \\ & \color{darkred}r_i \in \left[0,\frac{\min(\color{darkblue}{\mathit{sizeX}},\color{darkblue}{\mathit{sizeY}})}{2}\right]\end{align}\]

The objective is simply maximizing the area covered by all circles. Actually, I reformulated this into an objective that measures the percentage of the covered area. The pairwise no-overlap constraints are defined over only \(i \lt j\) to prevent checking a pair of circles twice. The symmetry-breaking constraint reduces the size of the problem. Note that we have both a non-convex quadratic objective and a non-convex no-overlap constraint.

The model is small: 75 variables and 424 constraints. However, it is somewhat of a nightmare for global solvers. It is extremely difficult to prove optimality. We can find an excellent starting solution by using a simple multi-start algorithm with a local NLP solver. In fact, in this experiment, the multi-start algorithm found the best solution. The global solvers did not improve on this. 

The results are:

----    156 VARIABLE z.L                   =       91.183  covered area (% of area)

----    162 PARAMETER results  

                   x           y           r

circle1      200.000      50.000      50.000
circle2       50.000      50.000      50.000
circle3      125.000      28.125      28.125
circle4      116.667      77.778      22.222
circle5      150.000      87.500      12.500
circle6       92.857       9.184       9.184
circle7      157.143       9.184       9.184
circle8      142.308      60.577       8.654
circle9        8.579      91.421       8.579
circle10     241.421       8.579       8.579
circle11     241.421      91.421       8.579
circle12       8.579       8.579       8.579
circle13      90.000      92.000       8.000
circle14     166.667      94.444       5.556
circle15     104.545      53.766       4.675
circle16      80.000       4.500       4.500
circle17     170.000       4.500       4.500
circle18       4.289      20.711       4.289
circle19     245.711      20.711       4.289
circle20       4.289      79.289       4.289
circle21     150.000      70.833       4.167
circle22      78.571      95.918       4.082
circle23     142.308      73.077       3.846
circle24     155.405      21.284       3.041
circle25       2.566      72.654       2.566

Conclusions

• These models are very difficult for global solvers. 
• This is a nice benchmark problem as we can scale it from very small to very large. In addition, humans can easily assess the quality of the solution by visual inspection of a plot. We can easily spot bad solutions. Even very small instances can be very tough to solve to proven optimality.
• Multi-start NLP is quite effective in finding good solutions. 
• GAMS support for multi-start NLP algorithms is limited and can be improved. We want to generate the model just once and use parallel execution. This use case would be a good candidate for the scenario solver Guss. It is inexplicable why this is not supported. 
• The optional symmetry-breaking constraint is beneficial when using global solvers. Global solvers behave more like MIP solvers, where such constraints can be highly effective. A tiny 4 circle problem yielded the following timings:
  
   Without symmetry constraint:
               313360 Nodes explored
              1968.00 Total time (wall s)
  
  With symmetry constraint:
                25759 Nodes explored
                86.00 Total time (wall s)
     
• However, the symmetry constraint hurts when using the multi-start algorithm. Local solvers take more time to become feasible and, in general, find worse local optima. For the local search, it is better to turn off the symmetry constraint. To illustrate, here are the results of 100 trials with 10 circles:
  
  With symmetry constraint:
  ----    112 PARAMETER zbest                =       83.686  best objective
  
  Without symmetry constraint:
  ----    112 PARAMETER zbest                =       86.096  best objective
  
  
• A solver like Baron has a built-in facility to perform a multi-start local search before it really starts working on the global problem. Here, I argue that it may be better to separate these phases.
• This is an interesting model. It is simple to explain, but there are a few angles that are not immediately obvious.

Appendix: GAMS model

$onText   Place 25 circles (size is endogenous) on a [0,250]x[0,100] square area such that circles don't overlap and the total area covered is maximized.   Very tough little non-convex NLP.               Obj MultiStart  91.183   (1015 trials, w/o symmetry constraint) Baron       91.183   (no improvement, 1 hour, with symmetry breaker)   Reference: https://yetanothermathprogrammingconsultant.blogspot.com/2024/05/another-very-small-but-very-difficult.html   $offText   *------------------------------------------------------------------- * size of problem *-------------------------------------------------------------------   Set i 'circles' /circle1*circle25/; alias (i,j);   set ij(i,j) 'compare circles i and j: only i<j'; ij(i,j) = ord(i) < ord(j);   scalars   sizeX 'size of area' /250/   sizeY 'size of area' /100/ ;   *------------------------------------------------------------------- * NLP Model *-------------------------------------------------------------------   positive variables    x(i) 'x-coordinate'    y(i) 'y-coordinate'    r(i) 'radius' ;   x.up(i) = sizeX; y.up(i) = sizeY; r.up(i) = min(sizeX,sizeY)/2;   variable z 'covered area (% of area)';   equations     cover  'calculate total area covered'     no_overlap(i,j) 'circles cannot overlap'     xlo(i) 'stay inside [0,250]'     xup(i) 'stay inside [0,250]'     ylo(i) 'stay inside [0,100]'     yup(i) 'stay inside [0,100]'     symmetry(i) 'symmetry breaker' ;     cover.. z =e= 100*sum(i, pi*sqr(r(i))) / (sizeX*sizeY); no_overlap(ij(i,j)).. sqr(x(i)-x(j)) + sqr(y(i)-y(j)) =g= sqr(r(i)+r(j)); xlo(i).. x(i) - r(i) =g= 0; xup(i).. x(i) + r(i) =l= sizeX; ylo(i).. y(i) - r(i) =g= 0; yup(i).. y(i) + r(i) =l= sizeY; symmetry(i-1) .. r(i) =l= r(i-1);     * improves logging of gap, but may slow things down * (nl obj becomes a constraint) z.up = 100;     model circlesw  'with symmetry' /all/ ; model circleswo 'w/o symmetry' /circlesw - symmetry/ ;     *------------------------------------------------------------------- * select algorithm *-------------------------------------------------------------------   $set MultiStart 1 $set Global     0     *------------------------------------------------------------------- * multistart approach using local NLP solver * best used without symmetry breaker * we get a solution:  91.183 @ trial1015 *-------------------------------------------------------------------   $ifThen %MultiStart% == 1   set k /trial1*trial1015/; option qcp = conopt; parameter best(i,*) 'best solution'; scalar zbest 'best objective' /0/; parameter trace(k) 'keep track of improvements'; circleswo.solprint = %solprint.Silent%; circleswo.solvelink = 5; * abort expensive solves option reslim = 10; loop(k,      x.l(i) = uniform(10,sizeX-10);      y.l(i) = uniform(10,sizeY-10);      r.l(i) = uniform(10,40);           solve circleswo maximizing z using qcp;      if(circleswo.modelstat=%modelStat.locallyOptimal% and z.l>zbest,          zbest = z.l;          best(i,'x') = x.l(i);          best(i,'y') = y.l(i);          best(i,'r') = r.l(i);          trace(k) = z.l;     );  );   display zbest,best,trace;   * ugly code to sort on r so we can use symmetry constraint in * subsequent solve. set rem(i) 'remaining'; rem(i) = yes; option strictSingleton = 0; singleton set cur(i) 'current'; scalar largest; loop(i,    largest = smax(rem,best(rem,'r'));    cur(rem) = best(rem,'r')=largest;    x.l(i) = best(cur,'x');    y.l(i) = best(cur,'y');    r.l(i) = best(cur,'r');    rem(cur) = no; ); z.l = zbest;   * reset to default for next solves circleswo.solprint = %solprint.On%;     $endIf     *------------------------------------------------------------------- * Global NLP solver * with or without symmetry breaker *-------------------------------------------------------------------   $ifThen %Global% == 1   option qcp=baron, reslim=3600; option threads=-1; solve circlesw maximizing z using qcp;   $endIf   *------------------------------------------------------------------- * Reporting *-------------------------------------------------------------------   display z.l;   parameter results(i,*); results(i,'x') = x.l(i); results(i,'y') = y.l(i); results(i,'r') = r.l(i); display results;

 

Notes: 

• When we find a solution with the multi-start local solver algorithm, we have the \(\color{darkred}r_i\) unsorted. To be able to use this solution as a starting point for the global solver we (assuming it uses the symmetry breaker constraint), we need to sort the solution so it is feasible w.r.t. to the symmetry constraint. 

Appendix 2: related example

In [1], a similar but slightly easier problem is posed. Given \(n\) circles, some of them overlapping, try to move them around a bit such that they are no longer overlapping. The objective is to minimize the total distance the circles are moved. The model can look like:

Non-convex Quadratic Model 2
\[\begin{align} \min&\sum_i  \left[ \left(\color{darkred}x_i-\color{darkblue}x_i^0\right)^2 +\left(\color{darkred}y_i-\color{darkblue}y_i^0\right)^2 \right] && &&\text{minimize squared distances}\\ &(\color{darkred}x_i-\color{darkred}x_j)^2 + (\color{darkred}y_i-\color{darkred}y_j)^2 \ge (\color{darkblue}r_i+\color{darkblue}r_j)^2 && \forall i\lt j && \text{no overlap} \end{align}\]

The objective is slightly changed: we dropped the \(\sqrt{.}\) to make it quadratic. (The square root function is a bit difficult: it is not defined for negative arguments, and the derivative does not exist at zero).

Let's generate a small data set:

----     53 PARAMETER d  reporting of data+results

                  x0          y0           r

circle1      -65.651      99.624      24.388
circle2       68.653      15.747      24.058
circle3       10.075      98.227      15.260
circle4      -39.772      52.450      16.004
circle5      -41.558     -73.862      33.565
circle6      -55.189      27.944      43.236
circle7      -30.034     -68.096      19.233
circle8       71.254     -49.984      36.629
circle9      -86.577      33.786      41.034
circle10       0.042     -12.929      22.146

Data: initial locations

After solving the problem, we should see:

Solution

This data set is small enough to prove global optimality. For larger problems, it may be better to use a multi-start NLP model or stop the global solver early.

References

1. Move circles from initial positions such that they do not overlap and total shift is minimized, https://or.stackexchange.com/questions/12052/move-circles-from-initial-positions-such-that-they-do-not-overlap-and-total-shif

GAMS Model

$onText      Problem  We have a number of overlapping circles. We want to move around  a bit (as little as possible) to make them non-overlapping.    This is a nonconvex problem.   $offText     *------------------------------------------------------------------- * size of problem *-------------------------------------------------------------------   set i /circle1*circle10/; alias(i,j);   set ij(i,j) 'i<j'; ij(i,j) = ord(i) < ord(j);   table bounds(*,*) 'limits for (x0,y0) and (x,y)'         min  max   x0   -100  100   y0   -100  100   r     10   50   x    -200  200   y    -200  200 ;   * (x0,y0) and r are data * (x,y) are variables      *------------------------------------------------------------------- * data *-------------------------------------------------------------------   Parameter    x0(i)  'prior x-coordinate'    y0(i)  'prior y-coordinate'    r(i)   'radius'    d(i,*) 'reporting of data+results' ; x0(i) = uniform(bounds('x0','min'),bounds('x0','max')); y0(i) = uniform(bounds('y0','min'),bounds('y0','max')); r(i) = uniform(bounds('r','min'),bounds('r','max'));   d(i,'x0') = x0(i); d(i,'y0') = y0(i); d(i,'r') = r(i); display d;   *------------------------------------------------------------------- * model *-------------------------------------------------------------------   variables   x(i) 'new x-coordinates'   y(i) 'new y-coordinates'   z 'objective: sum of squared distances' ;   * for global NLPs it is always advised to * to use proper bounds on all variables x.lo(i) = bounds('x','min'); x.up(i) = bounds('x','max'); y.lo(i) = bounds('y','min'); y.up(i) = bounds('y','max');     Equations    obj     'objective: (x,y) as close as possible to (x0,y0)'    no_overlap(i,j) 'circles i and j should not overlap (non-convex constraint)' ;   obj..    z =e= sum(i, sqr(x(i)-x0(i)) + sqr(y(i)-y0(i)));   no_overlap(ij(i,j))..    sqr(x(i)-x(j)) + sqr(y(i)-y(j)) =g= sqr(r(i)+r(j));   model move /all/;   *------------------------------------------------------------------- * select algorithm *-------------------------------------------------------------------   $set MultiStart 1 $set Global     0     *------------------------------------------------------------------- * multistart approach using local NLP solver *-------------------------------------------------------------------   $ifThen %MultiStart% == 1   set k /trial1*trial25/; option qcp = conopt; parameter best(i,*) 'best solution'; scalar zbest 'best objective' /INF/; parameter trace(k) 'keep track of improvements'; move.solprint = %solprint.Silent%; move.solvelink = 5; loop(k,      x.l(i) = uniform(bounds('x','min')+10,bounds('x','max')-10);      y.l(i) = uniform(bounds('y','min')+10,bounds('y','max')-10);           solve move minimizing z using qcp;      if(move.modelstat=%modelStat.locallyOptimal% and z.l<zbest,          zbest = z.l;          best(i,'x') = x.l(i);          best(i,'y') = y.l(i);          trace(k) = z.l;     );  );   display zbest,best,trace;   x.l(i) = best(i,'x'); y.l(i) = best(i,'y');     * reset to default for next solves move.solprint = %solprint.On%;   $endIf     *------------------------------------------------------------------- * Global NLP solver * with or without symmetry breaker *-------------------------------------------------------------------   $ifThen %Global% == 1   option qcp=scip, reslim=1000; option threads=-1; solve move minimizing z using qcp;   $endIf    *------------------------------------------------------------------- * reporting *-------------------------------------------------------------------   d(i,'x') = x.l(i); d(i,'y') = y.l(i); d(i,'sq.dist') = sqr(x.l(i)-x0(i)) + sqr(y.l(i)-y0(i)); display d;