Combinatorial Optimization Problem using IBM ILOG CPLEX

Question

Consider the combinatorial optimization problem described as below.

Let $D=(V,A)$ be a directed graph with $V$ the set of vertices and $A$ the set of arcs, i.e., $A=\{(i,j)\mid i,j\in V\}$. On each arc $a\in A$, there is a reward $r_a$ associated with it. The optimization problem is to form one or more cycles using the arcs, with each arc used in no more than one cycle. The cardinality of the cycle(s) cannot be more than a given number $K$. The objective is to maximize the sum of rewards for the arcs used in the cycles.

The major differences between this problem and the well-known Asymmetric Travelling Salesman Problem (ATSP) are as follows.

1. A solution that contains more than one "subtour" is considered feasible, however, the size of each subtour must be no more than $K$ vertices.

2. A feasible solution does not necessarily have to involve all vertices in $V$.

3. The objective is to maximize total reward of arcs used instead of to minimize the total cost of arcs used.

Let:

• $x_{i,j}\in\{0,1\}$ be a binary variable with $x_{i,j}=1$ indicating arc $(i,j)\in A$ is used in the solution and $0$ otherwise;

• $\Bbb P_K$ be the set of all simple paths with exactly $K$ vertices; and

• $\tau=(i_1,\ldots,i_K)\in\Bbb P_K$ be an arbitrary simple path with exactly $K$ vertices, i.e., the arcs in the path are $(i_1,i_2),(i_2,i_3),\ldots,(i_{K-1},i_K)$.

An integer linear programming (ILP) model is given as below. \begin{align}\max&\quad\sum_{(i,j)\in A}r_{ij}x_{ij}\tag1\\\text{s.t.}&\quad\sum_{(i,j)\in A}x_{ij}\le1,&\quad\forall i\in V\tag2\\&\quad\sum_{(j,i)\in A}x_{ji}=\sum_{(i,j)\in A}x_{ij},&\quad\forall i\in V\tag3\\&\quad x_{i_1,i_2}+x_{i_2,i_3}+\cdots+x_{i_{K-1},i_K}-x_{i_K,i_1}\le K-2,&\quad\forall(i_1,\ldots,i_K)\in\Bbb P_K\tag4\\&\quad x_{ij}\in\{0,1\},&\quad\forall(i,j)\in A\tag5.\end{align}

Can someone help with writing up constraint $(4)$ in CPLEX? I don't how to write this one constraint, can you do it for $K=4$?

My code is currently as follows.

int K =...;
int nbVertices =166;
range Vertices =1..nbVertices;
int Reward[Vertices][Vertices] = ...;
dexpr int totalReward = sum (i in Vertices, j in Vertices)Reward[I] 
[j]*Arc[i][j];
maximize totalReward;


subject to { 
forall(i in Vertices)
sum(j in Vertices)Arc[i][j] <=1; //constraint (2)
forall(i in Vertices) 
sum(j in Vertices)Arc[j][i] == sum(j in Vertices)Arc[i][j];//(3)
forall(i,j,k,l in Vertices: i!=j && j!=k && k!=l && l!=i) Arc[i][j]+
 Arc[j][k] + Arc[k][l] - Arc[k][i] <= K-2; //constraint (4)

Answer 1

In pseudocode, you'd have to do the following:

1. Create a binary $x_{ij}$ variable for every arc $(i,j)\in A$. Use the boolVar(.) method in the IloCplexModeler class. Personally I find it convenient to create a Map that maps an edge to its variable to quickly lookup the variable that corresponds to a given arc. If you use Java, you would have to define an Arc class, if you use Python, you could use a tuple.
2. Generate the set $P_k$ of all simple paths of length $k$. Each path is a sequence of arcs.
3. For every path in $P_k$, create a constraint of type 4. Iterate over the arcs in the path, and use your map to lookup the corresponding variables.

Step 3 would look something like this:

for(Path path : paths){
  IloLinearNumExpr expr = cplex.linearNumExpr();
  for(Arc arc : path){
    expr.addTerm(1.0, varMap.get(arc));
  }
  cplex.addLe(expr, K-2, "cycleConstraint");
}

The above, simplified, code snippet does not yet include the return arc which has a negative coefficient, but I hope you get the idea.

If you happen to use java, I would recommend to model the graph using the graph library JGraphT. Then you can easily iterate the desired paths.

Answer 2

Here (4) refers to a set of constraints, written for all possible paths with exactly $K$ vertices. The formulation is therefore non-compact which means the no. of these constraints grow exponentially as the problem size (no. of nodes) increases. In this case, there would be $\frac1n\cdot {}^nP_r$ constraints, where $n = |V|$, $r = K$ and $P$ represents permutations. Hence, it is difficult to add all such constraints in one go and also the LP model becomes unnecessarily bulky for CPLEX branch and bound.

These constraints are therefore often relaxed initially and added as and when relevant violations are found using a cutting plane approach (prior to branch and bound) and then branch and cut approach (during branch and bound). In both cases you need to implement a separation routine that will check the current LP relaxed solution for possible violations of a constraint from (4) and subsequently add that constraint to the model if a violation is found. When using the branch and cut approach in CPLEX, you will need cut callback (https://www.ibm.com/docs/en/icos/20.1.0?topic=legacy-cut-callback).