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';

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';
);