In[3]:= cntrs = {{1, 2}, {-1, 1}, {1, 3}};

In[4]:= cons = 
 Flatten @ 
  Table[r[i] + r[j] <= 
    Sqrt[(cntrs[[i]] - cntrs[[j]]).(cntrs[[i]] - cntrs[[j]])], {i, 
    2}, {j, i + 1, 3}]

Out[4]= {r[1] + r[2] <= Sqrt[5], r[1] + r[3] <= 1, 
 r[2] + r[3] <= 2 Sqrt[2]}

In[5]:= sln = 
 NMaximize[{r[1] + r[2] + r[3], Sequence @@ cons}, {r[1], r[2], r[3]}]

Out[5]= {3.03225, {r[1] -> 0.20382, r[2] -> 2.03225, r[3] -> 0.79618}}

In[7]:= Show[
 Graphics @ Table[Circle[cntrs[[i]], r[i]], {i, 3}] /. sln[[2]], 
 Graphics @ Table[Point[cntrs[[i]]], {i, 3}]]

It's a global optimization problem. Suppose you are given 3 centers. Your variables will be the radii of the 3 circles. Your objective function could be the sum of the 3 radii. Your constraints are the three constraints that the circles not overlap, that the sum of the radii of each pair be less than or equal to the distance between the centers of the pair of circles. You could do it with NMaximize.

The constraint that two circles not overlap is that the distance between the centers must be greater than or equal to the sum of their radii.