Many mathematical modelling problems I face can be boiled down to constrained, multi-dimensional, global optimisation problems. Typically I have some model, with several parameters and some noisy data. I describe the relationship between the goodness-of-fit of a model to some data and the parameter values in the model as a mathematical function (an objective function) which I can then maximise within feasible parameter ranges. See my last post for an example of fitting a sine wave model to some made-up data.

The more parameters there are in the model, the harder the optimisation problem (particularly if there are few data!). Computational optimisation can take up a lot of CPU time, but maybe that's ok? CPUs get faster all the time, so maybe we can just wait for our next computer and the problem will go away naturally? Actually, over the past few years it has become apparent that individual CPUs are not really getting any faster, but instead we are seeing increased performance from computers and other devices through increased number of cores (effective CPUs) per device rather than increased speed per CPU. This shift in computational architecture represents a significant general programming challenge: how can we split up our problems to be solved by several independent CPUs and merge the resulting pieces efficiently?

This general problem is relevant for optimisation too. I recently built a new desktop computer. It has 6 cores, and with hyperthreading is capable of emulating 12 cores. How can I use all of these CPUs effectively to solve problems? To explore this, over the last few days, I've been looking at some optimisation libraries (so far in Python and R) in the hope of finding a satisfactory answer to this question. I will continue investigating and outline the results at a later date, but for now I want to present the optimisation problem that I am using for testing.

## Packing circles into a rectangle

Consider N circles fully enclosed by a rectangle of width W and height H. The problem is to maximise the area enclosed by the circles while ensuring they do not overlap each other, or cross the box boundaries. Let each circle $i$ have centre coordinates $(x_i,y_i)$ and radius $r_i$. We can write down the objective function (the fractional area enclosed by the circles, or packing density), which depends on $x_i$,$y_i$ and $r_i$, as:

$f(\overline{x})=\frac{ \pi\sum_{i=1}^{N} r_{i}^{2} }{WH}$subject to $(x_i,r_i)$ being at least $r_i$ away from the rectangle edges ($4N$ linear inequality constraints):

$\begin{align} g_{1,i}(\overline{x})&=x_{i}+r_{i}-W\geq 0\\ g_{2,i}(\overline{x})&=y_{i}+r_{i}-H\geq 0\\ g_{3,i}(\overline{x})&=r_{i}-x_{i}\geq 0\\ g_{4,i}(\overline{x})&=r_{i}-y_{i}\geq 0\\ \end{align} \forall i=1,...,N $and subject to no overlap between any circles ($N(1-N/2)$ non-linear inequality constraints):

$g_{5,i,j}(\overline{x})=r_i+r_j-\sqrt{(x_i-x_j)^2+(y_i-y_j)^2 }\geq 0 \hspace{1em} \forall i,j \hspace{1em} i=1,...,N \hspace{1em} i\neq j$The idea here is to maximise $f$, by choosing $x_i$, $y_i$ and $r_i$ while obeying the constraints $g_1,...,g_5$. We can also consider a simpler variant where instead of choosing $x_i$, $y_i$ and $r_i$, we only choose $x_i$ and $y_i$, and "inflate" each circle with a deterministic algorithm to satisfy constraint $g_5$. This variant has only linear constraints and only $2N$ parameters (instead of $3N$), but has an algorithmic objective function $f$.

These examples are useful, since a global optimum solution (where $N$ is also allowed to vary) is known, giving us some idea of achievable packing densities. Also, it is trivial to increase the number of parameters to make the problem arbitrarily more difficult (simply increase the number of circles $N$).

For now, here are some examples of some (sub-optimally) packed circles: