# Extracting a sequence as a decision variable from a matrix (PuLP)

I am solving a sequencing problem x[i][j], i.e. a binary matrix indicating that product i is followed by j. Each product is to be exactly scheduled once each, and i must be before j if x[i][j] = 1, which is enforced by auxiliary variable p[i]. This decision variable shows the index of the i-th element or row in the sequence. For the example below, p[i], would be [0 2 1 4 5], read as entity/row 1 is in the 0th index in the sequence, entity or row 2 is at the 2nd index in the sequence, entity or row 3 at the 1st index of the sequence and so on. I wrote this optimization in PuLP / Python.

So far, this problem works, but for some reason, I want to combine/integrate it into another problem. For this new problem, I need the exact sequence as in x[i][j] as a decision variable.

So, using constraints, I want to establish the produced sequence as another decision variable. Going by the example below (codeblock), this would be an array like [1 3 2 4 5] when counting from 1. This decision variable I want to call o[i]. I could not figure out how to establish o[i] as a decision variable using only linear constraints on the given variables.

example is below in the code block

I have excluded my initial solution and costs function for privacy reasons!

# the decision variable is a succesion matrix of i by j, i.e. whether i is followed by j
x = LpVariable.dicts("transition", (range(N), range(N)), 0, 1, cat='Binary')

# example x[i][j], which is N by N
# [0 0 1 0 0]
# [0 0 0 1 0]
# [0 1 0 0 0]
# [0 0 0 0 1]
# [0 0 0 0 0]
#
# read as (row) i is succeeded by (column) j
# so, since the 1st column has no predecessor  (column sum 0) , it is the initial product
# thus, starting in row 1, it can be observed that the sucessor is column 3, so product 3
# so 1 -> 3, now do the same for row 3, 1 -> 3 -> 2
# the order in this example is 1 -> 3 -> 2 -> 4 -> 5

# auxiliary variable to enforce a single starting point
#u = LpVariable.dicts("position", range(N), 1, N)

# axuiliary variable to maintain position integrity
# p[i] variables (integer: position of entity i in the loop)
p = LpVariable.dicts("position", range(N), lowBound=0, upBound=N-1, cat='Integer')

prob.setObjective(costs)

# constraint 1: each product has exactly one sucessor
# we do not enforce this hard, since the last scheduled has no successor
for i in range(N):
prob += lpSum(x[i][j] for j in range(N)) <= 1

# constraint 2: each product is preceded by exactly one other
# we do not enforce this hard, since the first one does not have a predecessor
for j in range(N):
prob += lpSum(x[i][j] for i in range(N)) <= 1

# constraint 3: the total number of movements is at max N-1
# this means every product is selected once, except the start product
prob += lpSum(x[i][j] for i in range(N) for j in range(N)) == N - 1

# constraint 4: avoid circular dependencies, i.e. 1 -> 2 -> 1 is not allowed
# this means a sequence cannot loop back to a previously visited point
for i in range(N):
for j in range(N):
prob += x[i][j] + x[j][i] <= 1

# constraint 7: the start product has no predecessor
prob+= lpSum(x[i][zero_j] for i in range(N)) == 0

# constraint 8: the start product must have a successor
prob+= lpSum(x[zero_j][j] for j in range(N)) == 1

# constraint 9: the upper bound for the costs is the baseline value
# this baseline comes from warm_start
prob += costs <= baseline

# constraint 10: positional constraint
# if i is in the sequence, its position must be smaller than j if j is after i
for i in range(N):
for j in range(N):
prob += p[i] + 1 <= p[j] + M * (1-x[i][j])

# constraint 11
# the sum of assigned positions should be equal to the sum of range N
# i.e. each position is present once
prob += lpSum(p[i] for i in range(N)) == sum_n
• Search for MTZ constraints. Dec 28, 2023 at 15:23

Let's assume that $$i,j\in \lbrace 1,\dots,N\rbrace$$ and $$o_i$$ is an integer variable with domain $$\lbrace 1,\dots,N\rbrace$$ for each $$i.$$ You can add the constraints $$o_i + x_{i,j} - N(1-x_{i,j}) \le o_j \le o_i + x_{i,j} + N(1-x_{i,j})$$ for every index pair $$i\neq j.$$ If $$x_{i,j}=1,$$ this forces $$o_j = o_i + 1.$$ If $$x_{i,j}=0,$$ the corresponding constraints are nonbinding.

• Thanks, I implemented is as follows. for i in range(N): for j in range(N): prob += o[i] + x[i][j] - N*(1-x[i][j] <= o[j] <= o[i] + x[i][j] + N*(1-x[i][j]) However, it produced the same results as p[i], i.e. did not generate the exact produced sequence in x[i][j]. In essence, what am I seeking, is a decision variable which orders elements over range(N), in this case integers, according to values in p[i]. Dec 28, 2023 at 18:17

I was able to solve it as follows

Created $$z_{ij}$$, shaped similarly to $$x_{ij}$$

# zij is auxaillary to determine oi
z = LpVariable.dicts("auxiliary", (range(N), range(N)), cat='Integer')

Then, fill with 1's such that the column index multiplied with the '1' equals the position $$p$$ of the $$i$$-th element. Also enforce binary

for i in range(N):
prob += lpSum(j*z[i][j] for j in range(N)) == p[i]

for j in range(N):
prob += z[i][j] <= 1
prob += z[i][j] >= 0

In addition, make sure no values are redundantly placed

for i in range(N):
prob += lpSum(z[i][j] for j in range(N)) == 1

for j in range(N):
prob += lpSum(z[i][j] for i in range(N)) == 1

Then, the exact order can be deduced by taking the index of the row and multiply it with the 1.

for j in range(N):
prob += o[j] == lpSum(i*z[i][j] for i in range(N))