$title maximum area of an N-gon

$ontext

In this model, use a smarter choice of starting point (the regular polygon
centered at (0,1/2)) to improve the robustness of the formulation. In the
ngon.gms formulation, the naive starting point leads to failure for
N>8. Here we add a "buffer" between successive values of theta. Also, we
try conopt instead of minos.

It appears to work for N=12,16,20 but for N=20 the solution is
nonoptimal; the maximal area is smaller than for N=18. 

$offtext

option limrow=0, limcol=0, solprint=off;

set vertices/1*30/;
alias (vertices, i, j);

scalar N "number of vertices" /20/;

variable rho(vertices) "polar radius"
         theta(vertices) "polar angle"
         area;

scalar pi;
pi = 4.0 * arctan(1.0);

rho.lo(vertices) = 0.0;
rho.up(vertices) = 1.0;
theta.lo(vertices) = 0.0;
theta.up(vertices) = pi;

equations
        inscribe(i,j)       "ensure that each vertex pair no more than 1 apart"
        monotonic(vertices) "enforce monotonic increase in polar angles"
        objective           "area of polygon"
        mirror              "thetaN = pi-theta2";

inscribe(i,j)$(ord(i) < ord(j) and ord(j) <= N)..
        rho(i)*rho(i) + rho(j)*rho(j) - 2*rho(i)*rho(j)*cos(theta(i) - theta(j)) =l= 1;

* impose a minimal separation between the successive angles
monotonic(vertices)$(ord(vertices) < N)..
        theta(vertices) + pi/(20*N) =l= theta(vertices+1);

* ensure that thetaN = pi-theta2, as at the initial point
mirror..
         theta('2') =e= pi - sum(i$(ord(i)=N), theta(i));

objective..
        area =e=
          .5*sum(i$(ord(i)>=2 and ord(i)<=N-1), rho(i)*rho(i+1)*sin(theta(i+1)-theta(i)));


* Use better initial values: those from the regular polygon with vertex N
* at the origin
scalar xtemp, ytemp;
loop(i$(ord(i)>1),
   xtemp = .5 * cos(-pi/2 + ((ord(i)-1)/N)*2*pi);
   ytemp = .5 + .5 * sin(-pi/2 + ((ord(i)-1)/N)*2*pi);
   rho.l(i)   = sqrt(xtemp*xtemp + ytemp*ytemp);
   theta.l(i) = arctan(ytemp/xtemp);
);

* fix the first vertex position at the origin
rho.fx('1') = 0.0; theta.fx('1') = 0.0;

model ngon/inscribe,monotonic,mirror,objective/;

option nlp=conopt;
solve ngon using nlp maximizing area;

option N:0; option area:8;
display N, area.l;

* compare area to area of a circle of radius 1/2
scalar areaRatio;
areaRatio = area.l / (pi/4.0);
option areaRatio:8; display areaRatio;

display rho.l, theta.l;

* convert to Cartesian coordinates
set components/x,y/;
parameter cartesian(vertices,components);
loop(i$(ord(i) <= N),
  cartesian(i,'x') = rho.l(i) * cos(theta.l(i));
  cartesian(i,'y') = rho.l(i) * sin(theta.l(i));
);
option cartesian:8; display cartesian;