param N := 6;

var x{1..N};
param xmu1 := 0.7853981633;
param xmu2 := 2.356194491;
param pi := 3.141592654;
param xinc := 2*pi/30;

param phi{i in 1..31} := xinc*(i-1);
param x1{i in 1..31} := 0.4 + sin((2*pi)*((pi-phi[i])/(2*pi)-0.16));
param y1{i in 1..31} := 2 + 0.9*sin(pi-phi[i]);

var m{i in 1..31} = 2*x[1]*x[3]*sin(phi[i]);
var l{i in 1..31} = 2*x[3]*x[4]-2*x[1]*x[3]*cos(phi[i]);
var k{i in 1..31} = x[1]*x[1] - x[2]*x[2] + x[3]*x[3] + x[4]*x[4] - 2*x[4]*x[1]*cos(phi[i]);
var a{i in 1..31} = l[i]*l[i] + m[i]*m[i];
var b{i in 1..31} = 2*k[i]*l[i];
var c{i in 1..31} = k[i]*k[i] - m[i]*m[i];
var term{i in 1..31} = if (pi-phi[i] >= 0) then sqrt(abs(b[i]*b[i]-4*a[i]*c[i])) else -1*sqrt(abs(b[i]*b[i]-4*a[i]*c[i]));
var coss{i in 1..31} = (-b[i]+term[i])/(2*a[i]);
var sins{i in 1..31} = sqrt(abs(1-coss[i]*coss[i]));
var cosy{i in 1..31} = (x[4]+x[3]*coss[i]-x[1]*cos(phi[i]))/x[2];
var siny{i in 1..31} = (x[3]*sins[i]-x[1]*sin(phi[i]))/x[2];
var x1a{i in 1..31} = x[1]*cos(phi[i]) + x[5]*cosy[i] - x[6]*siny[i];
var y1a{i in 1..31} = x[1]*sin(phi[i]) + x[5]*siny[i] + x[6]*cosy[i];

minimize f:
sqrt( (sum {i in 1..31} ((x1a[i]-x1[i])^2 + (y1a[i]-y1[i])^2))/31.0);

subject to cons1:
-x[1] + x[2] + x[3] - x[4] >= 0;

subject to cons2:
-x[1] - x[2] + x[3] + x[4] >= 0;

subject to cons3:
-x[2]*x[2] - x[3]*x[3] + (x[4]-x[1])*(x[4]-x[1]) + 2*x[2]*x[2]*cos(xmu1) >= 0;

subject to cons4:
x[2]*x[2] + x[3]*x[3] - (x[4]+x[1])*(x[4]+x[1]) - 2*x[2]*x[2]*cos(xmu2) >= 0;

subject to cons5:
0.5 <= x[1] <= 3;

subject to cons6:
0 <= x[2];

subject to cons7:
0 <= x[3];

subject to cons8:
2 <= x[4] <= 10;

data;
var x:=
1	1
2	4.5
3	4
4	5
5	3
6	3;
solve;

display x;