# 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`.

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.