COP medium

Problem Cardinality Constrained Multi-cycle Problem (CCMcP)

The Cardinality Constrained Multi-cycle Problem (CCMcP) is a variation of the Kidney Exchange Problem (KEP). One can consider the CCMcP as an Asymmetric Travelling Salesman Problem (ATSP) with subtours (cycles) allowed (but of limited size $k$). Each arc has an associated weight. The arcs with negative weights cannot be part of subtours, and the objective is to maximize the sum of weights occurring along the arcs of the computed subtours (cycles). The problem is studied in “On the kidney exchange problem: cardinality constrained cycle and chain problems on directed graphs: a survey of integer programming approaches” by V. Mak-Hau, in Journal of Combinatorial Optimization 2017.

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.

weights = cp_array([
    [0, -1, 4, 2, 5, 6],
    [-1, 0, -1, 8, -1, -1],
    [3, 8, 0, -1, 9, 3],
    [5, -1, -1, 0, 6, 5],
    [4, -1, 9, -1, 0, 2],
    [7, 8, 8, 2, -1, 0]
])
n, k = len(weights), 3  # n is the number of nodes, k is the length limit of cycles (here, 3)
nMaxCycles = n//k + (1 if n%k != 0 else 0)

Note that we need to call cp_array() here on the two-dimensional list of integers (to specialize its type), because of the expression we will write later when posting the objective. However, the good news is that when data come from a file, the transformation operated by cp_array() is automatically performed.

Ultimately, we need to know how cycles are formed (what is the successor node of each node?). Hence, we introduce an array $x$ of variables. We also need a way to identify the cycles (to which cyle belongs each node?). We introduce a second array $y$ of variables.

# x[i] is the successor node of node i (in the cycle where i belongs)
x = VarArray(size=n, dom=range(n))

# y[i] is the cycle (index) where the node i belongs
y = VarArray(size=n, dom=range(nMaxCycles))

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 in x: ", x[0].dom)

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

Array x:  [x[0], x[1], x[2], x[3], x[4], x[5]]
Domain of any variable in x:  0..5
Array y:  [y[0], y[1], y[2], y[3], y[4], y[5]]
Domain of any variable in y:  0 1

To start, we need to impose that no two nodes can be succeeded by the same node (otherwise, we have not a structure in term of cycles). Also, we ensure that a node and its successor belongs to the same cycle (by means of a constraint [Element](/documentation/constraints/Element)).

satisfy(
    AllDifferent(x), 
    
    # ensuring correct cycles
    [y[i] == y[x[i]] for i in range(n)]
);

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

print(posted())

allDifferent(list:[x[0], x[1], x[2], x[3], x[4], x[5]])
element(list:[y[0], y[1], y[2], y[3], y[4], y[5]], index:x[0], condition:(eq,y[0]))
element(list:[y[0], y[1], y[2], y[3], y[4], y[5]], index:x[1], condition:(eq,y[1]))
element(list:[y[0], y[1], y[2], y[3], y[4], y[5]], index:x[2], condition:(eq,y[2]))
element(list:[y[0], y[1], y[2], y[3], y[4], y[5]], index:x[3], condition:(eq,y[3]))
element(list:[y[0], y[1], y[2], y[3], y[4], y[5]], index:x[4], condition:(eq,y[4]))
element(list:[y[0], y[1], y[2], y[3], y[4], y[5]], index:x[5], condition:(eq,y[5]))

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("x:", values(x))
    print("y:", values(y))

x: [0, 1, 2, 3, 4, 5]
y: [0, 0, 0, 0, 0, 0]

This is quite disappointing: cycles are just self-loops. Besides, all cycles have the same number (index) 0.

Let us prevent self-loops and let us avoid arcs with negative weights.

satisfy(
   # disabling infeasible arcs
   [x[i] != j for i in range(n) for j in range(n) if i == j or weights[i][j] < 0]
);

Let us run again the solver:

if solve() is SAT:
    print("x:", values(x))
    print("y:", values(y))

x: [2, 3, 0, 4, 5, 1]
y: [0, 0, 0, 0, 0, 0]

This is better. We have now two cycles (assuming that the ouptut of the solver is $[2,3,0,4,5,1]$ for x): (0,2,0) and (1,3,4,5,1). However, they have the same index, and the second one is too long.

This is the right time to limit the length of cycles with a constraint BinPacking (we indicate that the size of each arc is 1 in any bin/cycle):

satisfy(
   # each cycle contains at most k arcs
   BinPacking(y, sizes=1) <= k
);

Let us run again the solver:

if solve() is SAT:
    print("x:", values(x))
    print("y:", values(y))

x: [2, 3, 4, 5, 0, 1]
y: [1, 0, 1, 0, 1, 0]

This is far better. We have two cycles (assuming that the ouptut of the solver is $[2,3,4,5,0,1]$ for x and $[1,0,1,0,1,0]$ for y)): (0,2,4,0) and (1,3,5,1), both of size 3.

One can note that numbers associated with the cycles are intercheangeable. We don’t care if the first encountered one has number 0 and the second one number 1, or the reverse. This is why (to break some symmetries) we post a constraint Precedence on y:

satisfy(
   # tag(symmetry-breaking)
   Precedence(y)
);

Let us run again the solver:

if solve() is SAT:
    print("x:", values(x))
    print("y:", values(y))

x: [2, 3, 4, 5, 0, 1]
y: [0, 1, 0, 1, 0, 1]

One can observe that numbers used for cycles are now ordered (their first occurrences are ordered).

The last thing to do is to post the objective:

maximize(
    # maximizing the sum of arc weights of selected cycles
    Sum(weights[i][x[i]] for i in range(n))
);

We can run 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("x: ", values(x))
    print("y: ", values(y))
    print(f"objective value: {bound()}")

x:  [2, 3, 4, 5, 0, 1]
y:  [0, 1, 0, 1, 0, 1]
objective value: 38

It seems that the solution we got earlier was already the best one.

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 *

weights, k = data
n = len(weights)

# x[i] is the successor node of node i (in the cycle where i belongs)
x = VarArray(size=n, dom=range(n))

# y[i] is the cycle (index) where the node i belongs
y = VarArray(size=n, dom=range(n))

satisfy(
    AllDifferent(x),

    # ensuring correct cycles
    [y[i] == y[x[i]] for i in range(n)],

    # disabling infeasible arcs
    [x[i] != j for i in range(n) for j in range(n) if i == j or weights[i][j] < 0],

    # each cycle has k as maximum length
    BinPacking(y, sizes=1) <= k,

    # tag(symmetry-breaking)
    Precedence(y)
)

maximize(
    # maximizing the sum of arc weights of selected cycles
    Sum(weights[i][x[i]] for i in range(n))
)