3
$\begingroup$

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.

|improve this question
$\endgroup$
  • 2
    $\begingroup$ This recent thread has a coded example of a bilevel optimization. $\endgroup$ – Daniel Lichtblau Feb 3 '14 at 16:46
  • $\begingroup$ @DanielLichtblau That's really useful Q&A. But can we specify output argument for Nearest command like NMaximize? I'm not sure. $\endgroup$ – user2667048 Feb 3 '14 at 17:09
  • 2
    $\begingroup$ 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. $\endgroup$ – Daniel Lichtblau Feb 3 '14 at 17:19
  • $\begingroup$ Let me amend that. Could use nearest[x_?NumberQ, y_?NumberQ] := EuclideanDistance[{x, y}, nf[{x, y}][[1]]] $\endgroup$ – Daniel Lichtblau Feb 3 '14 at 17:36
  • 5
    $\begingroup$ 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. $\endgroup$ – Daniel Lichtblau Feb 3 '14 at 19:10
5
$\begingroup$

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.

|improve this answer
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – user2667048 Feb 3 '14 at 20:23
  • 3
    $\begingroup$ 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. $\endgroup$ – Daniel Lichtblau Feb 3 '14 at 22:36
  • $\begingroup$ Good interpretation. $\endgroup$ – user2667048 Feb 4 '14 at 6:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.