Bakery Problem

Use the PuLP library https://pypi.org/project/PuLP/ to solve the following problems. Documentation to PuLP can be found here: https://coin-or.github.io/pulp/main/index.html

Bakery problem (20 points)

Consider a small bakery with a single oven. It needs to schedule baking of n pastries, each of them having three requirements:

time when the preparations are done and pastry is ready for baking time needed for baking, i.e., for how long should it remain in the oven deadline: time when the custommer comes to pick up the pastry At each moment, only one kind of pastry can be present in the oven. Use an ILP to find a shortest baking schedule. Schedule, in this context, is a set of starting times s1, ..., sn denoting when should each pastry be put into oven. Note: these times need not be integral. However, integral variables will be useful to enforce that the periods when two different kinds of pastries are in the oven do not overlap.

Let us denote $e_1$, ..., $e_n$ the ending times of baking of each of the pastries, i.e., $e_i = s_i + baking time of pastry_i$. We need to make sure that for each two pastries $i,j$, one of the following needs to be true: $e_i ≤ s_j$ or $e_j ≤ s_i$. Obviously, they cannot hold at the same time and it depends on the precedence between i and j which one is true. Since we do not know the precedence in advance, which of these constraints should we include in the LP?

Big-M method

This name usually refers to an alternative way how to start the simplex method without knowledge of the initial basic feasible solution. We did not cover this in the class and I do not go into details of this here either. But the other meaning of Big-M is a method for switching someof the constraints on/off depending on the value of some binary variable.

Imagine, we have variable x which should be bounded by $10$ if and only if some binary variable $z$ is set to zero. Also, assume that there is no reason to increase x beyound some large number $M$ (e.g., because we are minimizing over $x$, or we know that no feasible solution can have $x > M$ for some other reasons). Then, we can write $x ≤ 10 + Mz$: if $z$ is $0$, this switches the constraint ON. If $z=1$, this constraint evaluates to $x ≤ 10 + M$ which, by choice of $M$ is satisfied by any reasonable solution to our LP and this effectively switches the costraint OFF. Usually, due to possible numerical issues, it is recommended to use $M$ as small as possible. You can check the following blog for more discussion of big-M: https://orinanobworld.blogspot.com/2011/07/perils-of-big-m.html.

You may check that there is a suitable choice of $M$ in our problem and use this approach in your solution.

Input

Text file containing a single line for each kind of pastry consisting of four numbers (integers) separated by spaces:

ID PRE DLN BAK

ID denotes the numerical ID of the pastry, PRE the time since midnight since when the pastry is ready for baking, BAK is the time it needs to spend in the oven, and DLN is the deadline when the pastry needs to be surely finished. All times are in seconds.

Output

List of starting times of each pastry. Should look like this:

s1: 23.0
s2: 4.0
s3: 25.0
s4: 72.0
s5: 34.0

...

Bakery problem visualization (10 points)

Use library matplotlib to visualize your solution suitably. I leave to your creativity how to do it, but it should be clear what are the moments when oven needs to be open, what pastry goes out and what should be put in. There are many other things to visualize: expected arrivals of custommers and times when each pastry is ready, critical preparations (which pastry needs special care to be prepared on time, otherwise it would delay the whole schedule, etc). The main criterion for evaluation of this will be clarity and information it provides.

# Import packages. You can import additional packages, if you want
# You can change the way they are imported, e.g import pulp as pl or whatever
# But take care to adapt the solver configuration accordingly.
from pulp import *
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import time
import matplotlib.patches as mpatches
import matplotlib.lines as mlines

# Use the following solver
solver = COIN_CMD(threads=8)
# at home, you can try it with different solvers, e.g. GLPK, or with a different
# number of threads.
# WARNING: your code when run in vocareum should finish within 10 minutes!!!

def bakery():
    # Input file is called ./bakery.txt
    input_filename = './bakery.txt'

    # Use solver defined above as a parameter of .solve() method.
    # e.g., if your LpProblem is called prob, then run
    # prob.solve(solver) to solve it.

    # Load input data
    with open(input_filename, "r") as f:
        data = [list(map(int, line.strip().split())) for line in f]

    n = len(data)

    # Create table
    df = pd.DataFrame(data, columns = ['ID', 'PRE', 'DLN', 'BAK'])

    # Define the problem
    prob = LpProblem('Scheduling', LpMinimize)

    # We estimate the big-M value in such a way that is bigger than the deadlines,
    # but not big enough to create problems to our code
    M = 27001

    # Dependent variable
    y = LpVariable('y', lowBound=0, cat='Continuous')
    # Starting time variables
    s = LpVariable.dicts('s', [i for i in range(n)], cat = 'Continuous')
    # Auxiliary binary variables
    b = LpVariable.dicts('z', [(i,j) for i in range(n) for j in range(n)], cat=('Binary'))

    # Objective function
    prob += y

    # Add the constraints
    for i in range(n):
        prob += s[i] >= df['PRE'][i]
        prob += s[i] <= (df['DLN'][i] - df['BAK'][i])
        prob += y >= (s[i] + df['BAK'][i])
        for j in range(n):
            if i>j:
                prob += s[i] >= (s[j] + df["BAK"][j] - M * b[i, j])
                prob += s[j] >= (s[i] + df["BAK"][i] - M + M * b[i, j])

    # Solve the ILP
    prob.solve(PULP_CBC_CMD(msg=False))

    # extract solution
    retval = {}
    for i in range(n):
        retval[f's_{i}'] = s[i].value()

    # Create a dictionary storing start time - end time of each pastry
    startend = {}
    for i in range(n):
        startend[i] = (retval[f's_{i}'], retval[f's_{i}'] + df.loc[i,'BAK'])

    # Define the figure and axes for the plot
    fig, ax = plt.subplots(figsize=(15, 8))

    # Set title
    ax.set_title('Baking Schedule', fontsize = 20)

    # Sort the rows of the dataframe by the 'PRE' column
    df = df.sort_values(by='PRE')

    # Create a list of tuples with starting time and row index for each pastry
    pastry_order = [(startend[i][0], i) for i in range(n)]

    # I reverse the order because I want the top bars to be those corresponding to the first pastries
    pastry_order.sort(reverse=True)

    # Create a list of categorical values for the y-axis
    y_values = [str(i) for i in range(n)]

    # Reorder the y-axis categorical values based on the sorted pastry order
    y_values = [y_values[i[1]] for i in pastry_order]

    # Initialize a flag to keep track of the current label position
    label_above = True

    # Loop over each pastry in the sorted order
    for _, i in pastry_order:
        # Get the starting and ending times for the pastry
        start_time, end_time = startend[i]
        # Get the baking time for the pastry
        baking_time = df.loc[i, 'BAK']
        prep_time = df.loc[i, 'PRE'] - retval[f's_{i}']

        # Add a horizontal bar for the pastry
        ax.barh(str(i), baking_time, left=start_time, height=0.5, color='#FFA8A8')

        # Add a horizontal bar for the preparation time
        ax.barh(str(i), prep_time, left=start_time, height=0.1, color='yellow')

        # Add a horizontal bar for the time from the end of baking to the deadline
        deadline_bar_start = end_time
        deadline_bar_width = df.loc[i, 'DLN'] - deadline_bar_start
        ax.barh(str(i), deadline_bar_width, left=deadline_bar_start, height=0.1, color='#B7F0B1')

        # Create red border around 'critical' preparations (those in which the pastry has to be baked immediately)
        if prep_time == 0:
            # Add a red border to the baking time bar
            ax.barh(str(i), baking_time, left=start_time, height=0.5, edgecolor='red', lw=2, fill=False)
        else:
            # Add a regular red baking time bar
            ax.barh(str(i), baking_time, left=start_time, height=0.5, color='#FFA8A8')

        # Creat black star after the preparations in which pastries must be sold immediately after being baked
        if deadline_bar_width == 0:
            plt.scatter(start_time + baking_time + 250, f'{i}', marker='*',  color='red')

        # Add a label for the pastry
        ax.text(start_time + baking_time / 2, str(i), f"{i}", ha='center', va='center')

        # Add a vertical line from the end of the baking time bar to the x-axis
        ax.plot([end_time, end_time], [0, n], linestyle='dashed', color='gray', alpha=0.5)

        # Add a label for the end time of the baking time bar
        end_time_hours = int(end_time // 3600)
        end_time_minutes = int((end_time % 3600) // 60)
        end_time_seconds = int(end_time % 60)
        if end_time_seconds == 0:
            end_time_str = f"{end_time_hours:02d}:{end_time_minutes:02d}"
        else:
            end_time_str = f"{end_time_hours:02d}:{end_time_minutes:02d}:{end_time_seconds:02d}"

        # Choose the label position based on the current flag value
        if label_above:
            label_ypos = -0.4
            va = 'top'
        else:
            label_ypos = n + 0.1
            va = 'bottom'

        ax.text(end_time, label_ypos, end_time_str, ha='center', va=va, fontsize=8)

        # Switch the label position flag for the next iteration
        label_above = not label_above

    # Set the y-axis label
    ax.set_ylabel('Pastry ID')

    # Set the x-axis limits and label
    ax.set_xlim(0, 28200)
    ax.set_xlabel('Time (hour:minute)')

    # Define the tick locations and labels
    x_ticks = [0, 3600, 7200, 10800, 14400, 18000, 21600, 25200]
    x_ticklabels = ['0:00', '1:00', '2:00', '3:00', '4:00', '5:00', '6:00', '7:00']

    # Set the tick locations
    ax.set_xticks(x_ticks)

    # Set the tick labels
    ax.set_xticklabels(x_ticklabels)

    # Create a list of y-axis tick locations
    y_ticks = list(range(n))

    # Set the y-axis tick locations and labels
    ax.set_yticks(y_ticks)
    ax.set_yticklabels(y_values)

    # Remove the y-axis ticks
    ax.tick_params(axis='y', which='both', length=0)

    # Create the legend handles
    prep_patch = mpatches.Patch(color='yellow', label='Ready to Bake')
    bake_patch = mpatches.Patch(color='#FFA8A8', label='Baking')
    dead_patch = mpatches.Patch(color='#B7F0B1', label='Ready to sell')
    crit_patch = mpatches.Patch(facecolor='none', edgecolor='red', label='Critical Preparations')
    star_handle = mlines.Line2D([], [], color='red', marker='*', linestyle='None', markersize=10, label='Sell immediately')

    # Add the legend
    ax.legend(handles=[prep_patch, bake_patch, dead_patch, crit_patch, star_handle], loc='lower center', bbox_to_anchor=(0.5, -0.15), ncol=5)

    # Write visualization to the correct file:
    visualization_filename = './visualization.png'

    # Save the plot
    plt.savefig(visualization_filename, dpi=300, bbox_inches='tight')

    # retval should be a dictionary such that retval['s_i'] is the starting
    # time of pastry i
    print(value(prob.objective))
    return retval