Simple Integer Optimization Problem: docplex CP model works but equivalent PuLP+CBC model is infeasible?

Question

I have an integer optimization problem with one constraint per decision variable and no objective function. It can be coded and solved using docplex, however I am struggling to implement an equivalent model in PuLP.

Problem

Given a series with indices 0...n, find integer values of the series such that the value for each index i corresponds to the number of occurrences of i in the series.

For example, a solution for n = 4 is:

0	1	2	3	4
2	1	2	0	0

For index i=0, the value is 2, which corresponds to the series "2, 1, 2, 0, 0" containing two occurrences of the value 0. Index i=1 has a value of 1, corresponding to the single occurrence of 1 in the series. And so on.

Docplex implementation works!

from docplex.cp.model import CpoModel

# Data
n = 5
NUMBERS = range(n)

# Model instance
mdl = CpoModel(name='find series')

# Decision variables
s = mdl.integer_var_list(n, 0, n-1, "series")

# Constraint: Value should equal the # of occurrences of the index in the series

for i in NUMBERS:
    mdl.add(sum(s[j] == i for j in NUMBERS) == s[i])

# Solve
msol = mdl.solve()

PuLP attempts do NOT work!

Not sure if it matters, but I'm using the default CBC PuLP solver in both attempts below.

Attempt 1

In this attempt, the if statement of the constraint seems to always evaluate to True, such that the constraints become X[i] = 5, which is infeasible since the maximum value of the decision variables is n = 4.

from pulp import *

# Data
n = 4
NUMBERS = range(n+1)

# Model instance
model = LpProblem("FindSeries")

# Decision variables
X = LpVariable.dicts('X', NUMBERS, lowBound=0, upBound=n, cat=LpInteger)

# Constraint: Value should equal the # of occurrences of the index in the series
for i in NUMBERS:
    model += lpSum([1 for j in NUMBERS if X[j]==i]) - X[i] == 0

# Solve
model.solve()

The constraints look like:

_C1: - X_0  = -5
_C2: - X_1  = -5
_C3: - X_2  = -5
_C4: - X_3  = -5
_C5: - X_4  = -5

Attempt 2

When the above failed, I tried implementing the constraint another way, more similar to the docplex code.

# Constraint: Value should equal the # of occurrences of the index in the series
for i in NUMBERS:
    model += lpSum([X[j]==i for j in NUMBERS]) == X[i]

This time, I'm totally confused about what is going on. The model constraints look like this:

_C1: 0 X_0 + X_1 + X_2 + X_3 + X_4 = 0
_C2: X_0 + 0 X_1 + X_2 + X_3 + X_4 = 5
_C3: X_0 + X_1 + 0 X_2 + X_3 + X_4 = 10
_C4: X_0 + X_1 + X_2 + 0 X_3 + X_4 = 15
_C5: X_0 + X_1 + X_2 + X_3 + 0 X_4 = 20

Questions

• Is it possible to solve this using PuLP + default CBC solver?
  • If not, why not?
  • If so,
    • How can I code the model correctly in PuLP?
    • Why did my PuLP attempts fail as they did?

Answer 1

A simple equivalent MIP model can be built as follows.

Introduce a binary variable $y_{i,j}\in\{0,1\}$ such that $y_{i,j}=1\Leftrightarrow x_i=j$. The problem can now be formulated as a straight MIP model (just add a dummy objective):

$$ \begin{aligned} &x_i = \sum_j j\cdot y_{i,j}&& \forall i && \text{extract value}\\ &\sum_j y_{i,j} = 1 && \forall i && \text{exactly one value in each row of $y$}\\ &x_i = \sum_j y_{j,i} && \forall i && \text{count occurrences}\\ &y_{i,j} \in \{0,1\} \end{aligned} $$

If you want, you can substitute out $x_i$ and recover those values during reporting.

Answer 2

It looks like you are confusing the variables of your problem ($x_i$) with the parameters (the number of occurrences $s_i$). The way I understand your problem is as follows:

The mathematical expression of your constraint is $$ \sum_{j}x_j = s_i \quad \mbox{for all }i $$ where $s_i$ denotes the number of occurrences of value $i$ in the series.

With PuLP:

for i in s:
   model += lpSum(x[j] for j in x) == s[i]

---

EDIT

I misunderstood the question and am editing according the comments below. So the constraint is $$ \sum_{j}|\{j : x_j=i\}| = x_i \quad \mbox{for all }i $$

This is not a linear constraint and thus cannot be implemented as is in PuLp. I suggest to use Erwin's MIP below, which can be solved with PuLp.