$title maximum area of an N-gon inscribed in unit circle

* model it with the origin roughly in the center of the figure
* (i.e. in the center of the inscribing circle).

option limrow=0, limcol=0;

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

scalar N "number of vertices (max 20)" /12/;

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

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

* bounds: ensure that radius is <= 1/2 for each point, and that all
* angles are between 0 and 2*pi
rho.lo(vertices) = 0.0;
rho.up(vertices) = 0.5;
theta.lo(vertices) = 0.0;
theta.up(vertices) = 2*pi;

equations
	monotonic(vertices) "enforce monotonic increase in polar angles"
	objective           "area of polygon";

monotonic(vertices)$(ord(vertices) < N)..
        theta(vertices) =l= theta(vertices+1);

* add in the area of the triangle between vertices 1 and N explicitly.
objective..
	area =e= 
	  .5*sum(i$(ord(i) <N), rho(i)*rho(i+1)*sin(theta(i+1)-theta(i))) +
          .5*sum(i$(ord(i) =N), rho('1')*rho(i)*sin(2*pi - theta(i)));

* set initial values: rho positioned on unit circle, angles equally spaced
loop(i$(ord(i)<N),
   rho.l(i)   = uniform(0,0.5);
   theta.l(i) = 2*pi*(ord(i)-1)/N*(1.0 + uniform(-0.1,0.1));
);

* fix the first vertex angle to be zero
theta.fx(i)$(ord(i) = 1) = 0.0;

model ngon/monotonic,objective/;

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

option N:0; option area:8; 
display N, area.l;
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;