# Combinatorial Optimization Problem using IBM ILOG CPLEX

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)

• It would help if you included your code for the other constraints. That way we could see what language you use (OPL or some programming language), and what variable names/data structures you use. May 25, 2021 at 6:24
• In future, please do not use images to show the problem and instead write it directly into the post (with MathJax), as it makes searching easier. I've done it for you this time as an example. May 26, 2021 at 6:48
• Hi TheSimpliFire, can you verify if I have written the constraints correctly? May 26, 2021 at 7:17
• Your code for constraint $(4)$ is not correct. The boolean string i!=j&&j!=k does not guarantee that i!=k, for example. May 26, 2021 at 8:01

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){
}
}


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.

• This is what I have written: int K =...;, int nbVertices =166;, range Vertices =1..nbVertices; int Reward[Vertices][Vertices] = ...; 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; forall(i in Vertices) sum(j in Vertices)Arc[j][i] == sum(j in Vertices)Arc[i][j]; 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 May 25, 2021 at 7:29

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).

• 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]; //Constraint (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[l][i] <= K-2; //Constraint 4 } May 25, 2021 at 7:48