How to formulate a MILP in CVX

I am new to CVX and will appreciate your help.
I have a set of points in 2d plane and a circle of some radius R. I would like to find the center of of the circle,i.e, (x,y) that includes most points. This problem can be formulated as MILP as following:

maximize sum (u_i)
subject to
u_i * [(x_i-x)^2+(y_i-y)^2]<=R for all i
xmin <= x <=xmax
ymin <= y <=ymax

where u_i is binary which indicates whether the point is included in the circle or not and " i " is index of the points.

I really do not know how to solve this problem in cvx

As I’m sure you know, MIDCP’s must still satisfy the disciplined convex programming rules, even on the binary variables. So the product between the u_i values and the distance values is problematic.

However, there is a solution. Put u_i on the right-hand side instead; something like this:

(x_i-x)^2 + (y_i-y)^2 <= R + (1 - u_i) * M

where M is chosen to be large enough so that when u_i is zero, that inequality is trivially satisfied. That said, you don’t want it to be any larger than necessary, either, to prevent any numerical issues. My suspicion is that 2*(xmax^2+ymax^2) would be enough.

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.

The problem setup is as following, I will drop n points uniformly over a square reqion.As the circle center does not have to be within the square, I dropped the xmax and y max constraint.

Mark_L_Stone, I dropped n points over a square region of 1000m by 1000m and run your code. I called MOSEK solver within CVX. I plotted the circle along with the points on square region. Is there any way to check if the solution I have is optimum, i.e., no other circle location can result in including more points?.

If CVX/MOSEK claims to have solved your problem, and M is not too small, then the solution reported by CVX should be an optimum solution (there’s no reason in general to think that the “optimal” circle center is unique though).

If the problem is infeasible, then M must be too small.

However, if the problem is feasible, it is possible that M could be too small and therefore the optimum reported by CVX not solve the problem you intended. This could happen for instance if there are 2 clusters of points, with a large separation between the clusters. If M is big enough to capture the maximum distances between points all is well. But if R + M is not big enough, the only feasible solution may be to choose the circle center in the vacant space in the middle between the clusters, potentially not capturing any points, whereas with a big enough M, the circle could be centered in one of the clusters, thereby capturing a bunch of points. Less extreme versions of a similar phenomenon, if M is too small, could result in CVX reporting a solution to your problem which does not solve the problem you intended.

I don’t know of a “direct” way of checking the correctness of the solution (relative to your intended problem, not relative to the problem you input to CVX). Perhaps someone else has some ideas on how to check.

Thank you Mark_L_Stone and Michael C. Grant for your valuable suggestions. CVX along with MOSEK solved the problem for one circle. Now, I want to extend it to multiple circles (3 circles) . The problem is as following

I have 3 circles with equal radius R and N points distributed on a 2D plane. I want to place these 3 circles such that they enclose the maximum number of points. Note that a point could be included in more than one circle. Let i be the index of points and j be the index of the circles, U_ij is a binary variable where U_ij=1 implies that point i is enclosed by circle j. The formulation I came up is
Assume we have 3 circles,

max V_i
subject to

        constraint 1)    U_ij>=V_i      for all   j=1,2,3
        constraint 2)   (x_i-x_ j)^2 + (y_i-y_ j)^2 <= R + (1 - U_ij) * M   for all i and j

where V_i is a binary variable given by V_i=min U_ij (minimization is taken over j) for example V_1=min{U_11, U_12, U_13}

The solver should give

  1. (x_ j,y_ j) for all 3 circles (center of the circles)
    
  2. Vi
    
  3. U_ij which indicates whether point i is enclosed bypoint j or not.

I think you are at the stage of just having modeling questions. You already know how to do big M using CVX’s MMIDCP capability.

I haven’t really checked what you’ve done. If a point is in more than one circle, do you want to count it only once? If so, make sure that’s what your formulation does.

yes Mark. If a point is in more than one circle, it should be counted once. That is why I introduced the binary variable V_ij. It is given by V_i=min U_ij (minimization is taken over j) for example V_1=min{U_11, U_12, U_13}. I wonder if CVX can handle this type of formulation

The problem should be accepted and solved by CVX/solver if there is enough memory and time. Therefore, it’s a question of whether your formulation is correct, to include the value of M. If M is too small, you will not get the right answer. If M winds up being big, there could be numerical challenges for the solver.

Hello Mark. How to prove that this problem is a SOC problem?. By the way, I do not have the following constraints in my modelling

xmin <= x_center <= xmax
ymin <= y_center <= ymax

If you want to learn about SOCPs, you can read http://web.stanford.edu/~boyd/cvxbook/ .

1 Like

The MOSEK modelling cook at https://mosek.com/resources/doc also has a lot of material about SOCP.

1 Like

I solved the MISOCP using cvx and mosek> I wonder what the complexity order of the algorithm.

I am working to solve the same problem, would you mind sharing the MOSEK code you used ?