Application Center - Maplesoft
  Submit your work
Application Center Applications Circle Packing in an Ellipse Preview

App Preview:

Circle Packing in an Ellipse

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application




Circle Packing in an Ellipse

Introduction

This application finds the optimum packing of various size circles in an ellipse, such that the area of the ellipse is minimized.

 

 

This is a difficult global optimization problem and demands strong solvers. This application uses Maple's Global Optimization Toolbox

 

Circle packing (and packing optimization in general) is characterized by a large optimization space and many constraints; for this application, 35 circles generates 666 constraint equations.

 

The number of circles can be increased to create an increasingly complex problem; Maple automatically generates the symbolic constraint equations.

 

Applications like this are used to stress-test global optimizers.

 

The constraints and ellipse parameterization are taken from "Packing circles within ellipses", Birgin et al., International Transactions in Operational Research , Volume 20, Issue 3, pages 365–389, May 2013.

Setup

> 

restart
with(GlobalOptimization)

Number of circles

> 

n := 35

Radius of circles

> 

for i to n do r[i] := 1+.1*i end do

Area of ellipse (assuming the ellipse is described by the equation x^2/a^2+y^2/b^2 = 1

> 

area := Pi*a*b

Decision Variables and Optimization Bounds

See the reference in the introduction for a description of the decision variables. The position of the circles can be recovered from si, ui, vi, a and b (where a and b are the semi-axis of circumscribing ellipse)

> 

bounds := seq(s[i] = 0 .. 1, i = 1 .. n), seq(u[i] = -50 .. 50, i = 1 .. n), seq(v[i] = -50 .. 50, i = 1 .. n), a = 0 .. 50, b = 0 .. 50

Constraints

See the reference in the introduction for a description of the constraints.

> 

cons1 := seq((u[i]/a)^2+(v[i]/b)^2 = 1, i = 1 .. n)

> 

cons2 := seq(((s[i]-1)/r[i])^2*((b^2/a^2)^2*u[i]^2+v[i]^2) >= 1, i = 1 .. n)

> 

cons3 := seq(seq((((1+(s[i]-1)*b^2/a^2)*u[i]-(1+(s[j]-1)*b^2/a^2)*u[j])^2+(s[i]*v[i]-s[j]*v[j])^2)/(r[i]+r[j])^2 >= 1, j = i+1 .. n), i = 1 .. n-1)

Hence the entire set of constraints

> 

cons := {cons1, cons2, cons3, a >= 1.5*b}

Optimization and Results

> 

soln := GlobalSolve(area, cons, bounds, timelimit = 20, populationsize = 1000)

Hence the area of the optimized circumscribing ellipse is

> 

soln[1]

1281.27644396036112

(5.1)

Recover the coordinates of the packed circles

> 

X := eval([seq(x[i] = (1+(s[i]-1)*b^2/a^2)*u[i], i = 1 .. n)], soln[2])

> 

Y := eval([seq(y[i] = s[i]*v[i], i = 1 .. n)], soln[2])


Plot the circles and the optimized ellipse

> 

colorSpread := ColorTools:-EvenSpread("SteelBlue", n)

> 

circs := seq(plottools:-disk([rhs(select(has, X, x[i])[]), rhs(select(has, Y, y[i])[])], r[i], color = colorSpread[i], thickness = 2), i = 1 .. n)
> 

boundingEllipse := plottools:-ellipse([0, 0], rhs(select(has, soln[2], a)[]), rhs(select(has, soln[2], b)[]), color = black, thickness = 5, color = "LightGray", scaling = constrained)
> 

plots:-display(circs, boundingEllipse, scaling = constrained, axesfont = [Calibri])