## Wednesday, May 15, 2024

### 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:

\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 thatcircles don't overlap and the total area covered is maximized. Very tough little non-convex NLP. ObjMultiStart 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 izbest,         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 solvescircleswo.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

1. Could we not apply combinatorial benders cuts as was done in the similar post here ? https://yetanothermathprogrammingconsultant.blogspot.com/2018/05/knapsack-packing-combinatorial-benders.html

1. I don't know how that would work as there are only continuous variables here.

2. What real life situation/example would require a model like this?

3. using GRG, i got a better solution:
x y r
125 71.875 28.125
4.289321881 79.28932188 4.289321881
133.3333333 22.22222222 22.22222222
4.289321881 20.71067812 4.289321881
150 3.125 3.125
145.4545455 46.23376623 4.675324675
92.85714286 90.81632653 9.183673469
245.7106781 79.28932188 4.289321881
102.9411765 50 2.941176471
114.2857143 4.081632653 4.081632653
107.6923077 39.42307692 8.653846154
100 29.16666667 4.166666667
200 50 50
83.33333333 5.555555556 5.555555556
8.578643763 91.42135624 8.578643763
170 95.5 4.5
100 12.5 12.5
157.1428571 90.81632653 9.183673469
245.7106781 20.71067812 4.289321881
107.6923077 26.92307692 3.846153846
50 50 50
160 8 8
241.4213562 8.578643763 8.578643763
8.578643763 8.578643763 8.578643763
241.4213562 91.42135624 8.578643763

1. Slightly better (0.01%). Multi-start is a heuristic, so that is to be expected. And, GRG2 is a good NLP solver.