Circle Packing and HTML reporting

Little example. Here, we try to pack \(n\) circles with a given radius \(r_i\) into a larger disc with an unknown radius \(R\). The goal is to minimize \(R\). The underlying model is simple: 

Packing of Circles
\[\begin{align} \min\> & \color{darkred}R \\ & \sum_c \left(\color{darkred}p_{i,c}-\color{darkred}p_{j,c}\right)^2 \ge \left(\color{darkblue}r_i+\color{darkblue}r_j\right)^2 & \forall i\lt j \\ & \sum_c \color{darkred}p_{i,c}^2 \le \left(\color{darkred}R-\color{darkblue}r_i\right)^2 & \forall i \\ & \color{darkred}R \ge 0\\ & c \in \{x,y\} \\ \end{align}\]

The first constraint says that inner circles cannot overlap. The second constraint states that inner circles should be inside the outer disc. These types of models are very difficult for global solvers. However, very good solutions can be found using a simple multistart NLP approach: solve several trials with a different random starting point.

Data

The data we need is a collection of radii of the circles we want to pack. The random data set was created by drawing radii from a log-normal distribution. This generates relatively few large circles and many small ones. Drawing from a uniform distribution may generate more challenging data sets (relatively fewer smaller circles). Our data set has 50 circles.

----     30 PARAMETER radius  lognormally distributed

circle1  2.324,    circle2  3.202,    circle3  3.427,    circle4  1.089,    circle5  1.886,    circle6  1.672,    circle7  2.232
circle8  4.337,    circle9  1.860,    circle10 2.719,    circle11 1.726,    circle12 3.228,    circle13 1.952,    circle14 2.132
circle15 6.971,    circle16 1.896,    circle17 1.420,    circle18 1.367,    circle19 3.189,    circle20 4.921,    circle21 1.527
circle22 3.877,    circle23 4.521,    circle24 1.893,    circle25 1.649,    circle26 2.727,    circle27 1.653,    circle28 1.727
circle29 3.786,    circle30 5.963,    circle31 8.221,    circle32 1.989,    circle33 2.437,    circle34 1.680,    circle35 2.576
circle36 2.384,    circle37 2.089,    circle38 1.619,    circle39 7.323,    circle40 4.109,    circle41 4.492,    circle42 1.866
circle43 3.967,    circle44 4.954,    circle45 4.144,    circle46 1.329,    circle47 7.765,    circle48 3.806,    circle49 2.958
circle50 6.200

Solve loops

A serial loop is the simplest. We just solve \(n\) optimization problems with a random starting point for the locations (x,y coordinates) of the circles. A noteworthy detail: it is important to set the initial value of \(R\) (the outer radius) to a nonzero value (the default is zero!). Without this, CONOPT (the local NLP solver I am using) is having quite some difficulties in finding a feasible solution.

An obvious improvement is to solve the subproblems in a parallel fashion. GAMS has a so-called scenario solver. Unfortunately, it does not have any capabilities to provide different initial values, so it is not very useful for this case. 

The GAMS tools for parallel solve loops are low-level. IMHO, the developers chose the wrong abstraction level, requiring us to use multiple loops. One extra facility I added is a limit on the number of jobs that can run simultaneously. In our case, this was set to the number of cores, but, of course, one can use any number. This can be useful in cases where we have large jobs and machines with relatively tight memory. A not uncommon situation, especially for laptops. The GAMS model below can be used as an illustration of how to set up such a parallel asynchronous solve loop in GAMS. 

It is noted that some solvers have built-in multistart capabilities. Baron and Knitro are examples.

Visualization

Of course, we need some visualization to assess the results. Here, I use simple HTML and SVG (Scalable Vector Graphics), which are easy to generate. An advantage: rendering HTML does not require installing special software. In some corporate environments, this can be a big plus. In addition, HTML is easy to generate (just write a text file) and debug (the text file is human-readable). It is quite low-tech, but that is an advantage, I think.

The results of my experiments look like:

Even looking at the single optimization job, CONOPT, the NLP solver I used, was doing a great job of packing the circles tightly. 

Conclusions

Here, we try to pack \(n=50\) circles of different sizes into one big outer circle. We minimize the radius of the outer circle such that all \(n\) inner circles don't overlap and are inside the outer circle. This is a difficult model for a global NLP solver. We use a simple multi-start heuristic using the Conopt local NLP solver. We compare and visualize the results of three experiments.

1. A single solve
2. A serial multi-start loop with 25 solves
3. A parallel multi-start loop with 250 solves

Conopt is doing a remarkably good job on this. The results are visualized using simple HTML+SVG.

Appendix: GAMS Model

$onText      Pack n circles of given radius    within a disc. Minimize the    size of this disc.      This problem has nonconvex constraints.       We have three approaches:       1. single solve       2. serial loop: 25 solves       3. parallel loop, 250 solves. No more than NumCores jobs at the time.          $offtext   $set numSerialJobs 25 $set numParallelJobs 250   *--------------------------------------- * data *---------------------------------------   set   i 'circles' /circle1*circle50/ ;   parameter radius(i) 'lognormally distributed'; radius(i) = exp(1+normal(0,0.5)); display radius;   *--------------------------------------- * additional sets *---------------------------------------   alias (i,j); sets    ij(i,j) 'compare i<j'    c 'coordinates' /x,y/    ; ij(i,j) = ord(i)<ord(j); display ij;   *--------------------------------------- * model *---------------------------------------   variables    r 'radius of outer circle'    p(i,c) 'position' ; r.lo = 0;   Equations   overlap(i,j) 'no-overlap'   inside(i)    'inside outer circle' ;   overlap(ij(i,j))..     sum(c, sqr(p(i,c)-p(j,c))) =g= sqr(radius(i)+radius(j));   inside(i)..     sum(c, sqr(p(i,c))) =l= sqr(r-radius(i));   model pack /all/;   option nlp=conopt; scalar t; t = timeElapsed; p.l(i,c) = uniform(-20,20); r.l = 100; solve pack minimizing r using nlp; display p.l,r.l;   abort$(pack.modelstat <> %modelStat.locallyOptimal%) "no solution";   parameter solutions(*,*,*); solutions('Single',i,c) = p.l(i,c); solutions('Single',i,'r') = radius(i); solutions('Single','outer','r') = r.l; solutions('Single','numjobs','statistics') = 1; solutions('Single','time','statistics') = TimeElapsed-t; display solutions;   *--------------------------------------- * multistart approach: serial version *---------------------------------------   set trials /trial1*trial%numSerialJobs%/;   scalar bestr /inf/;   * turn off info written to listing file pack.solprint = 2;     t = TimeElapsed; loop(trials,   p.l(i,c) = uniform(-10,10);   r.l = 100;   solve pack minimizing r using nlp;   continue$(pack.modelstat <> %modelStat.locallyOptimal%);   if (r.l<bestr,      bestr = r.l;      solutions(trials,'outer','r') = r.l;      solutions('Multistart Serial',i,c) = p.l(i,c);      solutions('Multistart Serial',i,'r') = radius(i);      solutions('Multistart Serial','outer','r') = r.l;   );   );   solutions('Multistart Serial','numjobs','statistics') = %numSerialJobs%; solutions('Multistart Serial','time','statistics') = TimeElapsed-t; display solutions;   *--------------------------------------- * multistart approach: parallel version *---------------------------------------   set    ptrials /ptrial1*ptrial%numParallelJobs%/    launched(ptrials) 'jobj that have started'    TimeIntervals '0.1 second intervals' /t1*t10000/ ; launched(ptrials) = no;   Parameter    h(ptrials)   'model handles'    sleepresult  'capture rerult of sleep function (not used)'    numRunning   'number of running jobs'    maxRunning   'max number of jobs running at the same time'    MaxWaitTime  'abandon collection loop after this much time' /7200/ ; maxRunning = numCores; h(ptrials)=0;   t = timeElapsed;   * * job submission loop * pack.solveLink = 3; loop(ptrials,    p.l(i,c) = uniform(-10,10);    r.l = 100;    loop(timeIntervals,       numrunning = sum(launched$(handleStatus(h(launched)) = 1), 1);       break$(numrunning < maxRunning);       sleepResult = sleep(0.1);         );    solve pack minimizing r using nlp;    h(ptrials) = pack.handle;    launched(ptrials) = yes; );   * * result collection loop * scalar StartCollectionTime; StartCollectionTime = TimeElapsed; bestr = inf; repeat    loop(ptrials$h(ptrials),       if(handleStatus(h(ptrials)) = 2,          pack.handle = h(ptrials);          execute_loadhandle pack;          if(pack.modelstat = %modelStat.locallyOptimal%,             if (r.l<bestr,                bestr = r.l;                solutions(ptrials,'outer','r') = r.l;                solutions('Multistart Parallel',i,c) = p.l(i,c);                solutions('Multistart Parallel',i,'r') = radius(i);                solutions('Multistart Parallel','outer','r') = r.l;             );            );                     display$handleDelete(h(ptrials)) 'trouble deleting handles';          h(ptrials) = 0;       );    );    break$(card(h)=0);    display$readyCollect(h, 1000) 'Problem waiting for next instance to complete'; until timeElapsed-StartCollectionTime > MaxWaitTime;   solutions('Multistart Parallel','numjobs','statistics') = %numParallelJobs%; solutions('Multistart Parallel','time','statistics') = TimeElapsed-t; display solutions;   *--------------------------------------- * global solver can do this in one solve * this is a difficult model for global * solvers. *--------------------------------------- $ontext option nlp=baron; solve pack minimizing r using nlp; $offtext   *--------------------------------------- * plot results *---------------------------------------   $set htmlfile 'result.html'   scalar r0 'radius of single job'; r0 = solutions('single','outer','r');   file f /%htmlfile%/; put f;   set run / Single, 'Multistart Serial','MultiStart Parallel'/;   put '<h1>Results of Circle Packing Model</h1>'/; put "<p>Pack ",card(i):0:0," circles with given radii in a container."     " We minimize the size of the outer black disc. </p>";     put "<table><tr>"/;   loop(run,    put '<td><svg width="400" height="300" viewBox="',         (-1.1*r0):0:3,' ',(-1.1*r0):0:3,' ',(2.2*r0):0:3,' ',(2.2*r0):0:3,'">'/;    put '<circle cx="0" cy="0" r="',solutions(run,'outer','r'):0:3,'"/>'/;       loop(i,         put '<circle cx="',solutions(run,i,'x'):0:3,                   '" cy="',solutions(run,i,'y'):0:3,                   '" r="',solutions(run,i,'r'):0:3,                   '" fill="darkorange"/>'/;    );       put '</svg><p style="text-align: center;"><b>',run.tl:0," Solve</b><br>",        " Number of optimization jobs: ",solutions(run,'numjobs','statistics'):0:0,"<br>",        " Time (seconds): ",solutions(run,'time','statistics'):0:1,"<br>",        " Length outer radius: ",solutions(run,'outer','r'):0:3,"<br>",        "</p></td>"/; );   put "</tr></table>"/; putclose;   executeTool 'win32.ShellExecute "%htmlfile%"'