This example uses the Convex-Concave procedure to solve the problem of packing \(N\) circles of equal radius into a square with side length \(l\).

This example relies on the cvxpy parser which is available here and the ECOS solver which is packaged with it.

The details of the problem can be found in Variations and Extensions of the Convex-Concave Procedure

You can download the notebook here.

In [1]: 
import cvxpy as cp, numpy as np, cvxopt, matplotlib.pyplot as plt, pickle
from IPython import display
In [2]: 
N = 41 # number of circles
l = 10 # side length of square

#algorithm parameters
delta = 1e-3 # end condition tolerance
Tau_inc = 1.5 # amount to increase tau
Tau_init = .005 # initial tau value
Tau_max = 1e6 # maximum tau value
max_iter = 200 # value to limit the maximum number of iterations attempted
np.random.seed(0)
In [3]: 
#define variables.
r = cp.Variable(1);
p = cp.Variable(2,N);
slack = cp.Variable(N*(N-1)/2)
slack_new = cp.Variable(2,2*N);

# variables for plotting.
circ = np.linspace(0,2*pi)
x_border = [0, l, l, 0, 0]
y_border = [0, 0, l, l, 0]
pi = np.pi
In [4]: 
prev_val = np.infty;
Tau = Tau_init
r_c = np.random.uniform(1)
p_c = l*np.random.uniform(size=(2,N))
for iteration in xrange(max_iter):   
    objective = cp.Minimize(-r+Tau*cp.sum_entries(slack)+Tau*cp.sum_entries(cp.sum_entries(slack_new)))
    constraints = []
    for i in xrange(N):
        constraints  += [r <= p[:,i]+slack_new[:,i], p[:,i] <= l-r+slack_new[:,i+N]]
    for i in xrange(N-1):
        for j in xrange(i+1,N):
            grad = np.matrix([[p_c[0,i]-p_c[0,j]],[p_c[1,i]-p_c[1,j]],[p_c[0,j]-p_c[0,i]],[p_c[1,j]-p_c[1,i]]])
            curval = np.linalg.norm(p_c[:,i]-p_c[:,j])
            curval *= curval
            constraints += [ cp.square(2*r)-curval-2*grad[0]*(p[0,i]-p_c[0,i])-2*grad[1]*(p[1,i]-p_c[1,i])-2*grad[2]*(p[0,j]-p_c[0,j])-2*grad[3]*(p[1,j]-p_c[1,j]) <= slack[((N-1)+(N-i))*i/2+j-i-1]]
    constraints += [slack >= 0]
    constraints += [slack_new >= 0]
    prob = cp.Problem(objective,constraints)
    result = prob.solve(solver=cp.ECOS,verbose=False)
    #solver error detected
    if(result == 'solver_error'):
        break
        
    #Plotting Code    
    #determine if constraints are violated, and if they are indicate that a circle should be red
    colors = [0]*N;        
    for i in xrange(N-1):
        if(slack_new.value[0,i] > delta or slack_new.value[0,N+i] > delta or slack_new.value[1,i] > delta or slack_new.value[1,N+i] > delta):
            colors[i] = 1
        for j in xrange(i+1,N):
            if(slack.value[((N-1)+(N-i))*i/2+j-i-1] > 1e-3):
                colors[i] = 1
                colors[j] = 1
        
    #plot the circles
    for i in xrange(N):
        if(colors[i] == 0):
            plt.plot(p.value[0,i]+r.value*np.cos(circ),p.value[1,i]+r.value*np.sin(circ),'b')
        else:
            plt.plot(p.value[0,i]+r.value*np.cos(circ),p.value[1,i]+r.value*np.sin(circ),'r')
    plt.plot(x_border,y_border,'g');
    plt.axes().set_aspect('equal')
    title = 'Iteration ' + repr(iteration) +' Tau is '+str.format('{0:.3f}',Tau)+'\n objective value is '+str.format('{0:.2f}',r.value*r.value*pi*N)+'\n violation is '+str.format('{0:.2f}',np.sum(slack.value)+np.sum(slack_new.value))
    plt.title(title)
    display.clear_output()
    plt.show();
    p_c = p.value
    r_c = r.value
    Tau = min(Tau*Tau_inc,Tau_max)
    if(np.abs(prev_val-result) <= delta and (np.sum(slack.value) < delta or Tau == Tau_max)):
            break
    prev_val = result