pythonlinear-programmingpulplinear-optimization

PuLP - error when constraining variables that are mutually exclusive

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?

Solution

  • I think you are looking for something like the below.

    3 things to point out.

    1. 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

    2. Several of these issues are "Big M" constraints.... do a little research if that is confusing

    3. 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]