# PuLP is ignoring constraints, and setting everything to 0 for minimization problem

I have a multi-objective optimization problem I am trying to solve with competing objectives. I am trying to model a network of industrial businesses which can share wastewater rather than sending it all to wastewater treatment.

I have set up one LpProblem() and want to alternate between minimizing the economic cost and minimizing the environmental impact. If I can achieve this I will use the optimal values obtained to create a weighted objective function to find a balance of the two. Where I am stuck is when I run the problem to minimize environmental impact, it finds the optimum to be 0 and ignores my water balance constraint which states that for each plant, all water going in must equal the demand of the plant and all water leaving must equal its available supply. I followed the method provided by @Reinderien and have done my best to produce a minimum reproducible example below.

from pulp import LpMinimize, LpProblem, LpStatus, LpStatusOptimal, lpSum, LpVariable, value, LpMaximize, lpDot

# input data
D =[100, 20, 80, 60, 150, 100]          # Demand in
S = [100, 20, 80, 60, 150, 100]         # Supply out (D[i] = S[i] no leakage)
Cin = [10, 40, 20, 30, 50, 35]          # Cin is the maximum contaminant concentration that can be accepted by plant j
Cout = [100, 250, 80, 200, 150, 130]    # Cout is the contaminant concentration of water exiting plant i
Cfresh = 1                              # concentration of contaminants in fresh water
Ctreat = [65, 115]                      # output concentration of treated wastewater from process k

plants = range(len(D))              # number of plants
treatments = range(len(Ctreat))     # number of treatment processes

# define the decision variables
a = LpVariable.matrix(name="TWW", lowBound=0, indices=(plants, plants))
e = LpVariable.matrix(name="UTWW", lowBound=0, indices=(plants, plants))
f = LpVariable.matrix(name="fresh", lowBound=0, indices=(plants))
w = LpVariable.matrix(name="waste", lowBound=0, indices=(plants))
# t[i][k] is the binary decision variable (1 if plant i uses treatment process k)
t = LpVariable.matrix(name="treatment", cat="Binary", indices=(plants, treatments))
# b[i][j][k] is the flow of treated water from i to j using process k
# this is supposed to be linearised so that PuLP can solve this problem
b = [[[a[i][j] if value(t[i][k]) == 1 else t[i][k] for k in treatments] for j in plants] for i in plants]

# defines the input flow into each plant
input_flow = [
fresh + lpSum(
a_row[j] + e_row[j]
for a_row, e_row in zip(a, e)
)
for j, fresh in zip(plants, f)
]

# defines the output flow from each plant
output_flow = [
waste + lpSum(a_row) + lpSum(e_row)
for a_row, e_row, waste in zip(a, e, w)
]


I have added my own input data and the additional variables I am working with. Specifically t[i][k] which is a binary variable that is set to 1 if plant i is using treatment process k. I use t to create a 3D array b[i][j][k] which takes the value of treated flow a[i][j] if the treatment process k is used. Otherwise it takes the value of t[i][k] which should be 0 if treatment k is not used. The idea behind this was to make sure the problem is still a linear problem.

Below, I calculate the contaminants as a product of the flow and concentration of each stream. These contaminants are used in the quality constraint and objective function. I also add the constraints and define the problem and objective function.

The quality constraint states that the contaminants entering a plant cannot exceed the maximum acceptable contaminant concentration Cin. It is taken from the case study I am following below.

Although there is a typo in this equation. The contaminants should be <= Cin[j]*D[j] so that the units are consistent.

# calculates the contaminants from the freshwater in, wastewater out and untreated water out
fresh_contam = lpDot(f, Cfresh)
waste_contam = lpDot(w, Cout)
untreated_contam = [
[
e[i][j] * Cout[i]
for j in plants
]
for i in plants
]

# calculate the sum of contaminants entering plant j from plant i
# turns a 3D array back into a 2D array
treated_contam = [
[
lpSum(
b[i][j][k] * Ctreat[k] for k in treatments
)
for j in plants
]
for i in plants
]

# quality constraint value that sums the contaminants entering plant j
input_contam = [
fresh + lpSum(
untreated_contam[i][j] + treated_contam[i][j]
for i in plants
)
for fresh, j in zip(fresh_contam, plants)
]

# define the problem and objective function
prob = LpProblem(name="treatment_flow)", sense=LpMinimize)
prob.setObjective(fresh_contam + waste_contam)

# ensures plants do not connect to themselves
for i in plants:
a[i][i] = 0
e[i][i] = 0

# the sum of all inputs to each destination plant must exactly equal its demand
for j, in_, in_target in zip(plants, input_flow, D):

# the sum of all outputs from each source plant must exactly equal its output flowrate
for i, out, out_target in zip(plants, output_flow, S):

# quality constraint ensures plant will only accept input water from sources as long as the max contaminant
# concentration Cin[j] is not exceeded
for j, contam, max_contam, demand in zip(plants, input_contam, Cin, D):



I'm unsure if the following constraint is necessary as my output didn't change with or without it

# additioanl constraint set regarding input_flow  and output_flow
for plant, in_, out in zip(plants, input_flow, output_flow):

# run problem and output results
prob.solve()
assert prob.status == LpStatusOptimal

print('Treated water flow rate (left to top):')
for arow in a:
print([value(x) for x in arow])
print('Untreated water flow rate (left to top):')
for erow in e:
print([value(x) for x in erow])
print("Freshwater input rate:", [value(x) for x in f])
print("Wastewater output rate:", [value(x) for x in w])
print("Total input flow:", [value(x) for x in input_flow])
print("Total output flow:", [value(x) for x in output_flow])
print("Binary variable if plant i uses treatment process k")
for trow in t:
print([value(x) for x in trow])


When I run the model, it is able to solve it but my main issue is that the freshwater in and wastewater out values are all 0 which is still wrong. These plants are in a system with input from the freshwater and output to the wastewater. Not every plant needs to be connected to the freshwater source or wastewater disposal site, but some need to be in order to bring water in and out of the system. Since there are no leaks in the system, the water in from freshwater should equal the water out through wastewater. Something is causing PuLP to disobey this water balance and I suspect it is a constraint I need to add but I am unsure as to what this constraint should be. It seems as if the plants are in a loop cycling their water together but ignoring the quality constraint as well.

• you can try to add constraints with $f_j$ and $w_i$ taking positive values and check to see if the solver is able to find a feasible solution. If it can't, it means that there are conflicting constraints somewhere. If it can, it means that your model is set up set to variables to $0$. Feb 1 at 7:59
• Well, your b is definitely a problem. That isn't whatsoever how binary variables work. I'll try to demonstrate a patched-up version. Feb 2 at 10:39

I don't think you should have binary constraints here at all. Infer "use" by a treatment-specific rate being non-zero. One obvious quirk is that, with no other constraints on available treatments, the treatment with the least contamination will always be used. Also, your treatments are all so low-contamination that no 'untreated' transfers are ever used; increase the treatment contamination to 200 to see the effect of mixed treated-untreated flows.

Reformulate to:

$$\min C_{untreat} \cdot w$$ $$s.t.$$ $$a_{iik} = 0 \quad \forall i, k$$ $$e_{ii} = 0 \quad \forall i$$ $$f_j + \sum_i \left(e_{ij} + \sum_k a_{ijk} \right) = d_j \quad\forall j$$ $$w_i + \sum_j \left(e_{ij} + \sum_k a_{ijk} \right) = s_i \quad\forall i$$ $$f_j C_{fresh} + \sum_i C_{treat} \cdot a_{ij} + C_{untreat} \cdot e_j \le C_{in(j)} \quad\forall j$$

## Demo code

This works fine:

import graphviz
import pulp

D = (100, 20, 80, 60, 150, 100)           # Demand into plant j
S = D                                     # Supply out of plant i (D[i] = S[i] no leakage)
Cin = (10, 40, 20, 30, 50, 35)            # contaminant concentration max accepted by plant j
Cuntreat = (100, 250, 80, 200, 150, 130)  # contaminant concentration of untreated water from plant i
Ctreat = (65, 115)                        # contaminant concentration of treated waste from process k
Cfresh = 1                                # contaminant concentration of fresh water

plants = range(len(D))              # indices of plants
treatments = range(len(Ctreat))     # indices of treatment processes

def make_vars() -> tuple[
list[list[list[pulp.LpVariable]]],  # a, treated water transfer
list[list[pulp.LpVariable]],        # e, untreated water transfer
list[pulp.LpVariable],              # f, fresh water in
list[pulp.LpVariable],              # w, waste water out
]:
a = pulp.LpVariable.matrix(
name='TWW',  # treated water transfer, (source, dest, treatment)
cat=pulp.LpContinuous, lowBound=0, indices=(plants, plants, treatments))
e = pulp.LpVariable.matrix(
name='UTWW',  # untreated water transfer, (source, dest)
cat=pulp.LpContinuous, lowBound=0, indices=(plants, plants))
f = pulp.LpVariable.matrix(
name='fresh',  # freshwater in, (dest)
cat=pulp.LpContinuous, lowBound=0, indices=plants)
w = pulp.LpVariable.matrix(
name='waste',  # wastewater out, (source)
cat=pulp.LpContinuous, lowBound=0, indices=plants)

# ensures plants do not connect to themselves
for i in plants:
a[i][i] = (0.,) * len(treatments)
e[i][i] = 0.

return a, e, f, w

def make_exprs(
a: list[list[list[pulp.LpVariable]]],
e: list[list[pulp.LpVariable]],
f: list[pulp.LpVariable],
w: list[pulp.LpVariable],
) -> tuple[
list[pulp.LpAffineExpression],  # total input
list[pulp.LpAffineExpression],  # total output
list[pulp.LpAffineExpression],  # total input contamination
pulp.LpAffineExpression,        # total waste contamination
]:
# input flow into each plant
input_flow = [
+ fresh
+ pulp.lpSum(
pulp.lpSum(a_slice[j]) + e_row[j]
for a_slice, e_row in zip(a, e)
)
for j, fresh in zip(plants, f)
]

# output flow from each plant
output_flow = [
+ pulp.lpSum(a_slice)
+ pulp.lpSum(e_row)
+ waste
for a_slice, e_row, waste in zip(a, e, w)
]

# total contaminants entering plant j
input_contam = [
+ fresh*Cfresh
+ pulp.lpSum(
pulp.lpDot(Ctreat, a[i][j])
for i in plants
)
+ pulp.lpDot(
Cuntreat,
[e[i][j] for i in plants],
)
for j, fresh in zip(plants, f)
]

# total wastewater contaminants exiting plant i
waste_contam = pulp.lpDot(Cuntreat, w)

return input_flow, output_flow, input_contam, waste_contam

def set_constraints(
prob: pulp.LpProblem,
input_flow: list[pulp.LpAffineExpression],
output_flow: list[pulp.LpAffineExpression],
input_contam: list[pulp.LpAffineExpression],
) -> None:
# the sum of all inputs to each destination plant must exactly equal its demand
for j, in_, in_target in zip(plants, input_flow, D):

# the sum of all outputs from each source plant must exactly equal its output flowrate
for i, out, out_target in zip(plants, output_flow, S):

# quality constraint ensures plant will only accept input water from sources as long as the max
# contaminant concentration Cin[j] is not exceeded
for j, contam, max_contam, demand in zip(plants, input_contam, Cin, D):

def dump_console(
a: list[list[list[pulp.LpVariable]]],
e: list[list[pulp.LpVariable]],
f: list[pulp.LpVariable],
w: list[pulp.LpVariable],
input_flow: list[pulp.LpAffineExpression],
output_flow: list[pulp.LpAffineExpression],
input_contam: list[pulp.LpAffineExpression],
) -> None:
print('Treated water rate (left to top):')
for treat in treatments:
print(f'  treatment {treat}:')
for arow in a:
print('    ', [pulp.value(x[treat]) for x in arow])
print()

print('Untreated water rate (left to top):')
for erow in e:
print([pulp.value(x) for x in erow])
print()

print('Fresh input rate:', [pulp.value(x) for x in f])
print('Fresh contamination:', [pulp.value(x)*Cfresh for x in f])
print('Total input rate:', [pulp.value(x) for x in input_flow])
print('Total input contamination:', [pulp.value(x) for x in input_contam])
print('Max input contamination:', [c*d for c, d in zip(Cin, D)])
print('Waste output rate:', [pulp.value(x) for x in w])
print('Waste contamination:', [pulp.value(x)*c for x, c in zip(w, Cuntreat)])
print('Total output rate:', [pulp.value(x) for x in output_flow])

def dump_graph(
a: list[list[list[pulp.LpVariable]]],
e: list[list[pulp.LpVariable]],
f: list[pulp.LpVariable],
w: list[pulp.LpVariable],
):
graph = graphviz.Digraph(
name='treatment_flow', format='svg',
graph_attr={
'rankdir': 'LR', 'overlap': 'false',
# 'nodesep': '1',
},
)

graph.node(name='fresh', label='Freshwater')
graph.node(name='waste', label='Wastewater')
for plant in plants:
graph.node(
name=str(plant),
shape='Mrecord',
label=
'{'
'{'
'<fresh_in> fresh|'
'<untreat_in> untreated|'
'<treat_in> treated'
'}|'
r'Plant \N|'
'{'
'<waste_out> waste|'
'<untreat_out> untreated|'
+ '|'.join(
f'<treat_{treatment}_out> treatment {treatment}'
for treatment in treatments
) +
'}'
'}'
)

for i, a_slice in enumerate(a):
for j, a_row in enumerate(a_slice):
for treatment, (contam, flow_var) in enumerate(zip(Ctreat, a_row)):
flow = pulp.value(flow_var)
if flow > 0:
graph.edge(
tail_name=f'{i}:treat_{treatment}_out',
label=f'{flow:.1f} ({contam*flow:.1f})',
)
for i, e_row in enumerate(e):
for j, flow_var in enumerate(e_row):
flow = pulp.value(flow_var)
if flow > 0:
graph.edge(
tail_name=f'{i}:untreat_out',
label=f'{flow:.1f}',
)
for j, flow_var in enumerate(f):
flow = pulp.value(flow_var)
if flow > 0:
graph.edge(
tail_name='fresh',
label=f'{flow:.1f} ({Cfresh*flow:.1f})',
)
for i, (flow_var, contam) in enumerate(zip(w, Cuntreat)):
flow = pulp.value(flow_var)
if flow > 0:
graph.edge(
tail_name=f'{i}:waste_out',
label=f'{flow:.1f} ({contam*flow:.1f})',
)

graph.view()

def main() -> None:
a, e, f, w = make_vars()
input_flow, output_flow, input_contam, waste_contam = make_exprs(a=a, e=e, f=f, w=w)
prob = pulp.LpProblem(name='treatment_flow', sense=pulp.LpMinimize)
prob.setObjective(pulp.lpSum(waste_contam))

set_constraints(
prob=prob,
input_flow=input_flow, output_flow=output_flow, input_contam=input_contam)

print(prob)
prob.solve()
assert prob.status == pulp.LpStatusOptimal

dump_console(
a=a, e=e, f=f, w=w, input_flow=input_flow,
output_flow=output_flow, input_contam=input_contam)
dump_graph(
a=a, e=e, f=f, w=w)

if __name__ == '__main__':
main()

treatment_flow:
MINIMIZE
100*waste_0 + 250*waste_1 + 80*waste_2 + 200*waste_3 + 150*waste_4 + 130*waste_5 + 0
SUBJECT TO
in_0: TWW_1_0_0 + TWW_1_0_1 + TWW_2_0_0 + TWW_2_0_1 + TWW_3_0_0 + TWW_3_0_1
+ TWW_4_0_0 + TWW_4_0_1 + TWW_5_0_0 + TWW_5_0_1 + UTWW_1_0 + UTWW_2_0
+ UTWW_3_0 + UTWW_4_0 + UTWW_5_0 + fresh_0 = 100
...

Optimal - objective value 27823.438
Optimal objective 27823.4375 - 23 iterations time 0.002
Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.01   (Wallclock seconds):       0.00

Treated water rate (left to top):
treatment 0:
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 20.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 60.0, 0.0]
[14.0625, 12.1875, 23.75, 27.1875, 0.0, 53.125]
[0.0, 0.0, 0.0, 0.0, 34.84375, 0.0]
treatment 1:
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

Untreated water rate (left to top):
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

Fresh input rate: [85.9375, 7.8125, 56.25, 32.8125, 35.15625, 46.875]
Fresh contamination: [85.9375, 7.8125, 56.25, 32.8125, 35.15625, 46.875]
Total input rate: [100.0, 20.0, 80.0, 60.0, 150.0, 100.0]
Total input contamination: [1000.0, 800.0, 1600.0, 1800.0, 7500.0, 3500.0]
Max input contamination: [1000, 800, 1600, 1800, 7500, 3500]
Waste output rate: [100.0, 0.0, 80.0, 0.0, 19.6875, 65.15625]
Waste contamination: [10000.0, 0.0, 6400.0, 0.0, 2953.125, 8470.3125]
Total output rate: [100.0, 20.0, 80.0, 60.0, 150.0, 100.0]


• I have taken your suggestions on and updated my post to provide more information on my problem. Apologies for not doing this originally. Feb 2 at 2:06
• @CalebO That's fine - see the second section in my answer and whether it answers your question. Feb 2 at 22:36
• Thank you this seems to solve my issue and is a very elegant solution. In my full problem there are costs associated with the different treatment processes so I'll see how that influences the choice in treatment processes. Feb 5 at 21:43