For Centroidal Voronoi Tessellations we have the following result:
Let $\rho(\mathbf{y})$ be the probability density function on the open set $\Omega \in \mathbf{R}^{\mathbf{N}}$ over which is found the sets of $k$ points $\{ \mathbf{z}_i \}_{i=1}^k$ and $k$ regions $\{ V_i \}_{i=1}^k$ that tessellate $\Omega$ . The functional $$ F( (\mathbf{z}_i,V_i,i=1,\cdots,k) )=\sum_{i=1}^{k} \int_{\mathbf{y} \in V_i} \rho(\mathbf{y}) |\mathbf{y}-\mathbf{z}_i|^2 \mathrm{d} \mathbf{y} $$ is minimized when the set over $\Omega$ is a Voronoi tessellation with $\{ \mathbf{z}_i \}_{i=1}^k$ as the centroids of the Voronoi regions $\{ V_i \}_{i=1}^k$ $$ \mathbf{z}_i=\dfrac{\int_{\mathbf{y} \in V_i}\mathbf{y}\rho(\mathbf{y})\mathrm{d} \mathbf{y}}{\int_{\mathbf{y} \in V_i}\rho(\mathbf{y})\mathrm{d} \mathbf{y}} $$
The question is, for a nonlinear target function $$ F( (\mathbf{z}_i,V_i,i=1,\cdots,k) )=\sum_{i=1}^{k} \int_{\mathbf{y} \in V_i} \rho(\mathbf{y}) e^{|\mathbf{y}-\mathbf{z}_i|}|\mathbf{y}-\mathbf{z}_i|^2 \mathrm{d} \mathbf{y} $$ how to minimize it with varying $\{\mathbf{z}_i \}_{i =1}^k$?