SeedRandom[1];
pts = RandomReal[{0, 12}, {100, 2}];
Needs["ComputationalGeometry`"];
dt = DelaunayTriangulation[pts];
dt // Column
toPairs[{m_, ns_List}] := Map[{m, #} &, ns];
edges = Flatten[Map[toPairs, dt], 1];
Graphics[GraphicsComplex[pts, {Line[edges],
  Red, PointSize[Large], Point[pts]}]]

g = GraphData[{"Heptahedral", 5}]

n = VertexCount[g];
pts = Table[PropertyValue[{g, i}, VertexCoordinates], {i, n}];
rads = Array[r, n];
xs = Array[x, n];
ys = Array[y, n];

f1 = Total[(r[#1] + r[#2] - EuclideanDistance[{x[#1], y[#1]}, {x[#2], y[#2]}])^2 & @@@ EdgeList[g]];
f2 = Total[MapThread[({#1, #2} - #3)^2 &, {xs, ys, pts}], -1];

rules = Last@NMinimize[f1 + 0.1 f2, Join[rads, xs, ys]];

circles = MapThread[Circle[{#1, #2}, #3] &, {xs, ys, rads} /. rules];
g2 = Fold[SetProperty[{##}, VertexCoordinates -> {x[#2], y[#2]} /. rules] &, g, Range[n]];

Show[Graphics[circles], g2]

Graphics`Region`RegionInit[];
dist[{a_, b_}] := Sqrt[(a - b).(a - b)];
n = 30;

(* get a random mesh *)
pts = RandomReal[1, {n, 2}];
m = DelaunayMesh[pts]["MeshObject"];

(* get edges and non-edges *)
edges = Sort /@ m["Edges"];
nonedges = Complement[Subsets[Range[n], {2}], edges];

(* variable lists *)
rads = Array[r, n];
pos = Array[{x[#], y[#]} &, n];

(* set up constraints for edges *)
cons = (Total@rads[[#]] == dist@pos[[#]]) & /@ edges;
cons = Join[cons, Thread[rads > 0]];

(* define function to maximize *)
f = Total[dist@pos[[#]] & /@ nonedges];

(* variable initializations *)
vars = Join[rads, Flatten[{pos, pts}, {2, 3}]];

(* do the maximization *)
{posv, radsv} = {pos, rads} /. Last@FindMaximum[{f, cons}, vars];

(* result *)
Graphics[{MapThread[Circle, {posv, radsv}],
  Opacity[0.5, Green], Line[posv[[#]] & /@ edges]}]