param N integer := 378;
param nn1{1..5};
param nn2{1..5};
param ny11 := nn1[1];
param ny22 := nn2[1];
param nu11 := nn1[2];
param nu22 := nn2[2];
param ny1 := nn1[3];
param ny2 := nn2[3];
param nu1 := nn1[4];
param nu2 := nn2[4];
param td1 := nn1[5];
param td2 := nn2[5];
param mo := max(ny11,ny22,td1+nu11,td2+nu22,ny1,ny2,td1+nu1,td2+nu2)+1;
param theta1{1..55};
param theta2{1..55};
set I = 1 .. N;
set J = N..2*N-1;
param Nu integer :=  189;
var a1v{i in 1..Nu};
var a2v{i in 1..Nu};
param k{i in 1..N} := i;
param mx := 1e-2;
param my := 8e-3;
param pi := 4*atan(1);
param alp :=exp(-15/75);
param pn{i in 1..Nu} := i*2*pi/N;
param bdy := 8.5e-3;
var alpha{1..Nu,1..2};
param ind{0..14} integer;
param t := 1e-2;
var phi{i in 0..14,0..2*N-1,1..5} := if i==0 then 1;
var x{i in -mo..2*N-1,j in 1..2};
var y{i in -mo..2*N-1, j in 1..2} := if i >= 0 then .01;

param lim;
set S1 := -lim..0;
set S2 := 1..lim;
set S := S1 union S2;
set W := {i in S,j in S:not(i==0 and j==0)};

s.t. weyl1{(m1,m2) in W}:
sum{n in N+1..2*N-1} cos(2*pi*(m1*(y[n,1]+bdy)/(2*bdy)+m2*(y[n,2]+bdy)/(2*bdy))) <= t;
s.t. weyl2{(m1,m2) in W}:
sum{n in N+1..2*N-1} sin(2*pi*(m1*(y[n,1]+bdy)/(2*bdy)+m2*(y[n,2]+bdy)/(2*bdy))) <= t;

s.t. bound51{j in N..2*N-1, i in 1..2}: -bdy<=y[j,i]<= bdy;
s.t. xstart{i in -mo..-1, j in 1..2}: x[i,j] =0;
s.t. pstart{i in 0..2*N-1,j in 1..4}: phi[0,i,j] = 1;
s.t. ystart{i in -mo..-1, j in 1..2}: y[i,j] =0;
s.t. inp1{i in 0..N-1, j in 1..2}:
x[i,j]= if j==1 then
sum{u in 1..Nu}(a1v[u]*cos(pn[u]*k[i+1]+alpha[u,j]))
else
sum{u in 1..Nu}(a2v[u]*cos(pn[u]*k[i+1]+alpha[u,j]));
s.t. inp2{i in N..2*N-1, j in 1..2}:
x[i,j] = x[i-N,j];
s.t. p1{i in 0..2*N-1,j in 1..ny11}:
phi[1,i,j]=y[i-j,1];
s.t. p2{i in 0..2*N-1,j in 1..ny22}:
phi[2,i,j]=y[i-j,2];
s.t. p3{i in 0..2*N-1,j in td1..td1+nu11}:
phi[3,i,j-td1+1]=x[i-j,1];
s.t. p4{i in 0..2*N-1,j in td2..td2+nu22}:
phi[4,i,j-td2+1]=x[i-j,2];
s.t. p5{i in 0..2*N-1,j in 1..ny1,l in 1..j}:
phi[5,i,(j-1)*j/2+l]=y[i-j,1]*y[i-l,1];
s.t. p6{i in 0..2*N-1,j in 1..ny2,l in 1..j}:
phi[6,i,(j-1)*j/2+l]=y[i-j,2]*y[i-l,2];
s.t. p7{i in 0..2*N-1,j in td1..td1+nu1,l in td1..j}:
phi[7,i,(j-td1)*(j-td1+1)/2+l-td1+1]=x[i-j,1]*x[i-l,1];
s.t. p8{i in 0..2*N-1,j in td2..td2+nu2,l in td2..j}:
phi[8,i,(j-td2)*(j-td2+1)/2+l-td2+1]=x[i-j,2]*x[i-l,2];
s.t. p9{i in 0..2*N-1,j in 1..ny1,l in td1..td1+nu1}:
phi[9,i,(j-1)*(nu1+1)+l-td1+1]=y[i-j,1]*x[i-l,1];
s.t. p10{i in 0..2*N-1,j in 1..ny1,l in td2..td2+nu2}:
phi[10,i,(j-1)*(nu2+1)+l-td2+1]=y[i-j,1]*x[i-l,2];
s.t. p11{i in 0..2*N-1,j in 1..ny2,l in td1..td1+nu1}:
phi[11,i,(j-1)*(nu1+1)+l-td1+1]=y[i-j,2]*x[i-l,1];
s.t. p12{i in 0..2*N-1,j in 1..ny2,l in td2..td2+nu2}:
phi[12,i,(j-1)*(nu2+1)+l-td2+1]=y[i-j,2]*x[i-l,2];
s.t. p13{i in 0..2*N-1,j in 1..ny1,l in 1..ny2}:
phi[13,i,(j-1)*ny2+l]=y[i-j,1]*y[i-l,2];
s.t. p14{i in 0..2*N-1,j in td1..td1+nu1,l in td2..td2+nu2}:
phi[14,i,(j-td1)*(nu2+1)+l-td2+1]=x[i-j,1]*x[i-l,2];
s.t. out1{i in 0..2*N-1,j in 1..2}:
y[i,j]=if j==1 then
sum{kk in 0..14,l in 1..ind[kk]} phi[kk,i,l]*theta1[sum{ll in 0..kk-1}ind[ll]+l]
else
sum{kk in 0..14,l in 1..ind[kk]} phi[kk,i,l]*theta2[sum{ll in 0..kk-1}ind[ll]+l];

s.t. move1{i in 0..N-2, j in 1..2}: -mx <= x[i+1,j] - x[i,j] <= mx;
s.t. move2{j in 1..2}: -mx <= x[0,j]-x[N-1,j] <= mx;
s.t. movey{i in N-1..2*N-2, j in 1..2}: -my <= y[i+1,j] - y[i,j] <= my;

data;
param nn1 :=
1 4
2 4
3 2
4 1
5 1
;
param nn2 :=
1 4
2 4
3 2
4 1
5 1
;

param theta1 :=
 1  2.03864867E-07
 2  0.285694268
 3  0.18635982
 4  0.089657647
 5  0.194695534
 6  0.0430581642
 7  0.127531047
 8  0.10198672
 9  0.00710628085
 10 -0.468482504
 11 -0.456244308
 12 -0.277500058
 13 -0.154069349
 14 -0.025268557
 15  0.561727289
 16  0.407211533
 17  0.215854217
 18  0.109830097
 19  0.00565726923
 20 -198.447103
 21  320.067051
 22 -130.013242
 23 -65.8084213
 24  77.5109037
 25 -21.3410784
 26  0.787890389
 27  7.39848218
 28 -19.3207814
 29  1.75715408
 30  7.3164134
 31 -10.2034947
 32  7.18411879
 33 -115.095362
 34 -9.34076031
 35  93.8935778
 36 -11.6909496
 37  62.4575882
 38  13.5410976
 39 -53.5132462
 40 -1.84210589
 41  61.3262926
 42  0.549318718
 43 -36.2658647
 44  3.1189285
 45 -28.802668
 46 -0.800808306
 47  16.789681
 48  232.67841
 49 -136.86058
 50 -183.824698
 51  108.914098
 52 -2.56941865
 53 -6.23845388
 54 -8.79138994
 55  25.648625
;
param theta2 :=
 1  2.7934064E-05
 2  0.103222348
 3  0.100604488
 4  0.0228412518
 5  0.0309463861
 6  0.142987314
 7  0.163457024
 8  0.159247778
 9  0.130155891
 10 -0.365105111
 11 -0.465492823
 12 -0.267075648
 13 -0.134845498
 14 -0.0296975095
 15  0.771148851
 16  0.494423193
 17  0.260581844
 18  0.107998863
 19 -0.00979548709
 20  14.5435202
 21 -13.2441078
 22 -0.440352568
 23  5.26672539
 24 -4.46454248
 25  0.69225083
 26 -0.112379484
 27  7.44615257
 28 -4.43107291
 29 -1.58503942
 30  3.09797088
 31 -5.81463511
 32  32.0417403
 33 -5.47940057
 34 -25.8930792
 35  9.41504905
 36 -42.62413
 37  25.9621262
 38  33.5940277
 39 -24.5028357
 40 -18.0154727
 41 -5.14320009
 42  11.4737285
 43  2.45151587
 44  23.8991518
 45 -5.085284
 46 -16.0534403
 47  2.93749645
 48 -33.1134629
 49  13.6296698
 50  22.7600102
 51 -9.2894652
 52  0.818777238
 53 -3.09977152
 54 -8.49067596
 55  12.3313026
;

let ind[0] := 1;
let ind[1] := ny11;
let ind[2] := ny22;
let ind[3] := nu11+1;
let ind[4] := nu22+1;
let ind[5] := ny1*(ny1+1)/2;
let ind[6] := ny2*(ny2+1)/2;
let ind[7] := (nu1+1)*(nu1+2)/2;
let ind[8] := (nu2+1)*(nu2+2)/2;
let ind[9] := ny1*(nu1+1);
let ind[10] := ny1*(nu2+1);
let ind[11] := ny2*(nu1+1);
let ind[12] := ny2*(nu2+1);
let ind[13] := ny1*ny2;
let ind[14] := (nu1+1)*(nu2+1);

let lim := 2;
repeat{
display lim;
option solver knitro;
option solver ipopt;
option knitro_options "feasible=1";
solve;
let lim := lim+4;
} until lim>=7;