Mixed Integer Programming - A complex example

You are given a row of 50 water buckets, each containing a specific concentration of salt. You also have a container that can hold exactly 3 buckets of water at a time. The task is to transfer a total of 30 buckets of water into an empty pool using only this container, which will require making 10 trips (since each trip can carry 3 buckets). 

The walking distance for each trip is determined by the distance between the nearest and farthest bucket out of the three buckets you choose to load into the container for that trip.

The goal is to minimize the total walking distances of the 10 trips, while the total concentration of salt in the pool remain below a specified threshold after transferring the 30 buckets into the pool.

Decision variables:

for each bucket, create an binary variable (0/1) yi to indicate if the bucket is selected

for the selection of 3 buckets to form a trip, this is tricky ... one way is to make every possible combination of 3 buckets a binary variable. 

z(i, j, k) = 0 / 1 indicates if bucket i, j and k are selected to form a trip.

The consequence is there will be many many z variables, depending on the available buckets. For a small number of buckets, e.g. 50, it creates 10000 or so variables, and a normal laptop still be able to handle it within a minute.

Objective function:

The total walking distance of all trips.

Given bucket i, j, k, the distance = max(i,j,k) - min(i,j,k) which can be expressed in a linear way by introducing again a lot of variables.

e.g. M >=x1, M>=x2, M>=x3

while minimize M would give you explicitly the maximum of x1,x2 and x3

Too much complexity and computation requirements.

Easy way is to pre-calculate the walking distance for every combination of i, j, k and store it as a dictionary or matrix. 

Then  sum(pre-calc walking distance(i,j,k) * z(i,j,k))  gives total walking distance, as z controls which triplet i,j,k is included.

Constraints:

Of course yi is 0 or 1

z(i,j,k) is also 0 or 1

The sum (yi) must == the selected buckets, i.e. 30

The sum (z(i,j,k)) must == the selected trips, i.e. 10

A tricky part is how to relate y to z? i.e. based on the selection y, what is the optimal trips that minimize the distance.

A triplet / trip z(i,j,k) is selected, only possible if yi, yj and yk are all selected, but it doesn't guarantee z(i,j,k) is selected. For example, the optimal solution may select g, h and i as one trip, while j, k, and h are selected for another trip. Therefore, z(i,j,k) == 1 means yi, yj, yk all = 1 but not vice versa.

To express that z(i,j,k) == 1 enforces yi, yj and yk to be 1.

 z(i,j,k)  <= yi

 z(i,j,k)  <= yj

 z(i,j,k)  <= yk

How to select the optimal i, j, k is up to the optimization algorithm, but the constraint here relates y to z.

However, notice that if a bucket i is selected by a triplet z(i,j,k)=1, then another variable z(i, k, m) must not be 1, because bucket i is not available for another triplet anymore. In another word, a bucket is used once and only once. Therefore, for every bucket yi, 

sum(i, ..., ...) + sum(..., i, ...) + sum(..., ..., i) <=1

Programmatically this can be done easily. refer to the code later.

The last one and an easy one is the concentration C.

sum( Ci * yi) < threshold

That's all. The code is below.

from ortools.linear_solver import pywraplp

import numpy as np

# Parameters

num_buckets = 52

selected_buckets = 30  #must be dividable by 3

trips = selected_buckets / 3 #each trip can take 3 buckets

C_max = 0.3  # Salt concentration threshold

# Dummy concentrations (assume some concentrations for the sake of example)

np.random.seed(42)  # Set a seed for reproducibility

C = np.random.rand(num_buckets)  # Concentration in each bucket

# Pre-calculate walk distances for each combination of 3 buckets, i.e. triplet (i, j, k)

# the distance between min and max buckets

# note the 3 loops is equivalent to get the unique combinations of i,j,k without repetition

walk_distance_dict = {}

for i in range(num_buckets):

    for j in range(i + 1, num_buckets):

        for k in range(j + 1, num_buckets):

            walk_distance_dict[(i, j, k)] = max(i, j, k) - min(i, j, k)

# Create the solver

# GLOP / CLP for LP, CBC for MIP, SCIP for MIP / Constraint programming,

# Xpress for LP, MIP, QP

# CP-SAT for constraint programming

solver = pywraplp.Solver.CreateSolver('SCIP')  # Use SCIP solver (available in OR-Tools)

# binary decision Variables to indicate if a bucket is selected or not

y = [solver.IntVar(0, 1, f'bucket_selected_{i}') for i in range(num_buckets)]  

# every possible combination (3 buckets i,j,k), a triplet

# z_ijk is binary to indicate if a triplet is selected

# the z_groups index the triplets by bucket id, as some constraints will apply later on the groups

z = {}

z_groups = {i:[] for i in range(num_buckets)}

for i in range(num_buckets):

    for j in range(i + 1, num_buckets):

        for k in range(j + 1, num_buckets):

            z[(i, j, k)] = solver.IntVar(0, 1, f'triplet_selected_{i}_{j}_{k}')

            z_groups[i].append(z[(i, j, k)]) #all triplets that contain bucket i

            z_groups[j].append(z[(i, j, k)]) #all triplets that contain bucket j

            z_groups[k].append(z[(i, j, k)]) #all triplets that contain bucket k

# Constraints

# total number of selected buckets

solver.Add(sum(y[i] for i in range(num_buckets)) == selected_buckets)

# Salt concentration must be below the threshold

solver.Add(sum(C[i] * y[i] for i in range(num_buckets)) <= C_max * selected_buckets)

# A triplet is selected, i.e. 1, only if all the 3 buckets are selected

# and enforce that z_ijk is 0 if any bucket is not selected

for i in range(num_buckets):

    for j in range(i + 1, num_buckets):

        for k in range(j + 1, num_buckets):

            # If a triplet (i, j, k) is selected, then all the three buckets must be selected

            solver.Add(z[(i, j, k)] <= y[i])

            solver.Add(z[(i, j, k)] <= y[j])

            solver.Add(z[(i, j, k)] <= y[k])

# The number of selected triplets equals to The total number of trips

solver.Add(sum(z[(i, j, k)] for i in range(num_buckets)

               for j in range(i + 1, num_buckets)

               for k in range(j + 1, num_buckets)) == trips)

#any bucket can be selected only once

# if a bucket is selected by one triplet, any other possible triplets can not select it again.

# so the sum of all triplets that contains i must be <= 1

for i, group in z_groups.items():

    solver.Add(sum(v for v in group) <= 1)    

# Objective function

# Minimize the total walk distance of the selected triplets

solver.Minimize(sum(walk_distance_dict[(i, j, k)] * z[(i, j, k)]

                   for i in range(num_buckets)

                   for j in range(i + 1, num_buckets)

                   for k in range(j + 1, num_buckets)))

# Solve the model

status = solver.Solve()

# Output the results

if status == pywraplp.Solver.OPTIMAL:

    print(f"Status: OPTIMAL")

    print("Selected Buckets:")

    for i in range(num_buckets):

        if y[i].solution_value() == 1:

            print(f"Bucket {i} is selected")

    print("\nTotal Walk Distance:")

    total_walk = solver.Objective().Value()

    print(f"Total Walk Distance: {total_walk}")

    print("\nWalk Distances per Selected Triplet:")

    for i in range(num_buckets):

        for j in range(i + 1, num_buckets):

            for k in range(j + 1, num_buckets):

                if z[(i, j, k)].solution_value() == 1:

                    print(f"Triplet ({i}, {j}, {k}): Walk Distance = {walk_distance_dict[(i, j, k)]}")

else:

    print("No optimal solution found.")

Status: OPTIMAL

Selected Buckets:

Bucket 4 is selected

Bucket 5 is selected

Bucket 6 is selected

...

Bucket 45 is selected

Bucket 46 is selected

Total Walk Distance: 21.0

Walk Distances per Selected Triplet:

Triplet (4, 5, 6): Walk Distance = 2

Triplet (8, 9, 10): Walk Distance = 2

Triplet (13, 14, 15): Walk Distance = 2

Triplet (17, 18, 19): Walk Distance = 2

Triplet (21, 22, 23): Walk Distance = 2

Triplet (26, 27, 28): Walk Distance = 2

Triplet (29, 31, 32): Walk Distance = 3

Triplet (36, 37, 38): Walk Distance = 2

Triplet (40, 41, 42): Walk Distance = 2

Triplet (44, 45, 46): Walk Distance = 2

ANother version that use itertools to replace the loops.

from ortools.linear_solver import pywraplp

import numpy as np

import itertools

# Parameters

num_buckets = 52

selected_buckets = 30  #must be dividable by 3

trips = selected_buckets / 3 #each trip can take 3 buckets

C_max = 0.3  # Salt concentration threshold

# Dummy concentrations (assume some concentrations for the sake of example)

np.random.seed(42)  # Set a seed for reproducibility

C = np.random.rand(num_buckets)  # Concentration in each bucket

# Pre-calculate walk distances for each combination of 3 buckets, i.e. triplet (i, j, k)

# the distance between min and max buckets

walk_distance_dict = {t: max(t) - min(t)  for t in itertools.combinations(range(num_buckets), 3)}

# Create the solver

# GLOP / CLP for LP, CBC for MIP, SCIP for MIP / Constraint programming,

# Xpress for LP, MIP, QP

# CP-SAT for constraint programming

solver = pywraplp.Solver.CreateSolver('SCIP')  # Use SCIP solver (available in OR-Tools)

# binary decision Variables to indicate if a bucket is selected or not

y = [solver.IntVar(0, 1, f'bucket_selected_{i}') for i in range(num_buckets)]  

# every possible combination (3 buckets i,j,k), a triplet

# z_ijk is binary to indicate if a triplet is selected

# the z_groups index the triplets by bucket id, as some constraints will apply later on the groups

z={t: solver.IntVar(0, 1, f'triplet_selected_{i}_{j}_{k}')  for t in itertools.combinations(range(num_buckets), 3)}

# put the triplets into groups by bucket number, this is for group-wise contraint later

z_groups = {i:[] for i in range(num_buckets)}

for (i,j,k) in z.keys():

    z_groups[i].append(z[(i, j, k)]) #all triplets that contain bucket i

    z_groups[j].append(z[(i, j, k)]) #all triplets that contain bucket j

    z_groups[k].append(z[(i, j, k)]) #all triplets that contain bucket k

# Constraints

# total number of selected buckets

solver.Add(sum(y[i] for i in range(num_buckets)) == selected_buckets)

# Salt concentration must be below the threshold

solver.Add(sum(C[i] * y[i] for i in range(num_buckets)) <= C_max * selected_buckets)

# A triplet is selected, i.e. 1, only if all the 3 buckets are selected

# and enforce that z_ijk is 0 if any bucket is not selected

for i,j,k in z.keys():

    # If a triplet (i, j, k) is selected, then all the three buckets must be selected

    solver.Add(z[(i, j, k)] <= y[i])

    solver.Add(z[(i, j, k)] <= y[j])

    solver.Add(z[(i, j, k)] <= y[k])

# The number of selected triplets equals to The total number of trips

solver.Add(sum(z.values()) == trips)

#any bucket can be selected only once

# if a bucket is selected by one triplet, any other possible triplets can not select it again.

# so the sum of all triplets that contains i must be <= 1

for i, group in z_groups.items():

    solver.Add(sum(group) <= 1)    

# Objective function

# Minimize the total walk distance of the selected triplets

solver.Minimize(sum(dis * trip for dis, trip in zip(walk_distance_dict.values(), z.values())))

# Solve the model

status = solver.Solve()

# Output the results

if status == pywraplp.Solver.OPTIMAL:

    print(f"Status: OPTIMAL")

    print("Selected Buckets:")

    for i in range(num_buckets):

        if y[i].solution_value() == 1:

            print(f"Bucket {i} is selected")

    total_walk = solver.Objective().Value()

    print(f"Total Walk Distance: {total_walk}")

    print("\nWalk Distances per Selected Triplet:")

    for i in range(num_buckets):

        for j in range(i + 1, num_buckets):

            for k in range(j + 1, num_buckets):

                if z[(i, j, k)].solution_value() == 1:

                    print(f"Triplet ({i}, {j}, {k}): Walk Distance = {walk_distance_dict[(i, j, k)]}")

else:

    print("No optimal solution found.")