PuLP - error when constraining variables that are mutually exclusive

Question

I'm using PuLP to try to solve an optimization problem, but keep getting "Status: Infeasible" when I try to set a constraint on the sum of the decision variables while they are also constrained by mutual exclusivity.

See my code below. I'm using binary variables y1, y2, and y3 to enforce mutual exclusivity on decision variables P1, P2, and P3, and y4 and y5 to enforce mutual exclusivity on P4 and P5. However, it seems like PuLP only mutually constrains y1, y2, y3, etc and not P1, P2, P3, etc.

from pulp import *

# Define the problem as a maximization problem
prob = LpProblem("AVB-Problem", LpMaximize)

# Define the decision variables
P1 = LpVariable("P1", 10, 20)
P2 = LpVariable("P2", 20, 30)
P3 = LpVariable("P3", 30, None)
P4 = LpVariable("P4", 15, 30)
P5 = LpVariable("P5", 30, None)

# Define the binary variables (used for mutual exclusivity rules)
y1 = LpVariable("y1", cat="Binary")
y2 = LpVariable("y2", cat="Binary")
y3 = LpVariable("y3", cat="Binary")
y4 = LpVariable("y4", cat="Binary")
y5 = LpVariable("y5", cat="Binary")

# Define the objective function
prob += 1*y1 + 1.5*y2 + 1.7*y3 + 2*y4 + 2.5*y5, "Z"

# Define the constraints
prob += P1 + P2 + P3 + P4 + P5 <= 50, "Total-Constraint"
prob += P1 <= y1*10000000000 + 10, "P1-Binary-Constraint"
prob += P2 <= y2*10000000000 + 20, "P2-Binary-Constraint"
prob += P3 <= y3*10000000000 + 30, "P3-Binary-Constraint"
prob += P4 <= y4*10000000000 + 15, "P4-Binary-Constraint"
prob += P5 <= y5*10000000000 + 30, "P5-Binary-Constraint"
prob += y1 + y2 + y3 == 1, "Mutual-Exclusivity-1"
prob += y4 + y5 == 1, "Mutual-Exclusivity-2"
# prob += y1 + y2 + y3 + y4 + y5 <= 2, 'Mutual-Exclusivity-3'

# Solve the problem
prob.solve(PULP_CBC_CMD(msg=1))

# Print the solution status
print(f"Status: {LpStatus[prob.status]}")

# Print the optimal solution and the optimal value of the objective function
for v in prob.variables():
    print(f"{v.name}: {v.varValue}")
print(f"Optimal Value of the Objective Function: {value(prob.objective)}")

# Get back the optimized P values
A1 = value(P1) * value(y1)
A2 = value(P2) * value(y2)
A3 = value(P3) * value(y3)
A4 = value(P4) * value(y4)
A5 = value(P5) * value(y5)

print(f"P1 Investment: {value(A1)}")
print(f"P2 Investment: {value(A2)}")
print(f"P3 Investment: {value(A3)}")
print(f"P4 Investment: {value(A4)}")
print(f"P5 Investment: {value(A5)}")

This returns the following error:

Status: Infeasible
P1: 10.0
P2: 20.0
P3: 30.0
P4: 15.0
P5: 30.0
y1: 0.0
y2: 1.0
y3: 0.0
y4: 1.0
y5: 0.0
Optimal Value of the Objective Function: 3.5
P1 Investment: 0.0
P2 Investment: 20.0
P3 Investment: 0.0
P4 Investment: 15.0
P5 Investment: 0.0

Note that PuLP respects mutual exclusivity on y1, y2, y3, etc but not on P1, P2, P3, etc. All of the P1, P2, P3, etc variables have non-zero value. I can "fix" this by changing the inequality P1 + P2 + P3 + P4 + P5 <= 50 to P1 + P2 + P3 + P4 + P5 <= 105, since 105 is the sum of the bounds on P1, P2, P3, P4 and P5, but I want to set a sum constraint of 50, not 105.

In essence, I want to do the following:

P1*y1 + P2*y2 + P3*y3 + P4*y4 + P5*y5 <= 50

However, this is non-linear and PuLP throws an error if I add that constraint. Is there a way to solve this with PuLP?

I guess it is the knapsack problem, but with the distinction that the P1, P2, P3, etc variables have lower bounds that cannot be 0. — wizardjoe, Apr 07 '23 at 14:52
What does exactly P1-3 mutually exclusive mean? They take value from range or zero, do they? Your objective depends on binaries. It's unclear why. What does your 'Mutual-Exclusivity-3' suppose to reflect? Please, write the original math formulation. It difficult to get it from your code. — Askold Ilvento, Apr 07 '23 at 15:01

score 2 · Accepted Answer · answered Apr 07 '23 at 15:04

I think you are looking for something like the below.

3 things to point out.

You should use constraints triggered by the selection variable to control the bounds... right now your hard-coded upper/lower bounds are fighting against your constraints. You only want to enforce those constraints if the selection occurs, so I re-wrote
Several of these issues are "Big M" constraints.... do a little research if that is confusing
I find it easier and more compact to index the variables instead of writing them out, so I modified.

Code:

from pulp import *

# Define the problem as a maximization problem
prob = LpProblem("AVB-Problem", LpMaximize)

M = 1000   # big M constraint value
lower_lims = [10, 20, 30, 15, 30]
upper_lims = [20, 30, M, 30, M]

selection_weights = [1, 1.5, 1.7, 2, 2.5]   # for objective

P = LpVariable.dicts('P', list(range(5)))  # NO bounds, we'll use constraints

select = LpVariable.dicts('y', list(range(5)), cat="Binary")


# Define the objective function
prob += sum(selection_weights[i] * select[i] for i in range(5)), "Z"

# Define the constraints
prob += sum(P[i] for i in range(5)) <= 50

# constrain the p values, based on selection variable
for i in range(5):
    prob += P[i] >= select[i] * lower_lims[i]
    prob += P[i] <= select[i] * upper_lims[i]

# mutual exclusivity constraints...
groups = [(0, 1, 2), (3, 4)]
for group in groups:
    prob += sum(select[i] for i in group) <= 1


# Solve the problem
prob.solve(PULP_CBC_CMD(msg=1))

# Print the solution status
print(f"Status: {LpStatus[prob.status]}")

# Print the optimal solution and the optimal value of the objective function
for v in prob.variables():
    print(f"{v.name}: {v.varValue}")
print(f"Optimal Value of the Objective Function: {value(prob.objective)}")

# Get back the optimized P values
for i in range(5):
    print(f'A{i}: {value(P[i]) * value(select[i])}')

Yields:

Status: Optimal
P_0: 0.0
P_1: 20.0
P_2: 0.0
P_3: 0.0
P_4: 30.0
y_0: 0.0
y_1: 1.0
y_2: 0.0
y_3: 0.0
y_4: 1.0
Optimal Value of the Objective Function: 4.0
A0: 0.0
A1: 20.0
A2: 0.0
A3: 0.0
A4: 30.0
[Finished in 97ms]

Great, this is awesome! I just realized that perhaps we don't even need to break out 2 different variables P and y, we could probably even just make P a boolean and optimize on P, then when constructing the constraint, multiply each P by their lower bound like ```(P1 lower bound) * P1 + (P2 lower bound) * P2 + ... <= 50``` — wizardjoe, Apr 07 '23 at 15:10
Well, it isn't clear... I assumed this was a piece of a larger problem and P might vary. If not, then as @Askold Ilvento points out, you might be taking the loooooong road to re-inventing a knapsack problem — AirSquid, Apr 07 '23 at 15:18
Would this still be valid if, instead of fixed selection weights, the selection weights are percentages of P? — wizardjoe, Apr 09 '23 at 04:11
Lay the math out on paper and see if you can make it linear. You are changing things around quite a bit here. Between the selections, proportions, and ranges, it can probably be made linear with some focus. — AirSquid, Apr 09 '23 at 19:04
I figured it out. I originally removed the y variable after the last time we worked on this, but added it back, and then set the objective function to P * weights, while adding y back to the upper/lower lim constraints. I guess the most confusing thing for me has been understanding the difference between bounds and constraints in PuLP and Pyomo. Thanks again for all of your help! — wizardjoe, Apr 09 '23 at 21:53

score 0 · Answer 2 · answered Apr 07 '23 at 14:06

0

You have defined all of your P variables with equal lower and upper bounds, which is odd/nonsensical. So they are fixed at that value.

That makes this constraint infeasible:

prob += P1 + P2 + P3 + P4 + P5 <= 50, "Total-Constraint"

As far as troubleshooting goes, try printing the model and reading through the math it produces, it is a good first T/S step and usually shines some light on things gone wrong:

print(prob)

AVB-Problem:
MAXIMIZE
1*y1 + 1.5*y2 + 1.7*y3 + 2*y4 + 2.5*y5 + 0.0
SUBJECT TO
Total_Constraint: P1 + P2 + P3 + P4 + P5 <= 50

P1_Binary_Constraint: P1 - 100000 y1 <= 10

P2_Binary_Constraint: P2 - 100000 y2 <= 20

P3_Binary_Constraint: P3 - 100000 y3 <= 30

P4_Binary_Constraint: P4 - 100000 y4 <= 15

P5_Binary_Constraint: P5 - 100000 y5 <= 30

Mutual_Exclusivity_1: y1 + y2 + y3 = 1

Mutual_Exclusivity_2: y4 + y5 = 1

VARIABLES
P1 = 10 Continuous
P2 = 20 Continuous
P3 = 30 Continuous
P4 = 15 Continuous
P5 = 30 Continuous
0 <= y1 <= 1 Integer
0 <= y2 <= 1 Integer
0 <= y3 <= 1 Integer
0 <= y4 <= 1 Integer
0 <= y5 <= 1 Integer

answered Apr 07 '23 at 14:06

AirSquid

10,214
2
7
31

I can make the upper bounds higher, or even None, but it won't make a difference. It still fails because the lower bounds force values on the variables, and adding them all up fails the 50 "Total-Constraint". I can set them to 0, but that leaves PuLP open to pick a number lower than my current lower bound, which is not what I want. For example, it could set P3 = 20, which is lower than 30, and it would set this all the time instead of setting P2 =20, since the reward in the objective function is higher. – wizardjoe Apr 07 '23 at 14:18
1

ahhh... ok perhaps I didn't read your question closely enough on the last part or it wasn't expressed well... the first portion seemed to inquire as to why this was infeasible. I see where you are going with this...stby, I'll modify the answer. It is doable. – AirSquid Apr 07 '23 at 14:30
Actually, I just noticed that if I set all the LB's to 0, PuLP actually leaves the P variables at 0! It seems like what's happening is that it only updates the y variables because the y variables are what are in the objective function. Perhaps the objective function needs to be rewritten with P1, P2, P3, etc, but no idea how to do that... – wizardjoe Apr 07 '23 at 14:36
Your current objective function is built upon the `y` variable, which is the selection variable, is that your intent? – AirSquid Apr 07 '23 at 14:44
It was, since it was easier to translate the rewards into coefficients that way. For example, I wanted to reward 1.5 if P2 took on some value greater than 20 (which, based on the binary constraint, would have set y2 = 1). But perhaps this is the wrong way to think about it. – wizardjoe Apr 07 '23 at 14:49

PuLP - error when constraining variables that are mutually exclusive

2 Answers2

Code:

Yields:

Linked