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