To refine our solutions obtained using heuristics and try to prove the optimality of the solutions, let us formulate the graph coloring problem as an ILP model. Remember that it may not be able to handle large instances. The model presented in this section and other exact algorithms are presented in Chapter 3 by Lewis (2021).
Let’s define the Sets considered in this approach:
- c: Colors
- north: Nodes (or vertices)
- my: Edges
A first question already arises: “How many colors should we consider?” In the worst case, all nodes are connected, so it must be considered c the same size as north. However, this approach could make our solutions even more difficult by unnecessarily increasing the number of decision variables and worsening the linear relaxation of the model. A good alternative is to use a heuristic method (such as DSatur) to obtain an upper limit for the number of colors.
In this problem, we have two groups of decision variables:
- X_{Yo, c}: These are binary variables that indicate that the node Yo is colored in c
- and_{c}: They are binary variables that indicate that color. c used
We must also create constraints to ensure that:
- Each node must be colored.
- If any node on an edge has a color, make sure the color is being used
- At most 1 node on each edge can be colored with a given color.
- Break the symmetry (this is not necessary, but might help)
Finally, our goal should minimize the total number of colors used, which is the sum of and. Summarizing we have the following equations.
Without further ado, let’s import piomo for the integer programming model.
import pyomo.environ as pyo
There are two approaches to modeling a problem in piomo: Abstract and Concrete Models. In the first approach, the algebraic expressions of the problem are defined before some data values are provided, while in the second, the model instance is created immediately as its elements are defined. More information about these approaches can be found in the library documentation or in the book by Bynum et al. (2021). Throughout this article, we will adopt the Concrete model formulation.
model = pyo.ConcreteModel()
So, we instantiate our Sets. Analyzing the iterables directly from dsatur
attributes N
and C
is valid, so they could be used in the initialize
keyword arguments. Alternatively, I’ll pass on the original. nodes and edges from our input data and create a range from the DSatur as our initialization for the colors.
model.C = pyo.Set(initialize=range(len(dsatur.C))) # Colors
model.N = pyo.Set(initialize=nodes) # Nodes
model.E = pyo.Set(initialize=edges)) # Edges
Next, we instantiate our decision variables.
model.x = pyo.Var(model.N, model.C, within=pyo.Binary)
model.y = pyo.Var(model.C, within=pyo.Binary)
And then our limitations.
# Fill every node with some color
def fill_cstr(model, i):
return sum(model.x(i, :)) == 1# Do not repeat colors on edges and indicator that color is used
def edge_cstr(model, i, j, c):
return model.x(i, c) + model.x(j, c) <= model.y(c)
# Break symmetry by setting a preference order (not required)
def break_symmetry(model, c):
if model.C.first() == c:
return 0 <= model.y(c)
else:
c_prev = model.C.prev(c)
return model.y(c) <= model.y(c_prev)
model.fill_cstr = pyo.Constraint(model.N, rule=fill_cstr)
model.edge_cstr = pyo.Constraint(model.E, model.C, rule=edge_cstr)
model.break_symmetry = pyo.Constraint(model.C, rule=break_symmetry)
You can try including other constraints that break the symmetry and see what works best with the available solver. In some experiments I’ve run, including a preference order using the total number of nodes assigned to a given color has been worse than ignoring it. Possibly due to the solver’s native symmetry breaking strategies.
# Break symmetry by setting a preference order with total assignments
def break_symmetry_agg(model, c):
if model.C.first() == c:
return 0 <= sum(model.x(:, c))
else:
c_prev = model.C.prev(c)
return sum(model.x(:, c)) <= sum(model.x(:, c_prev))
Finally, the goal…
# Total number of colors used
def obj(model):
return sum(model.y(:))model.obj = pyo.Objective(rule=obj)
And our model is ready to be solved! To do that I am using the HiGHS persistent solver, which is available at piomo in case spy It is also installed in your Python environment.
solver = pyo.SolverFactory("appsi_highs")
solver.highs_options("time_limit") = 120
res = solver.solve(model)
print(res)
For large instances, our solver may have difficulty trying to improve heuristic solutions. However, for a 32-node instance available on the code repositorythe solver was able to reduce the number of colors used from 9 to 8. It is worth mentioning that it took 24 seconds to complete the execution, while the DSatur The algorithm for the same instance took only 6 milliseconds (using pure Python, which is an interpreted language).