The problem of covering a region of the plane with a fixed number of minimum-radius identical balls is studied in the present work. An explicit construction of bi-Lipschitz mappings is provided to model small perturbations of the union of balls. This allows us to obtain analytical expressions for first- and second-order derivatives of the cost functional using nonsmooth shape optimization techniques under appropriate regularity assumptions. For regions defined as the union of disjoint convex polygons, algorithms based on Voronoi diagrams that do not rely on approximations are given to compute the derivatives. Extensive numerical experiments illustrate the capabilities and limitations of the introduced approach.