Solving a Resource Planning Problem with Mathematical Programming and Column Generation

Column generation for solving the graph coloring problem with Python and Gurobipy

I used to work in the oil and gas industry and one of the toughest problems we faced was resource planning and forecasting. My team and I constantly struggled with this aspect, often improvising and constantly putting out fires 🚒 . Just merely reacting to the situation, with the hope of finding a solution. But as Mark Wahlberg says in his movies "Hope ain’t a tactic".

"But improvisation will only bring you as far as the next crisis, and is never a substitute for thinking several steps ahead and planning to the end."

Robert Greene, The 48 Laws of Power

I always felt there had to be a better way to solve this problem, but I didn’t know it at the time 🤔. In fact acquiring the skills to solve such problems was one of my main motivations for continuing my studies. This article explores a specific application of mathematical optimization and how it could have made my job substantially easier back then. I guess that I am making this article as a gift 🎁 to my past-self.

I worked for a service company in the drilling and measurements division. The company provided tools and personnel to drill oil wells for various clients. Once a wellbore was drilled, the drilling tools were pulled out, sent to a workshop for maintenance and repairs, and then reused for another job. The start date for jobs was mainly determined by the client, so the crew and tools had to be onsite by that date. Based on past data, we could predict a rough estimate for the duration of different jobs. In the country where I worked, the oil and gas service industry was highly competitive. When a client gave us a job, we jumped into it without hesitation, accepting all jobs and figuring out the tool logistics later. Fortunately, most jobs were assigned at least a few weeks in advance. For example, a simplified month’s schedule for planned jobs might look like Figure 1 below (In reality we were dealing not just with 10 jobs per month but more like ~ 50).

Every single month we faced the same problem and always struggled to find the answer. Today, I am not surprised by why it was so difficult to solve; it turns out the problem is a hard combinatorial optimization problem. The problem is as follows: Given a set of jobs with their planned start and duration (thus end date), what is the minimum number of toolsets needed to cover them all? A toolset cannot move directly from one job to the next without a minimum service time, which, in my case, was three days. The goal is to find a lower bound for the number of tools needed for a given month so the maintenance team can plan ahead. Please note that a toolset can only be in a single job at the time.

Now it turns out that this problem is known in the literature as the "minimum vertex coloring problem" more popularly known as the "graph coloring problem" and at it has its origins on an older problem known as the map coloring problem [1].

The graph coloring problem has its roots in the famous four color theorem, first conjectured in 1852 by Francis Guthrie [2] , which posits that any map can be colored using only four colors such that no two adjacent regions share the same color.

More formally, the graph coloring problem is defined as follows: Given an undirected graph G(V,E) with nodes V and edges E, the objective is to assign colors to each node (vertex) such that no two adjacent nodes share the same color, and the total number of colors used is minimized. For example, the graph below in Figure 3, with six nodes, can be colored using only two colors.

At this point, you might be wondering what is this guy talking about❓ what is the relationship between coloring a map or a graph and solving the resource planning problem described above 🤔. It might not be evident at first glance, but both problems are essentially the same. My specific planning problem can be represented as a graph: each job is a node, and all jobs that overlap in time or start within three days of the finishing of another job are connected by an edge. Therefore, finding the minimum number of toolsets required is akin to finding the minimum number of colors needed to color this graph [3].

Moreover, there are numerous problems in various fields and industries that are equivalent or nearly equivalent to my planning problem. For instance:

• Exam planning: Given a set of courses such as Mathematics, Latin, Biology, etc., and students enrolled in some of these courses, the goal is to find the minimum number of timeslots needed for scheduling the final exams. This ensures that no student has conflicting exams (i.e., more than one exam in the same timeslot).
• Nurse and operating room planning: Given an amount of operations that are required to be scheduled in a given day, what is the minimum number of operating rooms or nurses that are required.
• Number of airplanes planning: Given a number of flights ✈️ within a day, each with a scheduled start time and duration, what is the minimum number of airplanes needed? considering that there is a required minimum service downtime after each flight.
• Number of train platform planning: Given a number of trains 🚆 arriving at a station at specific time slots, what is the minimum number of platforms required to handle the trains?
• Frequency assignment in radio stations: Given a set of radio stations 📻 within a certain distance of each other, what is the minimum number of radio frequencies needed to avoid signal conflicts between them?
• Computer register allocation problem: This problem involves assigning a limited number of CPU registers to a program’s variables during execution [4]. This problem is linked to graph coloring because it can be modeled as a graph where each node represents a variable and an edge connects two nodes if the corresponding variables are needed at the same time. Solving this involves coloring the graph with the minimum number of colors (registers) so that no two adjacent nodes (simultaneously needed variables) share the same color (register).
• Sudoku: Yes you are reading this correctly. Solving Sudoku is essentially solving a graph coloring problem because each Sudoku puzzle can be represented as a graph where each cell is a vertex, and edges connect vertices that share a row, column, or 3×3 subgrid. The goal is to assign a color (number 1–9) to each vertex (cell) such that no two adjacent vertices (cells sharing a row, column, or sub-grid) have the same color, analogous to ensuring no repeated numbers in those regions. For this specific setting we are not trying to minimize the number of colors used, what we want is to find a feasible solution with 9 colors.

Essentially, knowing how to solve the graph coloring problem enables you to solve many equivalent problems. Once you recognize one optimization problem, you begin to see them everywhere. It is like looking at the world with a brand new pair of glasses.

Before celebrating too early, I want to remind you that this is a very difficult combinatorial optimization problem. In fact, it is so challenging that it falls into the computational complexity category of NP-Hard. In this article, I will start by explaining the standard integer linear programming formulation for the graph coloring problem. However, due to its combinatorial complexity, this formulation might be prohibitive for solving large instances, even with the best available solvers. Therefore, I will introduce a powerful technique called column generation and how it can be applied to this problem. We will tackle this problem using Python 🐍 , Google Colab as our IDE, and Gurobi as our solver of choice. If you have never used Gurobi before, especially from Google Colab, I invite you to read one of my previous articles where I explain how to set up your Gurobi web license in your Colab environment. You can find the link to the article below:

Optimizing Project Schedules: Using Gurobipy and Cheche_pm warm-starts for Efficient RCPSP…

If you also need a tutorial and more information about linear programming (LP) or integer programming (IP) techniques I strongly recommend you to check these excellent articles written by Bruno Scalia C. F. Leite.

Linear programming: Theory and applications

A Comprehensive Guide to Modeling Techniques in Mixed-Integer Linear Programming

Article Index:

• Installing the libraries and setting up the environment
• Resource planning input data
• Standard ILP formulation for the graph coloring problem
• Column generation
• Column algorithm for the graph coloring problem
• Conclusions
• References

---

Installing the libraries and setting up the Colab environment

To use Gurobi in Google Colab, we first need to install it using the following code:

!pip install gurobipy

Once installed we can proceed to import the required libraries for this project.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from itertools import chain, combinations, permutations, product
import time
from tqdm import tqdm
import gurobipy as gp
from gurobipy import GRB
import networkx as nx
from pprint import pprint

We also need to initialize a gurobi session, so we need to create a params dictionary with all the relevant information of our gurobi license.

WLSACCESSID = '< Copy your WLSACCESSID here >'
WLSSECRET = '< Copy your WLSSECRET here >'
LICENSEID = '< Copy your LICENSEID here >'

params = {
"WLSACCESSID": WLSACCESSID,
"WLSSECRET": WLSSECRET,
"LICENSEID": LICENSEID,
}
env = gp.Env(params=params) #this line initializes the gurobi environment

---

Resource planning input data

For the remaining of the article we will continue to work with the example described by Figure 1. The data that generated this graph is presented in Table 1.

In this section, we will transform this data into a graph representation G(V,E), with nodes V representing the jobs and edges E representing the links between overlapping jobs or jobs with less than three days between them. To accomplish this, we will create an undirected graph using NetworkX with the following code:

data = {
    'Job': ['Job-1', 'Job-2', 'Job-3', 'Job-4', 'Job-5', 'Job-6', 'Job-7', 'Job-8', 'Job-9', 'Job-10'],
    'Start': ['2024-06-01', '2024-06-05', '2024-06-08', '2024-06-10', '2024-06-15', '2024-06-18', '2024-06-20', '2024-06-25', '2024-06-28', '2024-07-01'],
    'End': ['2024-06-10', '2024-06-12', '2024-06-15', '2024-06-20', '2024-06-25', '2024-06-28', '2024-07-05', '2024-07-08', '2024-07-10', '2024-07-15']
}

# Create DataFrame
df = pd.DataFrame(data)

# Convert Start and End to datetime
df['Start'] = pd.to_datetime(df['Start'])
df['End'] = pd.to_datetime(df['End'])

# Create the graph
G = nx.Graph()

# Add nodes
for job in df['Job']:
    G.add_node(job)

# Add edges based on overlapping dates or start date being less than 3 days after the end date of another job
for i in range(len(df)):
    for j in range(i + 1, len(df)):
        job1, job2 = df.iloc[i], df.iloc[j]
        if (job1['Start'] <= job2['End'] and job2['Start'] <= job1['End']) or 
           (abs((job1['End'] - job2['Start']).days) < 3 or abs((job2['End'] - job1['Start']).days) < 3):
            G.add_edge(job1['Job'], job2['Job'])

# Print the edges of the graph to check
print("Edges of the graph:")
print(G.edges())

We can print the Edges to verify that everything is correct. In order to generate a better view of the edges, lets use some pprint( ) instead of the regular print( ) function.

# Print the edges of the graph to check
print("Edges of the graph:")
pprint(list(G.edges()),width=60,compact=True)

>>>>>

Edges of the graph:
[('Job-1', 'Job-2'), ('Job-1', 'Job-3'), ('Job-1', 'Job-4'),
 ('Job-2', 'Job-3'), ('Job-2', 'Job-4'), ('Job-3', 'Job-4'),
 ('Job-3', 'Job-5'), ('Job-4', 'Job-5'), ('Job-4', 'Job-6'),
 ('Job-4', 'Job-7'), ('Job-5', 'Job-6'), ('Job-5', 'Job-7'),
 ('Job-5', 'Job-8'), ('Job-6', 'Job-7'), ('Job-6', 'Job-8'),
 ('Job-6', 'Job-9'), ('Job-7', 'Job-8'), ('Job-7', 'Job-9'),
 ('Job-7', 'Job-10'), ('Job-8', 'Job-9'),
 ('Job-8', 'Job-10'), ('Job-9', 'Job-10')]

After verifying that all edges comply with the data in Table 1 and the established business rules (jobs with less than three days between them are linked via an edge), let’s proceed to plot our graph G. After all, what’s the point of doing all this if we can’t generate some fancy visualizations, right?

plt.figure(figsize=(12, 8))
pos = nx.circular_layout(G)
nx.draw(G,pos= pos, with_labels=True,
               node_size=3000, font_size=12,node_color='lightgrey')
plt.show()

---

Standard ILP formulation for the graph coloring problem

Now that our data is ready, we can discuss integer programming formulations for the graph coloring problem. The literature offers a myriad of different integer programming formulations for this problem, including various families such as assignment formulations [5], representative formulations [6], partial ordering-based formulations [7], and set covering-based formulations [8]. For simplicity, in this section, we will cover only the assignment formulation by Mendez-Diaz and Zabala (2006) [5].

This is a purely binary programming formulation, with binary variables w_i and x_{i,v} ​. Here, iis the index of the potential different colors in set H. You might be wondering how many colors we assume initially; the simple answer is that we naively assume we have as many colors as there are nodes. The variable v is the index for the set of all nodes V.

• Constraint (1) Just ensures that each node has one unique color assigned to it.
• Constraint (2) Makes sure two adjacent nodes, do not share the same color.
• Constraints (3) and (4) are just there for symmetry braking and to increase the efficiency of the branch and bound algorithm.

To apply this formulation to our problem, you can use the code below. It creates the binary programming model and solves it.

E = list(G.edges())
V = list(G.nodes())

H,V = range(len(V)), range(len(V))

x_var = [(v,i) for v in V for i in H ]

model = gp.Model("Binary_ILP",env=env)

W = model.addVars(H,vtype=GRB.BINARY, name="w")
X = model.addVars(x_var ,vtype=GRB.BINARY, name="x")

model.setObjective(gp.quicksum(W[i] for i in H),GRB.MINIMIZE); 

for v in V:
  model.addConstr(gp.quicksum(X[(v,i)] for i in H) == 1 )

for i in H:
  for u,v in E:

    u = int(u.split("-")[1])-1
    v = int(v.split("-")[1])-1

    model.addConstr(X[(v,i)] + X[(u,i)] <= W[i] )

for i in H:
  model.addConstr(gp.quicksum(X[(v,i)] for v in V) >= W[i])

for i in range(1,len(V)):
  model.addConstr(W[i] <= W[i-1])

timelimit = 600
mipgap = 0.0

model.setParam('TimeLimit', timelimit)
model.setParam('MIPGap',mipgap)
model.optimize()

Now that we have found the optimal solution we can proceed to extract the color of each node by using the code below.

colors = ['lightskyblue','tomato','lime','lemonchiffon']
print(f"Total number of colors needed = {model.getObjective().getValue()}")
node_colors = []
for key in X.keys():

  if X[key].x >= 0.99:

    print(f" job-{key[0] +1} is colored with color: {colors[key[1]]}")
    node_colors.append(colors[key[1]])

This will produce the following output.

>>>>

Total number of colors needed = 4.0
 job-1 is colored with color: tomato
 job-2 is colored with color: lemonchiffon
 job-3 is colored with color: lime
 job-4 is colored with color: lightskyblue
 job-5 is colored with color: lemonchiffon
 job-6 is colored with color: lime
 job-7 is colored with color: tomato
 job-8 is colored with color: lightskyblue
 job-9 is colored with color: lemonchiffon
 job-10 is colored with color: lime

We can also proceed to generate to plot of the problem graph G with the coloring solution.

plt.figure(figsize=(12, 8))
pos = nx.circular_layout(G)
nx.draw(G,pos= pos, with_labels=True, node_size=3000, 
                    font_size=12,node_color=node_colors)
plt.show()

This is a great formulation, and for this small toy problem, Gurobi solved it incredibly fast 🔥 . However, if the number of nodes becomes too large, even the best commercial solvers might struggle to find a decent solution. For example, while the simple example of 10 jobs we’re using in this article is manageable, in my previous job, we often handled more than 40 jobs (nodes) in a single month. So, how challenging is a graph coloring problem with 40 nodes and 40 potential colors? We would have 40 binary variables for w_i​ and 1,600 binary variables for x_{i,v}​ (40 nodes times 40 colors). Now, let’s consider the solution space for this problem. How many possibilities do we have? To simplify, let’s focus on x_{i,v}​; each of these variables can be either 0 or 1. So the total combination of possibilities is 2¹⁶⁰⁰. Imagine trying to brute force the solution from that search space? (trying to explore each possibility)

2¹⁶⁰⁰ That number is just astronomically massive, it is in fact multiple times bigger than the number of observable atoms in the universe ( around 10⁸³ ) crazy don’t you think? 😵

Luckily for us, modern solvers do not use a brute force strategy to evaluate each possibility. However, the problem can still be challenging. That is why we are going to explore a powerful technique called column generation and apply it to solve this problem.

---

Column Generation

Column generation is a powerful mathematical optimization technique used primarily in large-scale linear programming problems, particularly in cases where the number of potential variables (columns) is prohibitively large (like in our resource planning problem). Think of column generation as starting with a small set of promising variables and gradually add more by identifying which new variables would improve the solution of the problem the most. The initial idea of column generation was introduced by Gilmore and Gomory [9] in 1961.

So Instead of solving the problem in its entirety (like in the previous section), which would be computationally challenging, we start with a manageable subset of variables and solve a so called "restricted master problem" (RMP).

Then using the solution of the RMP, we then generate additional variables that have the potential to improve the objective function, by solving a subproblem known as the pricing problem. This iterative process continues, dynamically adding the most promising variables, until no further improvements can be made, thus converging to the optimal solution.

In general column generation involves the following steps:

• (1) Solve the Restricted Master Problem (RMP): Begin by relaxing the original large-scale problem, which means solving a simpler, more manageable version of the problem. This initial step reduces computational complexity by not considering all variables at once, and by sometimes considering a different domain from the variables (example instead of binary or integer variables we will relax the problem by assuming continuous variables). Solving the RMP gives us an initial solution and the associated dual values for the constraints. Dual values (or shadow prices) represent the change in the objective function value per unit increase in the right-hand side of the constraints. They provide insight into how much we can improve our solution by relaxing or tightening each constraint. Think of duality as looking at the same problem from two different perspectives, the dual problem looks at the primal (original) problem from the perspective of the constraints.
• (2) Formulate the Pricing Subproblem: The pricing subproblem is formulated to identify new columns (variables) that can potentially improve the objective function of the RMP. The goal here is to find columns that have a negative reduced cost in a minimization problem (or a positive cost in a maximization problem).
• (3) Solve the Pricing Subproblem: Solve the pricing subproblem using the duals obtained from the RMP and searching for the best new column.
• (4) Calculate Reduced Costs: Reduced costs help us determine the potential benefit of adding a new column to the RMP. The reduced cost of a column is calculated as the difference between its original cost and the value of the dual variables corresponding to the constraints it satisfies. Mathematically the reduced cost of adding a new column can be calculated by the following expression.

• (5) Generate New Columns: If the pricing subproblem finds a columng with a negative reduced cost, then this column is added to the RMP. If no columns with negative reduced costs are found, this indicates that the current solution (with the current set of columns) is optimal.
• (6) Iterate: Repeat the process by solving the updated RMP with the new columns, recalculating dual values, and resolving the pricing subproblem. This iterative process continues until no further columns with negative reduced costs are identified, signaling convergence to the optimal solution.

Figure 7. showcases a full diagram of the entire process.

---

Column algorithm for the graph coloring problem

Now, armed with this powerful technique, let’s proceed to implement it for solving the graph coloring problem. We will follow the work of Mehrotra et al. (1996) [8]. Instead of using the assignment formulation by Mendez-Diaz and Zabala [5], which employs binary variables to determine if a node is colored by a given color, we will shift towards a set covering formulation.

In this approach, we consider independent sets S, with each set associated with a unique color. These sets can contain different nodes, but each node can only belong to one set. The objective is to find the minimum number of independent sets required to color a given graph G(V,E). With this formulation, the ILP for our master problem (MP) becomes:

Where y_s is a binary variable that is equal to 1 if the set _sis taken and 0 otherwise. The parameter `b{v,s}is binary and it is equal to 1 f the node _v_ is inside the set _s_ . Constraint (1) just guarantees that a node can only be part of one set _s._ If we relaxed the MP, then the variablesy_s` become continuous variables instead of binary variables. After defining the MP then the pricing subproblem can be defined as.

Where x_v is a binary variable that is equal to 1 if the node v is in the independent set and equal to 0 otherwise. d_v just represent the dual values obtained from the solution of the RMP.

To simplify the starting point of this column generation problem, we will initialize with some easy-to-build independent sets called singleton sets. Singleton sets contain only one node each, so we will start with a number of singleton sets equal to the number of nodes in our problem graph. We will create a class that initializes the data structures needed for our master and pricing problems.

class InitialPatternGenerator:

  def __init__(self,Graph):
    columns = ['independentSetCost','independentSet']
    self.G = Graph
    self.nodes = list(self.G.nodes())
    self.edges = list(self.G.edges())
    self.n_Nodes = len(self.nodes)
    IndepndentSet = pd.DataFrame(index = range(self.n_Nodes),columns = columns)
    self.IS_DF = IndepndentSet

  def generateBasicSingletons(self):

    self.IS_DF['independentSetCost'] = 1
    self.IS_DF['independentSet'] = [np.where(self.IS_DF.index == j,1,0) for j in range(self.n_Nodes)]
    return self.IS_DF

  def generate_dictionaries(self):

    name_to_number = dict()
    number_to_name = dict()

    for i,v in enumerate(self.nodes):
      name_to_number[v] = i
      number_to_name[i] = v

    self.name_to_number_dict = name_to_number
    self.number_to_name_dict = number_to_name

  def get_numeric_v_and_e(self):

    edges = []

    for u,v in self.edges:

      u = self.name_to_number_dict[u]
      v = self.name_to_number_dict[v]

      edges.append((u,v))

    nodes = [self.name_to_number_dict[x] for x in self.nodes]

    return nodes, edges

The class above called InitialPatternGenerator generates the different data structures that are required, from the problem graph G. Specifically the method generateBasicSingletons produces a dataframe with the singleton independent sets for each node v . The code below, prepares all of our initial problem inputs and parameters.

IS_generator = InitialPatternGenerator(G)
IS_generator.generate_dictionaries()

nodes,edges = IS_generator.get_numeric_v_and_e()

IS_DF = IS_generator.generateBasicSingletons()

Now that we have all the required inputs and parameters we will proceed to implement the class for the master problem by using the code below.

class MasterProblem:

  def __init__(self,IS_DF,nodes,edges,env):

    self.independentSetCost = IS_DF['independentSetCost'].values
    self.independentSet = IS_DF['independentSet'].values
    self.nodes = nodes
    self.edges = edges
    self.n_Nodes = len(nodes)
    self.model = gp.Model("MasterProblem",env=env)
    self.independentSetIndex = IS_DF.index.values

  def build_model(self):
    self.GenerateVariables()
    self.GenerateConstraints()
    self.GenerateObjective()
    self.model.update()

  def GenerateVariables(self):
    self.independentSetUseVar = self.model.addVars(self.independentSetIndex,
                                            vtype = GRB.BINARY,name = "independentSetUseVar")
  def GenerateConstraints(self):
    for i in range(self.n_Nodes):
      self.model.addConstr(gp.quicksum(self.independentSetUseVar[p]*self.independentSet[p][i] for p in self.independentSetIndex) >= 1, f'C_{i}')

  def GenerateObjective(self):
    self.model.setObjective(gp.quicksum(self.independentSetUseVar[p] for p in self.independentSetIndex),sense = GRB.MINIMIZE)

  def solveRelaxModel(self,silence=True):

    self.relaxedModel = self.model.relax()
    if silence:
      self.relaxedModel.setParam('LogToConsole',0)
    self.relaxedModel.optimize()

  def getDuals(self):

    return self.relaxedModel.getAttr("Pi",self.relaxedModel.getConstrs())

  def addColumn(self,objective,new_IS):

    ctName = f'independentSetUseVar[{len(self.model.getVars())}]'
    new_column = gp.Column(new_IS,self.model.getConstrs())
    self.model.addVar(vtype = GRB.INTEGER,lb = 0, obj = objective, 
                              column = new_column,name = ctName)
    self.model.update()

  def solveFinalModel(self,time_limit = 3600,mipgap = 0.05,silence = True):

    self.model.setParam('TimeLimit',time_limit)
    self.model.setParam('MIPGap',mipgap)
    if silence:
       self.model.setParam('LogToConsole',0)

    self.model.optimize()

  def extractSolution(self,silence = True):

    all_vars = self.model.getVars()
    values = self.model.getAttr("X", all_vars)
    names = self.model.getAttr("VarName", all_vars)

    solution = list()

    if not silence:
      print(f"Total number of colors needed = {self.model.getObjective().getValue()}")

    for name,x,pattern in zip(names,values,NEW_independent_sets):
      if x >= 0.99:
        indices = [index for index, value in enumerate(pattern) if value == 1.0]
        solution.append(indices)
    output = list()
    for c,i_S in enumerate(solution):
      for x in i_S:
        output.append((x,c))
        if not silence:
          print(f"job-{x+1} assigned to color {c+1}")

    return output

The MasterProblem class defined above contains the following methods:

• The constructor: Initializes the class with the input data, from a list of nodes, edges, our gurobi environment and the singleton pattern dataframe.
• buildomel: This method builds the entire model by calling other methods related to the creation of variables, objective function and constraints of the formulation of the MP.
• GenerateVariables: This method adds the variables to the model.
• GenerateConstraints: This method adds the constraints to the model.
• GenerateObjective: This method adds the objective function to the model.
• solveRelaxModel: This method relaxes the master problem (allows the variable y_s to become continuous instead of binary.
• getDuals: This method obtains the dual values from the solution of the relaxed problem.
• addColumn: This method adds a new column to the master problem.
• solveFinalModel: This method is reserved to be used at the end of the optimization process once all the relevant columns have been found.
• extractSolution: This method extracts the solution from the independent set columns, after solving the master problem without relaxing it.

With the master problem successfully implemented we can now proceed to implement the class of the pricing subproblem by following the snippet of code below.

class SubProblem:

  def __init__(self,IS_DF,nodes,edges,env,duals):

    self.independentSetCost = IS_DF['independentSetCost'].values
    self.independentSet = IS_DF['independentSet'].values
    self.nodes = nodes
    self.edges = edges
    self.n_Nodes = len(nodes)
    self.duals = duals
    self.model = gp.Model("SubProblem",env=env)
    self.nodesIndex = IS_DF.index.values

  def buildModel(self):
    self.GenerateVariables()
    self.GenerateConstraints()
    self.GenerateObjective()
    self.model.update()

  def GenerateVariables(self):
    self.nodePatternVar = self.model.addVars(self.nodesIndex,
                                                 vtype = GRB.BINARY,name = "nodePatternVar")
  def GenerateConstraints(self):

    for u,v in self.edges:

      self.model.addConstr(self.nodePatternVar[u] + self.nodePatternVar[v] <= 1)

  def GenerateObjective(self):
    self.model.setObjective(gp.quicksum(self.nodePatternVar[p]*self.duals[p] for p in self.nodesIndex),sense = GRB.MAXIMIZE)

  def getNewPattern(self):
    return self.model.getAttr("X",self.model.getVars())

  def solveModel(self,time_limit = 3600,mipgap = 0.05,silence = True):

    self.model.setParam('TimeLimit',time_limit)
    self.model.setParam('MIPGap',mipgap)
    if silence:
       self.model.setParam('LogToConsole',0)
    self.model.optimize()

  def getObjectiveValue(self):
    return self.model.getObjective().getValue()

Very similar as with the class of the master problem, the class pricing SubProblem defined above, has the following methods:

• The constructor: Initializes the class with the input data, from a list of nodes, edges, our gurobi environment, the singleton pattern dataframe and a list of the dual values coming from the master problem.
• The methods buildModel, GenerateVariables, GenerateConstraints and GenerateObjective work exactly the same as with the master problem.
• solveModel: This method solves the pricing subproblem
• getNewPattern: This method generates a new column from the solution of the pricing subproblem.
• getObjectiveValue: This method retrieves the objective value from the solution of the pricing subproblem.

Now that all the classes are implemented we can finally the proceed to build the main loop for solving our problem. We will be following the exact procedure described by Figure 7 , and following the code below.

master = MasterProblem(IS_DF,nodes,edges,env)
master.build_model()

NEW_independent_sets = [x.tolist() for x in IS_DF['independentSet'].values]

modelImprovable = True
counter = 0

for i in range(20):

  print(f"Iteration: {counter}")

  master.solveRelaxModel()
  duals = master.getDuals()

  subproblem = SubProblem(IS_DF,nodes,edges,env,duals)
  subproblem.buildModel()

  subproblem.solveModel(120,0.0,True)

  modelImprovable = 1 - subproblem.getObjectiveValue() < 0

  print("gain",1 - subproblem.getObjectiveValue())

  if modelImprovable:

    new_independent_set_cost = 1

    new_independent_set = subproblem.getNewPattern()
    print(new_independent_set)

    NEW_independent_sets.append(new_independent_set)

    master.addColumn(new_independent_set_cost,new_independent_set)

    counter+=1
    print("-"*120)
  else:
    break

master.solveFinalModel(3600,0.00)

After few iterations the column generation algorithm halts and we obtain the optimal solution that only requires 4 colors (independent sets), to color the problem graph.

master.extractSolution(silence = False)

Total number of colors needed = 4.0
job-3 assigned to color 1
job-7 assigned to color 1
job-4 assigned to color 2
job-8 assigned to color 2
job-2 assigned to color 3
job-6 assigned to color 3
job-10 assigned to color 3
job-1 assigned to color 4
job-5 assigned to color 4
job-9 assigned to color 4

[(2, 0), (6, 0), (3, 1), (7, 1),
 (1, 2), (5, 2), (9, 2), (0, 3),
 (4, 3), (8, 3)]

We can then proceed to plot the solution into the problem graph G.

colors = ['lightskyblue','tomato','lime','lemonchiffon']
node_colors = [None for x in range(len(nodes))]
for x,c in problem_solution:
  node_colors[x] = colors[c]

plt.figure(figsize=(12, 8))
pos = nx.circular_layout(G)
nx.draw(G,pos= pos, with_labels=True, node_size=3000, font_size=12,node_color=node_colors)
plt.show()

Conclusions

In conclusion, this article covered multiple applications of the graph coloring problem and its relevance to resource planning problems. We discussed the assignment binary programming formulation of Mendez-Diaz and Zabala, the set covering formulation of the same problem, and the implementation of a powerful column generation algorithm published by Mehrotra et al. This powerful technique is particularly useful for large problem instances (for smaller problem instances the assignment formulation of Mendez-Diaz and Zabala is more than enough). I hope this article serves as a helpful resource for anyone facing similar challenges in their current jobs. I would have loved 💙 to have had this knowledge back then, but in any case it’s never too late to learn 😄 . You can find the entire notebook with the code developed for this article in the link to my GitHub repository jsut below.

Youtube_optimization/Column_Generation_graph_coloring_problem_gurobipy_medium_example_for_sharing.ip…

I sincerely hope you found this article useful and entertaining. If so, I’d love to hear your thoughts! Please feel free to leave a comment or show your appreciation with a clap 👏 . And if you’re interested in staying updated on my latest articles, consider following me on Medium. Your support and feedback are what drive me to keep exploring and sharing in this fascinating field. Thank you for taking the time to read, and stay tuned for more insights in my next article!

---

References

• [1] Brooks, R. L. (1941), "On coloring the nodes of a network", _Proceedings of the Cambridge Philosophical Society_, 37 (2): 194–197
• [2] Francis Guthrie. (June 10, 1854), "Tinting Maps", [The Athenaeum](https://en.wikipedia.org/wiki/Athenaeum(Britishmagazine)): 726
• [3] Marx, Dániel (2004), "Graph colouring problems and their applications in scheduling", Periodica Polytechnica, Electrical Engineering, vol. 48, pp. 11–16
• [4] Chaitin, G. J. (1982), "Register allocation & spilling via graph colouring", _Proc. 1982 SIGPLAN Symposium on Compiler Construction_, pp. 98–105
• [5] Méndez-Díaz, Isabel, and Paula Zabala. "A branch-and-cut algorithm for graph coloring." Discrete Applied Mathematics 154.5 (2006): 826–847.
• [6] Campêlo, Manoel, Ricardo Corrêa, and Yuri Frota. "Cliques, holes and the vertex coloring polytope." Information Processing Letters 89.4 (2004): 159–164.
• [7] Jabrayilov, A., & Mutzel, P. (2018). New integer linear programming models for the vertex coloring problem. In LATIN 2018: Theoretical Informatics: 13th Latin American Symposium, Buenos Aires, Argentina, April 16–19, 2018, Proceedings 13 (pp. 640–652). Springer International Publishing.
• [8] Mehrotra, Anuj, and Michael A. Trick. "A column generation approach for graph coloring." informs Journal on Computing 8.4 (1996): 344–354.
• [9] Gilmore, Paul C., and Ralph E. Gomory. "A linear programming approach to the cutting-stock problem." Operations Research 9.6 (1961): 849–859.