# Find Chebyshev center of a region by optimization

Let's consider a region like this:

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 3 '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 3 '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 3 '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 3 '14 at 17:36
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 3 '14 at 19:10

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]


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 3 '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 3 '14 at 22:36
Good interpretation. –  user2667048 Feb 4 '14 at 6:25