The distance function used here is planar distance, and motion from a point to the host is straight-line motion, not the shortest path inside the polygon. This distinction is relevant only for nonconvex polygons, and can lead to disconnected Voronoi regions, as in the second snapshot. For the continuous version of the problem for the square, it is not known what the optimal solution is, even when

is, say, 7. Using ILP, as done here, and perhaps following up the computation with some local optimization, allows one to investigate possible solutions. But the big picture is known in that, under certain conditions, the optimal solution consists almost entirely of regular hexagons, as discussed in [2]. That paper also shows that the general problem is NP-complete. For more information on the continuous version, see [1].
To set up the problem as an ILP, variables

and

are introduced, where all subscripts run from 1 to

and

are the given points. The idea is that a
value of 1 indicates that

is a host, and an

value of 1 indicates that

is the nearest host to

. The objective function is

, where

gives the distance from

to

. The constraints are:
• all variables are integers;
•

,

;
•

;
•

, for each

.
The last constraint ensures that each point goes to one host; the next-to-last constraint ensures that a point goes to a host iff the
value for that host-point gets a value of 1. Thus the use of ≤ in the
sum constraint is a convenience; the travel distance always decreases, so the sum will in fact be exactly

.
[1] S. P. Fekete, J. S. B. Mitchell, and K. Beurer, "On the Continuous Fermat–Weber Problem,"
Operations Research, 53(1), 2005 pp. 61–76. doi:
10.1287/opre.1040.0137.
[2] C. H. Papadimitriou, "Worst-Case and Probabilistic Analysis of a Geometric Location Problem,"
SIAM Journal of Computing, 10(3), 1981 pp. 542–557. doi:
10.1137/0210040.