Delete minimum amount of overlapping circle?

Is there any way to delete the minimum amount of overlapping circle to get maximum amount of non-overlapping circle?

Illustration of problem:

I want to delete the overlapping circle, but keep the result maximum. For the problem above, the result should be this:

I've tried to hard code it by creating each circle combination then calculate number of intersects and disjoints. However, this takes too long. Is there's any quick solution for this kind of problem?

Please try this it should work ;) You need to install pulp. The script is divided into 3 parts:

• parametrization of your problem
• random generation of circles (I made random radii, but you can use fixed radii if necessary)
• solution with pulp

from pulp import LpMaximize, LpProblem, LpStatus, lpSum, LpVariable
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import rand

#%% PARAMETRIZE the PROBLEM
N = 10 # number of circles
X = 15 # size on X
Y = 15 # size on Y

Rmin = 1 # min radius
Rmax = 2 # max radius

#%% GENERATE RANDOM CIRCLES

cx = rand(N)*(X-2*Rmax) + Rmax
cy = rand(N)*(Y-2*Rmax) + Rmax
r  = rand(N)*(Rmax-Rmin) + Rmin

plt.figure(1)
plt.clf()
for i in range(N): plt.gca().add_artist(plt.Circle((cx[i], cy[i]), r[i], alpha=0.7))
plt.axis('image')
plt.xlim(0,X)
plt.ylim(0,Y)


#%% GET SOLUTION
model = LpProblem(name="small-problem", sense=LpMaximize)
var = [LpVariable(name="x%u"%i, lowBound=0,upBound=1,cat='Integer') for i in range(N)]

# Add objective function to model
model += lpSum(var)

# Add constraints
for i in range(N):
    for j in range(i+1,N):
        dij = np.sqrt((cx[i]-cx[j])**2 + (cy[i]-cy[j])**2)
        if dij < (r[i]+r[j]):
            model += (var[i]+var[j]<=1,"intersec %u - %u"%(i,j))


# Solve it
status = model.solve()
print(f"status: {model.status}, {LpStatus[model.status]}")
print(f"objective: {model.objective.value()}")
for v in model.variables():
    print(f"{v.name}: {v.value()}")
    if v.value():
        i = int(v.name[1])
        plt.gca().add_artist(plt.Circle((cx[i], cy[i]), r[i],fc='none',ec='C1',lw=2))