$ontext
# Maximum area for unit-diameter polygon (variation of vanderbei's model)
# one point fixed at (0,0), others constrained in ([epsilon,1],[-1,1])
# note that 0th and Nth point are at (0,0)
# Since we bound the x coordinates away from 0,
# vanderbei's ordering constraint is equivalent to my idea of
# y_i/x_i < y_{i+1}/x_{i+1} for all pairs (i,i+1) where neither = anchor

$offtext

option seed = 105;
$if not setglobal N $setglobal N 20
set i /0*%N%/;
alias (i,j);

scalar pi;  pi= 4*arctan(1.0);
scalar epsilon; epsilon = 1e-4;

variables x(i), y(i), area;
* place start points approximately on circle, radius 1/2
x.l(i) = .4 * cos(2*pi*(ord(i)-1)/%N% + pi) + uniform(0.45,0.55);
y.l(i) = .4 * sin(2*pi*(ord(i)-1)/%N% + pi) + uniform(0,0.1);

x.lo(i) = epsilon; x.up(i) = 1;
y.lo(i) = -1; y.up(i) = 1;

set anchor(i) /0,%N%/;
x.fx(anchor) = 0;
y.fx(anchor) = 0;
set vertpairs(i,j);
vertpairs(i,j) = yes$(ord(j) gt ord(i) and (not anchor(j)));

equation obj, ordered(i), diam(i,j);

* use cross product to evaluate area of triangles (see e.g. Shaum's vector anal)
obj..
  area =e= 0.5 * sum(i$(not sameas(i,'%N%')), x(i)*y(i+1) - x(i+1)*y(i));

* y(i)/x(i) is ordered
ordered(i)$((not anchor(i)) and (not anchor(i+1)))..
  x(i)*y(i+1) - x(i+1)*y(i) =g= 0;

diam(i,j)$(vertpairs(i,j))..
  sqr(x(i) - x(j)) + sqr(y(i) - y(j)) =l= 1;

model pgon /obj,ordered,diam/;
pgon.holdfixed = 1;

* option nlp=snopt;
option nlp=conopt;
* option nlp=minos5;
solve pgon using nlp max area;

option area:8;
display area.l;
display x.l;
display y.l;

parameter theta(i);
if {(pgon.modelstat <= 2),
  theta('0') = -INF;
  theta(i)$[not anchor(i)] = arctan(y.l(i)/x.l(i));
  theta('%N%') = +INF;
  display theta;
  loop {i$[not anchor(i)],
    abort$[theta(i) > theta(i+1)] "unordered!";
  }
};