+ +
+

Scheduling as a graph coloring problem#

+
+

Preamble: Install Pyomo and a solver#

+

This cell selects and verifies a global SOLVER for the notebook. If run on Google Colab, the cell installs Pyomo and HiGHS.

+

Then, we set SOLVER to use the HiGHS solver via the Pyomo SolverFactory and verify that SOLVER is available.

+
+
+
import sys
+ 
+if 'google.colab' in sys.modules:
+    !pip install pyomo >/dev/null 2>/dev/null
+    !pip install highspy >/dev/null 2>/dev/null
+ 
+solver = 'appsi_highs'
+ 
+import pyomo.environ as pyo
+SOLVER = pyo.SolverFactory(solver)
+ 
+assert SOLVER.available(), f"Solver {solver} is not available."
+
+
+
+
+
+
+

Problem description#

+

In a small university, there is a single lecture hall that is suitable for exams. Consider the problem of choosing the time slots for each exam among a finite number of available time slots during the same week. Assume that the lecture hall is large enough to accommodate any number of exams in parallel. Since many courses are taken by the same students, exams for two courses that share at least one student cannot be scheduled in the same time slot. Every time slot the lecture hall is allocated for an exam has some indirect costs related to the hiring of temporary staff to invigilate the exams. The goal is thus to use as few time slots as possible while obeying the above restriction.

+
+
+

Mathematical formulation#

+

The scheduling problem above can be modeled as a graph coloring problem.

+

Consider an undirected graph \(G=(V,E)\) consisting of a set \(V\) of nodes and a collection \(E\subseteq V \times V\) of edges. Each node \(v \in V\) represents a course and two nodes are connected by an edge \(e \in E\) whenever at least one student takes the exams for both the corresponding courses.

+

Finding a conflict-free exam schedule is equivalent to finding a graph coloring of \(G\), i.e., an assignment of colors to the vertices of \(G\) such that no two adjacent vertices share the same color. In particular, in the formulation above we want to use the smallest number of colors possible.

+

We can formulate this problem using MILO as follows. Consider \(K\) possible colors (time slots) and define the following decision variables:

+
    +

  • \(w_{k} \in \{0, 1\}\), \(k = 1, \ldots, K\) is equal to 1 if the \(k\)-th color (time slot) is used;

    • +

  • \(x_{ik} \in \{0, 1\}\), \(i \in V\), \(k = 1,\dots,K\) is equal to 1 if the node (course) \(i\) is assigned the \(k\)-th color (time slot).

    • +
+

Using these variables, we can formulate the minimum graph coloring problem as:

+
+\[\begin{split} +\begin{align*} + \min \quad & \sum_{k=1}^K w_k \\ + \text{s.t.} \quad + & x_{ik} + x_{jk} \leq 1 && \forall \, (i, j) \in E, \, \forall \, k = 1,\dots,K \\ + & x_{ik} \leq w_k && \forall \, i \in V, \, \forall \, k = 1,\dots,K\\ + & \sum_{i \in V} x_{ik} = 1 && \forall \, k = 1,\dots,K\\ + & w_{k} \in \{ 0, 1 \} && \forall \, k = 1,\dots,K\\ + & x_{ik} \in \{ 0, 1 \} && \forall \, i \in V, \, \forall \, k = 1,\dots,K. +\end{align*} +\end{split}\]
+

The optimal solution of this problem is known as the chromatic number of the graph \(G\) and corresponds to the minimum number of time slots needed to schedule all exams.

+

In the rest of the notebook, for simplicity, we will use the graph terminology. The interested reader can easily interpret variables, constraints, and other considerations into the scheduling context.

+

The following cell defines the graph coloring problem in Pyomo. The function graph_coloring besides the graph G takes as optional argument the maximum number Kmax of colors to use in the model. Note that this might result in an infeasible model if Kmax is too small. The other optional argument catch_infeasibility is used precisely to diagnose the model infeasibility without raising errors.

+
+
+
import networkx as nx
+import matplotlib.pyplot as plt
+from time import perf_counter as pc
+
+
+def graph_coloring(G, Kmax=7, catch_infeasibility=False):
+    m = pyo.ConcreteModel("Graph Coloring")
+
+    m.colors = pyo.RangeSet(0, Kmax - 1)
+    m.nodes = pyo.Set(initialize=list(G.nodes))
+    m.edges = pyo.Set(initialize=list(G.edges), dimen=2)
+
+    m.x = pyo.Var(m.nodes, m.colors, domain=pyo.Binary)
+    m.w = pyo.Var(m.colors, domain=pyo.Binary)
+
+    @m.Constraint(m.edges, m.colors)
+    def edge_constraint(m, i, j, k):
+        return m.x[i, k] + m.x[j, k] <= 1
+
+    @m.Constraint(m.nodes)
+    def node_constraint(m, i):
+        return pyo.quicksum(m.x[i, k] for k in m.colors) == 1
+
+    @m.Constraint(m.nodes, m.colors)
+    def used_color_constraint(m, i, k):
+        return m.x[i, k] <= m.w[k]
+
+    @m.Objective(sense=pyo.minimize)
+    def number_used_colors(m):
+        return pyo.quicksum(m.w[k] for k in m.colors)
+
+    # solve the model with different options depending on whether we want to catch infeasibility
+    if catch_infeasibility:
+        # solve the model avoid raising an exception (by not loading the solution)
+        results = SOLVER.solve(m, load_solutions=False)
+        return results, m
+    else:
+        # solve the model and load the solution
+        SOLVER.solve(m, load_solutions=True)
+        return m
+
+
+def get_coloring(m):
+    coloring = {}
+    for i in m.nodes:
+        for k in m.colors:
+            if round(m.x[i, k]()) == 1:
+                coloring[i] = k
+    return coloring
+
+
+
+
+

Let us now solve this problem for a small example. We generate a random graph using the Erdos-Renyi \((n,p)\) model. Since we are interested in a connected graph, we keep generating random graphs changing the random seed until we obtained a connected one.

+
+
+
n = 14
+p = 3 / n
+
+seed = 1
+G = nx.gnp_random_graph(n, p, seed=seed)
+while not nx.is_connected(G):
+    seed += 1
+    G = nx.gnp_random_graph(n, p, seed=seed)
+
+pos = nx.spring_layout(G)
+nx.draw(G, pos=pos, with_labels=False)
+plt.show()
+
+
+
+
+../../_images/012c0ac7fc1ef68266b235563008020d2d0d8956166b1c2a0f4249e483bb05b6.png +
+
+

If we have a reasonable upper bound on the number of colors necessary, we can set it using the optional argument Kmax. Otherwise, the default for the model is Kmax=10.

+
+
+
Kmax = 4
+m = graph_coloring(G, Kmax)
+print(f"The chromatic number of G is: {int(m.number_used_colors())}")
+
+# Plot the graph with the coloring
+tab20 = plt.get_cmap("tab20")
+colors = [tab20(i) for i in range(0, 20, 4)]
+coloring = get_coloring(m)
+nx.draw(
+    G, pos=pos, with_labels=False, node_color=[colors[coloring[i]] for i in G.nodes]
+)
+
+
+
+
+
The chromatic number of G is: 3
+
+
+../../_images/8b81c2a184c076c22ffa803d1d8ea294a5a37b812824dd1cffab7845b1aed61c.png +
+
+

Alternatively, especially if we do not have a good guess for the upper bound of colors which are necessary, instead of choosing a very large Kmax which would make the MILO much larger and slower to solve, we can take an alternative approach.

+

Let us start by creating a much larger connected graph with the same strategy as before.

+
+
+
n = 40
+p = 10 / n
+seed = 1
+G = nx.gnp_random_graph(n, p, seed=seed)
+while not nx.is_connected(G):
+    seed += 1
+    G = nx.gnp_random_graph(n, p, seed=seed)
+
+pos = nx.spring_layout(G)
+nx.draw(G, pos=pos, with_labels=False, node_size=100, width=0.25)
+plt.show()
+
+
+
+
+../../_images/eea2e649c26e8415e05af2dc7aea168b29d2ba3be8955fc81e9bd4beba2661f2.png +
+
+
+
+
start_time = pc()
+Kmax = 2
+while True:
+    results, m = graph_coloring(G, Kmax, catch_infeasibility=True)
+    if results.solver.termination_condition == pyo.TerminationCondition.optimal:
+        m = graph_coloring(G, Kmax)
+        print(f"The chromatic number of G is: {Kmax}")
+        print(f"Time elapsed: {pc() - start_time:.2f} seconds")
+        break
+    else:
+        print(f"{Kmax} colors are not enough")
+        Kmax += 1
+
+# Plot the graph with the coloring
+tab20 = plt.get_cmap("tab20")
+colors = [tab20(i) for i in range(0, 20, 4)]
+coloring = get_coloring(m)
+nx.draw(
+    G,
+    pos=pos,
+    with_labels=False,
+    node_size=100,
+    width=0.25,
+    node_color=[colors[coloring[i]] for i in G.nodes],
+)
+
+
+
+
+
2 colors are not enough
+3 colors are not enough
+4 colors are not enough
+The chromatic number of G is: 5
+Time elapsed: 6.99 seconds
+
+
+../../_images/d122e190c38f57c64f87900ee69e5f9ddc764fe94eb3cfc6f5154678d1b51bfa.png +
+
+

This approach is much faster than solving the same MILO with a large Kmax, say 10, which results in a much larger model, both in terms of number of variables and constraints.

+
+
+
start_time = pc()
+m = graph_coloring(G, Kmax=10)
+print(f"The chromatic number of G is: {int(m.number_used_colors())}")
+print(f"Time elapsed: {pc() - start_time:.2f} seconds")
+
+
+
+
+
The chromatic number of G is: 5
+Time elapsed: 75.88 seconds
+
+
+
+
+

The networkx library has a greedy graph coloring algorithm which we can use to find a good upper bound for the number of colors needed. Note, however, that this is not guaranteed to be the minimum number of colors needed.

+
+
+
start_time = pc()
+nxcoloring = nx.algorithms.coloring.greedy_color(G)
+Kmax = len(set(nxcoloring.values()))
+print(
+    f"Upper bound found by networkx: Kmax={Kmax}. This will be used as the parameter of the optimization model."
+)
+m = graph_coloring(G, Kmax)
+print(f"The chromatic number of G is: {int(m.number_used_colors())}")
+print(f"Time elapsed: {pc() - start_time:.2f} seconds")
+
+
+
+
+
Upper bound found by networkx: Kmax=6. This will be used as the parameter of the optimization model.
+The chromatic number of G is: 5
+Time elapsed: 7.24 seconds
+
+
+
+
+
+
+ + + + +