Looking at the results

Within Mosel, the values of the variables and named constraints can be obtained using the getsol, getslack and similar functions. A simple report lists just the area and the positions of the vertices:

writeln("Area = ", getobjval)
forall (i in 1..N-1) do
  writeln("V", i, ": r=", getsol(rho(i)), " theta=", getsol(theta(i)))
end-do

This produces the following result for the case N=5:

Area = 0.657166
V1: r=0.616416 theta=0.703301
V2: r=1 theta=1.33111
V3: r=1 theta=1.96079
V4: r=0.620439 theta=2.58648