Creating the model is too slow. Solving it is not the problem here.
I have looked at a few similar questions before, but they don't have the problem from the size that I am giving it.
The constraint that is very slow to add to the model looks something like this:
x out of {0, 1} (binary)
y_0 = 0
y_i = y_(i-1) + x_i; N > i > 0
in other words a variable vector x=[1, 0, 1, 1] would give me
y_1 = 1
y_2 = 1
y_3 = 2
It's actually more complex, but I simplified it here to this minimum reproducible example. The meaning is: I don't want the value y to go above a certain threshold for all i.
I used pulp. Here is an example python script:
from pulp import LpProblem, LpVariable, LpMaximize, LpBinary, PULP_CBC_CMD, lpSum
N = 5000
max_y_val = 5
model = LpProblem(sense=LpMaximize)
print("defining X")
x = [LpVariable(name=f"x_{i}", cat=LpBinary) for i in range(N)]
print("defining Y")
y = [x[0]]
for i in range(1, N):
print(i)
y.append(y[i-1] + x[i])
print("add y to model")
for i in range(N):
print(i)
model += y[i] <= max_y_val
print("objective")
model += lpSum(x) # ...
print("start optimizer")
status = model.solve(PULP_CBC_CMD(msg=False, timeLimit=10))
print("score:", model.objective.value())
print("x:", *[x[i].value() for i in range(N)])
print()
With N=1000 it still goes quite fast. At N=5000 it becomes very slow when adding y to the model.
I rewrote the same code in rust (using good_lp library) to see if the issue is with python, but the problem persists. So most likely the issue is the way I am writing my problem. Thus it is not a language specific question.
My suspicion is, that it might be using a matrix M * x = y. The width of the matrix with higher numbers of N would become N + 1. Probably something that looks like a triangular matrix that only has elements on the left_bottom corner.
matrix x y
1 0 0 0 0 0 0
1 1 0 0 0 1 1
1 1 1 0 0 * 0 = 1
1 1 1 1 0 1 2
1 1 1 1 1 1 3
Do you know how I can do this with LP solvers in general? I feel like I missing some knowledge here.
