knp.gms : Kissing Number Problem using Variable Neighborhood Search

**Description**

Determining the maximum number of k-dimensional spheres of radius r that can be adjacent to a central sphere of radius r is known as the Kissing Number Problem (KNP). The problem has been solved for 2 (6), 3 (12) and very recently for 4 (24) dimensions. Here is a nonlinear (nonconvex) mathematical programming model known as the distance formulation for the solution of the KNP. We solve the problem by using the Variable Neighbourhood Search Algorithm. http://en.wikipedia.org/wiki/Kissing_number_problem Kucherenko, S, Belotti, P, Liberti, L, and Maculan, N, New formulations for the Kissing Number Problem. Discrete Applied Mathematics, 155:14, 1837--1841, 2007. http://doi.org/10.1016/j.dam.2006.05.012 Keywords: nonlinear programming, kissing number problem, variable neighborhood search, global optimization, sphere packing, mathematics

**Reference**

- Kucherenko, S, Belotti, P, Liberti, L, and Maculan, N, New formulations for the Kissing Number Problem. Discrete Applied Mathematics 155, 14 (2007), 1837-1841.

**Large Model of Type :** NLP

**Category :** GAMS Model library

**Main file :** knp.gms

```
$title Kissing Number Problem using Variable Neighborhood Search (KNP,SEQ=321)
$onText
Determining the maximum number of k-dimensional spheres of radius
r that can be adjacent to a central sphere of radius r is known
as the Kissing Number Problem (KNP). The problem has been solved
for 2 (6), 3 (12) and very recently for 4 (24) dimensions. Here
is a nonlinear (nonconvex) mathematical programming model known
as the distance formulation for the solution of the KNP. We solve
the problem by using the Variable Neighbourhood Search Algorithm.
http://en.wikipedia.org/wiki/Kissing_number_problem
Kucherenko, S, Belotti, P, Liberti, L, and Maculan, N,
New formulations for the Kissing Number Problem.
Discrete Applied Mathematics, 155:14, 1837--1841, 2007.
http://doi.org/10.1016/j.dam.2006.05.012
Keywords: nonlinear programming, kissing number problem, variable neighborhood search,
global optimization, sphere packing, mathematics
$offText
$eolCom //
$if not set dim $set dim 4
$if not set nspheres $set nspheres 24
Set
k 'dimension' / k1*k%dim% /
i 'spheres' / s1*s%nspheres% /;
Alias (i,j);
Variable
x(i,k) 'position of the center of the i-th sphere around the central sphere'
z 'feasibility indicator';
Equation
eq1(i) 'sphere centers have distance 2 from the center sphere'
eq2(i,j) 'minimum pairwise sphere separation distance';
eq1(i).. sum(k, sqr(x(i,k))) =e= 4;
eq2(i,j)$(ord(i) < ord(j)).. sum(k, sqr(x(i,k) - x(j,k))) =g= 4*z;
Model kissing / all /;
Scalar
lo / -2 /
up / 2 /;
x.lo(i,k) = lo;
x.up(i,k) = up;
x.l(i,k) = uniform(lo,up);
Parameter
p(i,k) 'center points of best solution'
bestobj 'feasibility indicator of best solution' / 0 /
bestbnd 'best bound on optimal value' / inf /
maxnk 'major iteration limit (search space)' / 20 /
maxns 'minor iteration limit (random starts)' / 5 /
nk 'major iteration' / 1 /
ns 'minor iteration';
kissing.solveLink = %solveLink.CallScript%;
solve kissing max z using nlp;
* Store solution as best solution
if(kissing.modelStat = %modelStat.locallyOptimal% or
kissing.modelStat = %modelStat.optimal% or
kissing.modelStat = %modelStat.feasibleSolution%,
bestobj = z.l;
p(i,k) = x.l(i,k);
else
* Do not start VNS, if we have no solution
maxnk = 0;
);
* Store dual bound, if available
bestbnd$(kissing.objEst <> na) = min(bestbnd, kissing.objEst);
* Variable Neighborhood Search Algorithm
option solPrint = off, limRow = 0, limCol = 0;
while(nk <= maxnk and bestobj < 1 and bestbnd >= 1 and kissing.solveStat <> %solveStat.userInterrupt%,
ns = 1;
repeat
// Sample a new point in the neighborhood of best point
x.l(i,k) = uniform(p(i,k) - nk*(p(i,k) - lo)/maxnk, p(i,k) + nk*(up - p(i,k))/maxnk);
solve kissing max z using nlp;
// in case we have no solution make sure z.l is small enough to avoid update of bestobj
z.l$(kissing.modelStat <> %modelStat.optimal% and
kissing.modelStat <> %modelStat.feasibleSolution% and
kissing.modelStat <> %modelStat.locallyOptimal%) = bestobj - 1;
// update dual bound
bestbnd$(kissing.objEst <> na) = min(bestbnd,kissing.objEst);
ns = ns + 1;
until(ns = maxns + 1) or (z.l > bestobj) or (bestbnd < 1) or (kissing.solveStat = %solveStat.userInterrupt%);
if(z.l <= bestobj,
// enlarge neighborhood and do minor iterations again
nk = nk + 1;
else
// update best bound, recenter neighborhood, and restart with small neighborhood
bestobj = z.l;
p(i,k) = x.l(i,k);
nk = 1;
);
display bestbnd, bestobj;
);
if(bestobj >= 1,
display 'KNP(%dim%) >= %nspheres%';
elseIf bestbnd < 1,
display 'KNP(%dim%) < %nspheres%';
elseIf nk > maxnk,
display 'Most likely: KNP(%dim%) < %nspheres%';
elseIf maxnk = 0,
display 'Could not solve initial NLP';
else
display 'VNS was interrupted';
);
```