PyCSP3
COP  medium

# Problem Mario

From Amaury Ollagnier and Jean-Guillaume Fages, in the context of the 2013 Minizinc Competition: This models a routing problem based on a little example of Mario’s day. Mario is an Italian Plumber and his work is mainly to find gold in the plumbing of all the houses of the neighborhood. Mario is moving in the city using his kart that has a specified amount of fuel. Mario starts his day of work from his house and always ends to his friend Luigi’s house to have the supper. The problem here is to plan the best path for Mario in order to earn the more money with the amount of fuel of his kart. From a more general point of view, the problem is to find a path in a graph:

• path endpoints are given (from Mario’s to Luigi’s)
• the sum of weights associated to arcs in the path is restricted (fuel consumption)
• the sum of weights associated to nodes in the path has to be maximized (gold coins)’’

Finding the Best Path for Mario. Image from pngimg.com

To build a COP (Constraint Optimization Problem) model, we need first to import the library PyCSP$^3$:

from pycsp3 import *


Then, we need some data. Here, the data come from the JSON file “Mario_easy-2.json”. Note that fuels[i][j] indicates the amount of fuel needed to go from house i to house j, and golds[i] indicates the quantity of gold that can be collected at house i.

marioHouse = 0
luigiHouse = 1
fuelLimit = 2000
fuels = cp_array([
[0, 221, 274, 808, 13, 677, 670, 921, 943, 969, 13, 18, 217, 86, 322],
[0, 0, 702, 83, 813, 679, 906, 246, 335, 529, 719, 528, 451, 242, 712],
[274, 702, 0, 127, 110, 72, 835, 5, 161, 430, 686, 664, 799, 523, 73],
[808, 83, 127, 0, 717, 80, 31, 71, 683, 668, 248, 826, 916, 467, 753],
[13, 813, 110, 717, 0, 951, 593, 579, 706, 579, 101, 551, 280, 414, 294],
[677, 679, 72, 80, 951, 0, 262, 12, 138, 222, 146, 571, 907, 225, 938],
[670, 906, 835, 31, 593, 262, 0, 189, 558, 27, 287, 977, 226, 454, 501],
[921, 246, 5, 71, 579, 12, 189, 0, 504, 221, 483, 226, 38, 314, 118],
[943, 335, 161, 683, 706, 138, 558, 504, 0, 949, 393, 721, 267, 167, 420],
[969, 529, 430, 668, 579, 222, 27, 221, 949, 0, 757, 747, 980, 589, 528],
[13, 719, 686, 248, 101, 146, 287, 483, 393, 757, 0, 633, 334, 492, 859],
[18, 528, 664, 826, 551, 571, 977, 226, 721, 747, 633, 0, 33, 981, 375],
[217, 451, 799, 916, 280, 907, 226, 38, 267, 980, 334, 33, 0, 824, 491],
[86, 242, 523, 467, 414, 225, 454, 314, 167, 589, 492, 981, 824, 0, 143],
[322, 712, 73, 753, 294, 938, 501, 118, 420, 528, 859, 375, 491, 143, 0]
])
golds = [0, 0, 40, 67, 89, 50, 6, 19, 47, 68, 94, 86, 34, 14, 14]
nHouses = len(golds)


Remark. Because it is not possible to use a special object as an index of a list of integers in Python, we transform the list of integers fuels into a more specific list with the function cp_array(). This will be useful later when we shall post some constraints Element. However, when the data is directly loaded from a file (which is the usual case), all lists of integers have the specific type of list returned by cp_array(), and so, it is very rare to need to call this function explicitly.

Remark. Note that it is possible to load the data from a local JSON file with the statement:

or a distant JSON file by using an URL.

We start our COP model by introducing an array $x$ of 15 variables, one variable per house. These variables will indicate what is the circuit followed by Mario.

# x[i] is the house succeeding to the ith house (itself if not part of the route)
x = VarArray(size=nHouses, dom=range(nHouses))


We can display the structure of the array, as well as the domain of the first variable (note that all variables have the same domain).

print("Array x: ", x)
print("Domain of any variable: ", x[0].dom)

Array x:  [x[0], x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10], x[11], x[12], x[13], x[14]]
Domain of any variable:  0..14


Concerning the constraints, we start with a single simple constraint Intension: ensuring that after visiting Luigi’s house, Mario returns to his house.

satisfy(
# Mario's house succeedes to Luigi's house
x[luigiHouse] == marioHouse
);


Then, we post a constraint Circuit.

satisfy(
# Mario must make a tour (not necessarily complete)
Circuit(x)
);


We can display the internal representation of the two posted constraints; this way, although a little bit technical, we can see that the arguments to the two constraints are correct.

print(posted())

intension(function:eq(x[1],0))
circuit(list:[x[0], x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10], x[11], x[12], x[13], x[14]])


Interestingly, by calling the function solve(), we can check that the problem is satisfiable (SAT). We can also display the found solution. Here, we call the function values() that collects the values assigned to a specified list of variables.

if solve() is SAT:
print(values(x))
print("Consumed fuel: ", sum(fuels[i][x[i].value] for i in range(nHouses)))
print("Collected gold: ", sum((x[i].value != i) * golds[i] for i in range(nHouses)))

[1, 0, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]
Consumed fuel:  221
Collected gold:  0


It seems that the first found solution is correct, but not very interesting. There are many loops (when we have x[i]=i) indicating that the circuit is very short (and so, the quantity of collected gold is very small).

To incitate Mario at following a longer circuit, one can post the objective of collecting the maximal quantity of gold.

maximize(
# maximizing collected gold
Sum((x[i] != i) * golds[i] for i in range(nHouses) if golds[i] != 0)
);


We can run again the solver, with this optimization task. Note that we need to check that the status returned by the solver is now OPTIMUM.

if solve() is OPTIMUM:
print(values(x))
print("Consumed fuel: ", sum(fuels[i][x[i].value] for i in range(nHouses)))
print("Collected gold: " , bound())

[2, 0, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 1]
Consumed fuel:  7075
Collected gold:  628


The good news is that the quantity of collected gold is impressive (everything that could be collected has been collected). The bad news is that the “solution” is not really valid because we have not yet taken fuel limitation into account. Actually, the consumed fuel is really too high.

satisfy(
# we cannot consume more than the available fuel
Sum(fuels[i][x[i]] for i in range(nHouses)) <= fuelLimit
);


Technically, the constraint is a rather complex expression involving a main sum and some subexpresions that correspond to partial forms of constraints Element. This is why if we display, by calling posted(-1), the 5 first posted constraints during the last call to satisfy(), we can see that some auxiliary variables, called $aux_gb$, have been automatically introduced so as to form complete constraints Element. Note that le and eq respectively stands for ‘less than or equal to’ and ‘equal to’.

print(posted(-1)[:5])

sum(list:[aux_gb[0], aux_gb[1], aux_gb[2], aux_gb[3], aux_gb[4], aux_gb[5], aux_gb[6], aux_gb[7], aux_gb[8], aux_gb[9], aux_gb[10], aux_gb[11], aux_gb[12], aux_gb[13], aux_gb[14]], condition:(le,2000))
element(list:[0, 221, 274, 808, 13, 677, 670, 921, 943, 969, 13, 18, 217, 86, 322], index:x[0], condition:(eq,aux_gb[0]))
element(list:[0, 0, 702, 83, 813, 679, 906, 246, 335, 529, 719, 528, 451, 242, 712], index:x[1], condition:(eq,aux_gb[1]))
element(list:[274, 702, 0, 127, 110, 72, 835, 5, 161, 430, 686, 664, 799, 523, 73], index:x[2], condition:(eq,aux_gb[2]))
element(list:[808, 83, 127, 0, 717, 80, 31, 71, 683, 668, 248, 826, 916, 467, 753], index:x[3], condition:(eq,aux_gb[3]))


We can run again the solver.

if solve() is OPTIMUM:
print(values(x))
print("Consumed fuel: ", sum(fuels[i][x[i].value] for i in range(nHouses)))
print("Collected gold: " , bound())

[13, 0, 14, 1, 2, 10, 3, 5, 12, 6, 4, 7, 11, 8, 9]
Consumed fuel:  1890
Collected gold:  628


Finally, we give below the model in one piece. Here the data is expected to be given by the user (in a command line).

from pycsp3 import *

marioHouse, luigiHouse, fuelLimit, houses = data
fuels, golds = zip(*houses)  # using cp_array is not necessary since intern arrays
# have the right type (for the constraint Element)
nHouses = len(houses)

# s[i] is the house succeeding to the ith house (itself if not part of the route)
s = VarArray(size=nHouses, dom=range(nHouses))

satisfy(
# we cannot consume more than the available fuel
Sum(fuels[i][s[i]] for i in range(nHouses)) <= fuelLimit,

# Mario must make a tour (not necessarily complete)
Circuit(s),

# Mario's house succeeds to Luigi's house
s[luigiHouse] == marioHouse
)

maximize(
# maximizing collected gold
Sum((s[i] != i) * golds[i] for i in range(nHouses) if i not in {marioHouse, luigiHouse})
)