Solving PuzzlOR "Electrifying" puzzle with constraint logic programming

PuzzlOR's problem for December 2014–February 2015 was "Electrifying": http://puzzlor.com/2014-12_Electrifying.html

The puzzle asks to find the best placement of the 3 generators to power all the houses. Each house connects to one (the nearest) of the generators, and the cost of the connection is proportional to the Euclidean distance between the house and the generator. The best placement is the one that minimizes the total cost.

Recently Isaac Slavitt blogged about how to solve the puzzle with brute force and simulated annealing, and Jean-Francois Puget demonstrated an integer programming approach with CPLEX.

I want to show how the puzzle can be solved with constraint logic programming and ECLiPSe CLP.

Constraint logic programming is a logic programming extension that includes concepts from constraint satisfaction.

ECLiPSe CLP is an open-source Prolog-based system with a good support for modeling and solving problems with constraint logic programming. Constraint programming is not particularly well-suited for solving the "Electrifying" puzzle because the distances are reals, not integers. Some constraint programming systems and libraries can't work with reals, but ECLiPSe's library ic supports real interval arithmetic.

Here is my complete program to solve the puzzle using ECLiPSe (https://github.com/kit1980/sdymchenko-com/blob/master/electrifying-eclipse/electrifying.ecl):

:- lib(ic).
:- lib(branch_and_bound).

solve(HouseXs, HouseYs, K, GenXs, GenYs, Cost) :-
    dim(HouseXs, [N]),
    MaxX #= max(HouseXs), MaxY #= max(HouseYs),
    dim(GenXs, [K]), dim(GenYs, [K]),
    GenXs :: 1..MaxX, GenYs :: 1..MaxY,
    ( for(I, 1, N), foreach(Di, Distances), param(HouseXs, HouseYs, GenXs, GenYs, K) do
        ( for(J, 1, K), fromto(1.0Inf, Dprev, Dcurr, Di), param(I, HouseXs, HouseYs, GenXs, GenYs) do
            Dcurr $= min(Dprev, sqrt(sqr(HouseXs[I] - GenXs[J]) + sqr(HouseYs[I] - GenYs[J])))
        )
    ),
    Cost $= sum(Distances),
    bb_min(labeling([](GenXs, GenYs)), Cost, bb_options{delta:0.1}).

main :-
    HouseXs = [](2, 2, 3, 3, 4, 4,  5, 5, 5, 5, 6, 6, 7, 7, 8, 8, 8, 9, 9, 10),
    HouseYs = [](2, 5, 6, 8, 2, 10, 2, 3, 6, 8, 1, 5, 2, 8, 5, 7, 8, 4, 7, 3),
    solve(HouseXs, HouseYs, 3, GenXs, GenYs, Cost),
    ( foreacharg(X, GenXs), foreacharg(Y, GenYs) do
        LetterCode is 0'A + Y - 1,
        char_code(Letter, LetterCode),
        write(X), write(Letter), nl
    ).

The heart of the program is the solve predicate. It has 3 input parameters: HouseXs, HouseYs – arrays of X and Y coordinates of the houses, and K – the number of generators. The 3 output parameters are GenXs, GenYs, and Cost – the X and Y coordinates of the generators and the cost of the whole system. N is the number of houses. MaxX and MaxY are used to constraint the maximum coordinates of the generators.

There are two nested ECLiPSe loop constructs. The param parts just define what variables are non-local to the loop.

In the outer loop, for(I, 1, N), foreach(Di, Distances) means that as I iterates from 1 to N, Di are collected to the Distances – the list of distances from every house to its nearest generator.

The inner loop iterates over generators from 1 to K to define Di as the smallest of all Euclidean distances from I-th house to J-th generator. fromto(1.0Inf, Dprev, Dcurr, Di) is similar to the following imperative sequence: initialize Dprev with infinity, on each loop iteration define Dcurr as minimum of Dprev and the distance expression and then redefine Dprev as Dcurr, and after the loop define Di as Dcurr.

See http://eclipseclp.org/doc/bips/kernel/control/do-2.html for more info on ECLiPSe loops.

Note that Distances initially contains not numbers, but expressions for the constraint solver, because the optimal placement of the generators in not known yet. After the search finishes, all Distances elements are instantiated to concrete real numbers.

Cost is constrained to be the sum of all the values in Distances.

The last line of the solve predicate uses bb_min from the branch and bound library to find the minimum value of the Cost that can be obtained by labeling (assigning concrete values) to the X and Y coordinated of the generators (GenXs and GenYs). bb_min doesn't guarantee that the value it finds is the actual minimum, but it guarantees that there doesn't exist a solution that costs delta or more less; the value of 0.1 for delta is low enough to obtain the optimal generator placement for the given problem instance.

The main predicate defines the coordinates of the houses, runs solve, and outputs the results.

Here are the last several lines of the program's output:

Found a solution with cost 36.818603706820426__36.818603706820511
Found a solution with cost 36.407112649913259__36.407112649913365
Found a solution with cost 35.62746370694807__35.627463706948156
Found no solution with cost 0.0 .. 35.527463706948154
5B
4G
8F

The cost and coordinates of the generators are the same as in Isaac Slavitt's and Jean-Francois Puget's blog posts.

The program is not very efficient: the running time on my laptop is about 23 seconds. There are several ways to make it faster:

• use search heuristics (i.e. change labeling([](GenXs, GenYs) to bb_min(search([](GenXs, GenYs), 0, most_constrained, indomain_middle, complete, []) – this improves the time to 19 seconds)
• precompute the distances between all possible pairs of the grid cells – the constraint solver should do much better work in this case
• use the precomputed distances to formulate the puzzle as an integer linear programingm problem using ECLiPSe's eplex library.