LPL - A Mathematical Modeling Language

model snake "6-Snakes Problem";
  set i,j    "set of notches";
      e{i,j} "i and j are neighours";
      h      "set of snakes";
      k      "Length of snakes (ordered set of position)";
  parameter J{h}   "6 starting positions (heads of snake)";
  binary variable z{i,h,k};
  constraint
    A{h}: z[J,h,1] = 1;
    B{i}: sum{h,k} z = 1;
    C{h,k}: sum{i} z = 1;
    --D{i,h,k|k<#k}: sum{j in e} z[j,h,k+1] >= z[i,h,k]; //-Sl(1,1000);
    D{i,h,k|k<#k}: z[i,h,k] -> or{j in e} z[j,h,k+1];
  solve;

model data;
    set p:=1..14 "vertical lines";
        q:=1..12 "horizontal lines";    
    parameter a{p}:=Min(6,9.5-Abs(7.5-p)) "number of notches at p";
      n:=sum{p} a "total number of notches"; 
      b{p}:=[3 2 2 1 1 0 0 1 1 0 0 1 3 4] "first y-position of the notches at p";
      X{p};
    {p|p>1} (X[p]:=X[p-1] + if(p%2,2/Sqrt(3),1/Sqrt(3)));
    
    i:=1..n;  
    parameter c:=1; d; x{i} "x-coordinate"; y{i} "y-coordinate";
    {p} (d:=1, while(d<=a,(x[c]:=X,y[c]:=b+2*d,d:=d+1,c:=c+1)));
    e{i,j|i<>j}:=Sqrt((x[i]-x[j])^2+(y[i]-y[j])^2) <= 1.2;
    
    h:=[1..6]; 
    --J{h}:=[1 2 3 70 71 72];  //no solution? , is infeasible 
    J{h}:=[1 29 34 39 44 50]; //ok
    --J{h}:=[1 29 34 39 44 50];
    k:=1..12;
  end
  
  model output friend data;
    parameter PI:=3.14159;
              xa{h}:=x[39]-.3+7.6*Sin(PI/6+2*PI*(h-1)/6);
              ya{h}:=y[39]+.5+7.6*Cos(PI/6+2*PI*(h-1)/6); 
              W{i,h}:=exist{k} z;
              Z{h,k}:=argmax{i} z;
    Draw.Scale(30,30);
    Draw.DefFont('Verdana',8);
    --Draw.DefLine('line.dash',15,3);       //draw grid
    --{q} Draw.Line('#dash',-1,q+1,#q,q+1);
    --{p} Draw.Line('#dash',X,0,X,#p);
    --Draw.Text('(0,0)',-.75,1.9);
    {e[i,j]} Draw.Line(x[i],y[i],x[j],y[j]);
--    {h,k|k<#k} Draw.Line(x[Z[h,k]],y[Z[h,k]],x[Z[h,k+1]],y[Z[h,k+1]],0,5,1);
    {i} Draw.Circle(i&'',x,y,.3,1,0);
    {h} Draw.Circle(J&'',x[J],y[J],.5,h+2,0);
    {i,h|J<>i and W} Draw.Circle(i&'',x,y,.3,h+2,0);
    {h} Draw.Line(xa,ya,if(h=#h,xa[1],xa[h+1]),if(h=#h,ya[1],ya[h+1]),0,3);
  end

end

Mathematical Modeling with LPL : Explain and run a Model