This would be fairly straightforward to formulate in CVX, except for the constraint
u_i * [(x_i-x)^2+(y_i-y)^2] <= R
which is not DCP-compliant. I believe that can be worked around with a big M approach. Choose a numerical value of M which is big enough so that when M is added to the right-hand side, the constraint must always be satisfied for any plausible value of optimal x and y, so it would be related to the maximum distance between the input points. I leave the details of making sure your M is big enough to you.
This is not actually an MILP, because the constraint is nonloinear, but it is a mixed integer Second Order Come Problem, or mixed integer quadratically constrained linear program.
I’ll change the notation a bit from yours, and let x and y be the input vectors containing the 1st and 2nd coordinates respectively for each input data point. I have then made x_cenrter and y_center to be CVX (optimization) variables constituting the optimal center. Also note that I have changed R to be the radius, rather than the square of the radius as you had used it.
Note that as you have defined the problem, your constraints involving xmin, xmax, ymin,ymax seem like they shouldn’t be part of your problem. I have included them should you actually want to impose these constraints, but you can remove them, which I think would then solve your originally stated problem. If you do impose the xmin, xmax, ymin, ymax constraints with values which actually constrain the otherwise optimal solution, then they may need to be taken into account when determining what value of M is large enough,.
I have not actually tested this, because I do not have an integer solver installed in CVX, but I believe CVX will accept the following code, and should solve it, given enough run-time and memory, and presuming you don’t ruin the “numerics” by choosing an excessively large value of M. Nonetheless, I am not providing any guarantees, and leave it to you to check the correctness and determine the value of M to use.
cvx_begin
variable u(n) binary
variables x_center y_center
maximize(sum(u))
for i=1:n
norm([x(i) y(i)] - [x_center y_center]) <= R + M * (1-u(i))
end
xmin <= x_center <= xmax
ymin <= y_center <= ymax
cvx_end
Edit: It appears that mcg provided essentially the same solution to mine a few seconds before me. Nevertheless, I provided some complementary details to those he provided.