$title Max Min Location of Points in Unit Square (MAXMIN,seq=263)
$onText
This test problem locates points in the unit square such that the
distance between any two points is maximized. For certain number
of points we know optimal arrangements. This knowledge is also
used to find a lower bound on the objective by looking for perfect
square arrangements (suggested by S Dirkse).
Several formulations are presented which serve as good examples to
investigate the performance of different solution approaches. The
problem was originally proposed by Dick van Hertog and has been implemented
by Janos Pinter and used extensively by LGO with 13 and 20 points.
E. Stinstra, D. den Hertog, H.P. Stehouwer, A. Vestjens,
Constrained Maximin Designs for Computer Experiments,
Technometrics, 2002. (under revision)
Janos Pinter, LGO - Users Guide, Pinter Consulting Services, Halifax,
Canada, 2003.
Keywords: nonlinear programming, discontinuous derivatives, mathematics,
maximized minimum distance, circle packing in a square
$offText
$eolCom //
$if not set points $set points 13
Set
d 'dimension of space' / x, y /
n 'number of points' / p1*p%points% /
low(n,n) 'lower triangle';
Alias (n,nn);
low(n,nn) = ord(n) > ord(nn);
Variable
point(n,d) 'coordinates of points'
dist(n,n) 'distance between all points'
mindist;
Equation
defdist(n,n) 'distance definitions'
mindist1(n,n) 'minimum distance formulation 1'
mindist1a(n,n) 'minimum distance formulation 1 without dist'
mindist2 'minimum distance formulation 2'
mindist2a 'minimum distance formulation 2 without dist';
defdist(low(n,nn)).. dist(low) =e= sqrt(sum(d, sqr(point(n,d) - point(nn,d))));
mindist1(low).. mindist =l= dist(low);
mindist1a(low(n,nn)).. mindist =l= sqrt(sum(d, sqr(point(n,d) - point(nn,d))));
mindist2.. mindist =e= smin(low, dist(low));
mindist2a.. mindist =e= smin(low(n,nn), sqrt(sum(d, sqr(point(n,d) - point(nn,d)))));
Model
maxmin1 / defdist, mindist1 /
maxmin2 / defdist, mindist2 /
maxmin1a / mindist1a /
maxmin2a / mindist2a /;
Scalar p; // Pinter's
p = 0;
loop((n,d), // original
p = round(mod(p,10)) + 1; // nominal
point.l(n,d) = p/10; // point 0.1,.2, ... 1.0, 0.1, ...
);
point.lo(n,d) = 0;
point.up(n,d) = 1;
point.l (n,d) = uniform(0,1);
dist.l(low(n,nn)) = sqrt(sqr(point.l(n,'x') - point.l(nn,'x')) + sqr(point.l(n,'y') - point.l(nn,'y')));
point.fx('p1',d) = 0; // fix one point
Parameter bnd 'lower bound on objective';
bnd = 1/max(ceil(sqrt(card(n)))-1,1);
display bnd;
option limCol = 0, limRow = 0;
if(card(n) > 9, option solPrint = off;);
* for experimentation we will combine different model version
* with different bounds and starting points
*
* dist.lo(low) = -inf;
* dist.lo(low) = 0;
* dist.lo(low) = 0.01;
* dist.lo(low) = bnd/2;
* dist.lo(low) = bnd;
*
* solve maxmin1 max mindist using nlp;
* solve maxmin1a max mindist using nlp;
* solve maxmin2 max mindist using dnlp;
* solve maxmin2a max mindist using dnlp;
* maxmin2 and maxmin2a without bounds are well suited for LGO
* maxmin1a with bounds is well suited for conopt3 (bounds 200 point is ok)
solve maxmin1a max mindist using dnlp;
display bnd,mindist.l, point.l;