PuLP Conditional Constraint Sum (Based on Variable)

1
votes

I am trying to optimize a problem where the decision variable (X) is a NxN binary matrix. I am finding trouble making one of the constraints of the problem work.

The first constraint implies that the sum of X in each row must be == 1. (Covered)

The second constraint (the one I cannot make to work) requires that for those columns where the diagonal is == 1, the sum of X must be >= 2. I have generated the following constraint in PuLP:

for j in W:
    prob +=  sum(X[i][j] for i in W if X[j][j] >= 1) >= 2

PuLP indicates the solution status is 'Infeasible'. ¿What am I getting wrong? ¿This constraint cannot be implemented in PuLP?

Essentially, an example solution matrix that covers the previous requirements would be:

[0,0,0,1,0]
[0,0,0,1,0]
[0,0,0,1,0]
[0,0,0,0,1]
[0,0,0,0,1]
1
pythonoptimizationconstraintslinear-programmingpulp

1 Answers

1
votes

This is a linear-programming framework and therefore only linear functions can be used. Arbitrary python-expressions like:

  • if
  • abs
  • min

are not passed to the solver, but are evaluated a-priori (the modelling-framework / solver do not see this) usually resulting in garbage.

These kind of things must be linearized. Some frameworks can provide some of those automatically, but that is a very tough task in general (as most linearizations are impossible without model-specific assumptions!).

(cvxpy for example supports re-formulation of abs and min while being compatible with some rule-set)

Here you will need to do this yourself.

If i understood the task correctly, it looks like:

  • binary-matrix
  • every row sums to 1
  • every col sum is unconstrained, except for a column, where there is a 1 at the diagonal position -> then this column-sum has a lower-bound of 2

Example solution:

D!    D!

1  0  0  0   sum = 1 
1  0  0  0   sum = 1
0  0  1  0   sum = 1
0  0  1  0   sum = 1

sums
2  0  2  0

We can linearize the expression as:

for column in columns:
    (1-diag_elem_of_column) * bigM + sum(column) >= 2

    <-> 

    (1-diag_elem_of_column) * 2 + sum(column) >= 2

This is a classic big-M formulation, where there is some binary indicator variable (diagonal-element) activating/deactivating some linear expression using some large constant big-M. This large constant makes those approaches dependent on assumptions, but in your case this large constant does not need to be very large (and tightening is important for the integrality-gap -> better for the solver).

The idea is:

  • Enforce sum(column) >= 2 except when there is no diagonal-element hit
  • When there is a diagonal-element hit -> (1-diag_elem_of_column) * 2 = 0 -> the sum is constrained to be >= 2 (as we need to get 2 non-zeros in the remaining variables)
  • When there is no diagonal-element hit -> (1-diag_elem_of_column) * 2 = 2 >= 2 (sum of remaining elements is free -> no restriction as we already satisfied the minimum of 2)

Code-wise this could look like (untested):

# assumption: square-matrix
# assumption: X[row][col] storage-order

N = 5

for diag_element_index, col in enumerate(N):
  prob += (1 - X[diag_element_index][col]) * 2 + sum([X[row][col] for row in range(N)]) >= 2