Find Chebyshev center of a region by optimization

Question

Let's consider a region like this:

Region

The aim is to find the "Chebyshev center" of the region. By using the polygons, one can approximate the perimeter of region by an $n$-gon ($n$ is large enough). Then by saving the coordinates of the perimeter in two vectors, say $x$ and $y$, and using Linear programming, compute the minimal-radius ball enclosing the region. But I want to try another way, an optimization!

I have two data sets, produced by MATLAB, (get them from here) and I want to find the Chebyshev center by optimization (not linear programming). From here I learnt how to find the nearest distance of an interior point from it's boundary:

X = "KP" /. imp // Flatten;
Y = "KI" /. imp // Flatten;
pts = Transpose[{X, Y}];
Nearest[pts, {3, 1.5}] 
(* {3,1.5} is a sample interior point *)

If I denote the nearest distance by $d$, the aim is to maximize $d$ (or minimize $-d$) over box inside region. (Like $x \in [1.5 \; 4.5] , y \in [0.5 \; 2]$.) The Chebyshev center will be obtained by maximizing $d$. But I can't define to Mathematica how treat with FindMinimum and Nearest.

This recent thread has a coded example of a bilevel optimization. — Daniel Lichtblau, Feb 03 '14 at 16:46
@DanielLichtblau That's really useful Q&A. But can we specify output argument for Nearest command like NMaximize? I'm not sure. — user2667048, Feb 03 '14 at 17:09
You would have to use Nearest[...][[1]] because Nearest returns a list. Also your inner function might need to be done as nearest[x_?NumberQ,y_?NumberQ] := nf[{x,y}][[1]] where nf is a NearestFunction created from your data points. Reason for this is to avoid any attempt by the outer optimization to use symbolic computations since the inner optimizer is in effect a black box. — Daniel Lichtblau, Feb 03 '14 at 17:19
Let me amend that. Could use nearest[x_?NumberQ, y_?NumberQ] := EuclideanDistance[{x, y}, nf[{x, y}][[1]]] — Daniel Lichtblau, Feb 03 '14 at 17:36
@DanielLichtblau Can I use this for nf function besides them?: nf[x_, y_] := Nearest[{pts, {x, y} }, 1.5 <= x <= 4.5, 0.5 <= y <= 2] — user2667048, Feb 03 '14 at 17:55
No. Look up Nearest. It doesn't take arguments like 1.5<=x<=4.5. General remark: Go to basic documentation before you go to MSE. — Daniel Lichtblau, Feb 03 '14 at 19:10
@DanielLichtblau Ok I found this: nf = Nearest[pts]; nearest[x_?NumberQ, y_?NumberQ] := EuclideanDistance[{x, y}, nf[{x, y}][[1]]]; FindMinimum[{-nearest[x, y], 1.5 <= x <= 4.5, 0.5 <= y <= 2}, {x, y} ]. I run the code and got result! Is this correct? If yes please answer to this question! — user2667048, Feb 03 '14 at 19:36
I don't know if it is correct. I don't know what result you got. I'm not even sure what point set you are using. The one I got from your link is not closed (has a horizontal leg and two that are not quite vertical, but nothing to close off the top). — Daniel Lichtblau, Feb 03 '14 at 19:38
@DanielLichtblau and I don't know what you mean by this sentence: 'has a horizontal leg and two that are not quite vertical, but nothing to close off the top' — user2667048, Feb 03 '14 at 19:44
It was an artifact of my not doing a good job plotting it. I now see the top. I'll post a regular response soon. — Daniel Lichtblau, Feb 03 '14 at 19:58
(1) "I want to find the chebyshev center by optimization (not linear programming)." I don't understand this statement. Linear programming is a form of optimization. (2) Never mind all that. There is a linear time algorithm to find the smallest circle containing a set of points. Just apply it to the vertices of your polygon, and the center of the resulting circle is your desired result. No need for complicated numerical methods. — , Feb 03 '14 at 20:35

score 5 · Accepted Answer · answered Feb 03 '14 at 20:16

Import and reconstruct the point set.

ptsdata = Import["/tmp/PI(1).mat", "LabeledData"];
pts = Transpose[({"KP", "KI"} /. ptsdata)[[All, 1, All]]];

Here is the region.

ListPlot[Transpose[{xcoords, ycoords}], PlotRange -> All, 
 AspectRatio -> Automatic]

enter image description here

Create the nearest function and the function we will maximize.

nf = Nearest[pts];
nearest[x_?NumberQ, y_?NumberQ] := 
 EuclideanDistance[{x, y}, nf[{x, y}][[1]]]

We add some constraints that pretty much force the point to stay on the inside. There are better ways to do this but that would take more work.

FindMaximum[{nearest[x, y], 1 <= x <= 6, 
  0 <= y <= 3, (y - 3) <= -1/3 (x - 3)^2}, {x, 3}, {y, 2}]

During evaluation of In[464]:= FindMaximum::eit: The algorithm does not converge to the tolerance of 4.806217383937354`*^-6 in 500 iterations. The best estimated solution, with feasibility residual, KKT residual, or complementary residual of {9.82931122747*10^-7,0.00262578662189,2.24578696996*10^-7}, is returned. >>

(* Out[464]= {1.42730059484, {x -> 3.36291430575, y -> 1.42730059479}} *)

NMaximize will give a slightly better result.

I tried the NMaximize, It gives better result. But when I compare the result to one obtained by Linear programming, I see a bit difference. — user2667048, Feb 03 '14 at 20:23
Does your linear programming look at the segments between polygon vertices? If so, that could alter slightly the distance-to-boundary values. The approach using Nearest will only give distance to vertices, not to segments. — Daniel Lichtblau, Feb 03 '14 at 22:36

Find Chebyshev center of a region by optimization

1 Answers1