Application Center - Maplesoft

App Preview:

Packing Circles into a Triangle

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

Learn about Maple
Download Application




Circle Packing in a Triangle

Introduction

This application finds the best packing and largest radius of equal-sized circles, such that they fit in a pre-defined triangle.

 

One solution, as visualized by this application, is illustrated below.

 

 

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, 20 circles generates 310 constraint equations.

 

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

 

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

Setup

restart
with(GlobalOptimization)

Number of circles

n := 15

Decision Variables and Optimization Bounds

The decision variables are the coordinates (x__i, y__i) of the centers of the circles, and their radius..

vars := [seq(x[i], i = 1 .. n), seq(y[i], i = 1 .. n), r]

Search bounds for the decision variables

bounds := seq(vars[i] = -3 .. 3, i = 1 .. 2*n), r = 0 .. 3

Vertices of the triangle

pts := [[1.3, 2.3], [-.4, .9], [.8, .8]]

Constraints

No two equal-radius circles overlap if the distance between their origins is twice their radius

cons1 := seq(seq(2*r <= sqrt((x[i]-x[j])^2+(y[i]-y[j])^2), j = i+1 .. n), i = 1 .. n-1)

Distance squared from the origin of a circle to any of the lines of the triangle must be greater than the radius squared

dist := proc (pts, x, y) options operator, arrow; {r^2 <= (x*(-pts[1, 2]+pts[2, 2])+pts[2, 1]*(pts[1, 2]-y)+pts[1, 1]*(-pts[2, 2]+y))^2/((pts[1, 1]-pts[2, 1])^2+(pts[1, 2]-pts[2, 2])^2), r^2 <= (x*(-pts[1, 2]+pts[3, 2])+pts[3, 1]*(pts[1, 2]-y)+pts[1, 1]*(-pts[3, 2]+y))^2/((pts[1, 1]-pts[3, 1])^2+(pts[1, 2]-pts[3, 2])^2), r^2 <= (x*(-pts[2, 2]+pts[3, 2])+pts[3, 1]*(pts[2, 2]-y)+pts[2, 1]*(-pts[3, 2]+y))^2/((pts[2, 1]-pts[3, 1])^2+(pts[2, 2]-pts[3, 2])^2)} end proc

cons2 := seq(dist(pts, x[i], y[i])[], i = 1 .. n)

The centre of each circle must be in the triangle.

circleCenters := proc (pts, x, y) options operator, arrow; {(pts[3, 1]*(-y+pts[2, 2])+pts[2, 1]*(y-pts[3, 2])+x*(-pts[2, 2]+pts[3, 2]))*signum(pts[3, 1]*(pts[1, 2]-pts[2, 2])+pts[1, 1]*(pts[2, 2]-pts[3, 2])+pts[2, 1]*(-pts[1, 2]+pts[3, 2])) <= 0, (pts[3, 1]*(y-pts[1, 2])+x*(pts[1, 2]-pts[3, 2])+pts[1, 1]*(-y+pts[3, 2]))*signum(pts[3, 1]*(pts[1, 2]-pts[2, 2])+pts[1, 1]*(pts[2, 2]-pts[3, 2])+pts[2, 1]*(-pts[1, 2]+pts[3, 2])) <= 0, -(pts[2, 1]*(y-pts[1, 2])+x*(pts[1, 2]-pts[2, 2])+pts[1, 1]*(-y+pts[2, 2]))*signum(-pts[3, 1]*(-pts[1, 2]+pts[2, 2])-pts[2, 1]*(pts[1, 2]-pts[3, 2])-pts[1, 1]*(-pts[2, 2]+pts[3, 2])) <= 0} end proc

cons3 := seq(circleCenters(pts, x[i], y[i])[], i = 1 .. n)

Collect all constraints

allcons := {cons1, cons2, cons3}

Optimization and Results

soln := GlobalSolve(r, allcons, bounds, maximize, timelimit = 25)

[.114143907837003108, [r = HFloat(0.11414390783700311), x[1] = HFloat(0.9308910993660217), x[2] = HFloat(-0.11064229721940184), x[3] = HFloat(0.6304712810716314), x[4] = HFloat(0.0655799701223353), x[5] = HFloat(0.21535102492190517), x[6] = HFloat(0.5942495451549262), x[7] = HFloat(0.7865092074471123), x[8] = HFloat(0.6676670927268661), x[9] = HFloat(0.858700153406567), x[10] = HFloat(0.4542490137298946), x[11] = HFloat(0.24180223746407262), x[12] = HFloat(0.7704690394892866), x[13] = HFloat(0.9528059147764201), x[14] = HFloat(0.4180245048058099), x[15] = HFloat(0.44356463126393225), y[1] = HFloat(1.553628027895338), y[2] = HFloat(0.9904264132537809), y[3] = HFloat(1.3422425417473725), y[4] = HFloat(1.1355506334175647), y[5] = HFloat(0.9632603030753388), y[6] = HFloat(1.570919926685674), y[7] = HFloat(1.12048235213861), y[8] = HFloat(0.9255672974249253), y[9] = HFloat(1.337055190016974), y[10] = HFloat(1.197118321583588), y[11] = HFloat(1.2806748535813481), y[12] = HFloat(1.716047514072701), y[13] = HFloat(1.8644880939143764), y[14] = HFloat(1.4257990737451316), y[15] = HFloat(0.9690806695246915)]]

(5.1)

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

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

boundingTriangle := plottools:-polygon(pts, thickness = 5, color = "LightGray")

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

NULL