To refine our solutions obtained using heuristics and attempt to prove the optimality of the solutions, let us formulate the graph coloring problem as an ILP model. Keep in mind that it might not be able to handle large instances though. The model presented in this section and other exact algorithms are presented in *chapter 3* by Lewis (2021).

Let’s define it *Sets *taken into account in this approach:

*VS*: Colors*NOT*: Nodes (or vertices)*E*: Edges

A first question already arises: “How many colors should we consider?” “. In the worst case, all nodes are connected, so one must consider *VS* the same size as *NOT*. 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 get an upper limit for the number of colors.

In this problem, we have two groups of decision variables:

*X_*{*I*,*vs*}: Are binary variables indicating this node*I*is colored in*vs**y_*{*vs*}: Are binary variables indicating this color*vs*is used

We must also create constraints to ensure that:

- Each node must be colored
- If a node on an edge has a color, make sure the color is used
- At most 1 node on each edge can be colored with a given color
- Break the symmetry (this is not required, but it can help)

Finally, our goal should minimize the total number of colors used, which is the sum of *Yes*. In summary, we have the following equations.

Without further ado, let’s import *pyomo* for the entire programming model.

`import pyomo.environ as pyo`

There are two approaches to modeling a problem in *pyomo*: *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. You can learn more about these approaches in the library documentation or in the book by Bynum et al. (2021). Throughout this article, we will adopt the *Concrete* formulation of the model.

`model = pyo.ConcreteModel()`

Next, let’s instantiate our *Sets*. Parse iterables directly from `dsatur`

attributes `N`

And `C`

is valid, we could therefore use them in the `initialize`

keyword arguments. Alternatively, I will forward the original *knots *And *edges *from our input data and create a range from the *DSatur* like our initialization for 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 constraints.

`# 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 symmetry breaking constraints and see what works best with your available solver. In some experiments I’ve done, including a preference ordering using the total number of nodes assigned to a given color has proven to be worse than ignoring it. Perhaps 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 objective…

`# 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 this, I use the HiGHS persistent solver, available in *pyomo* in case *spying* 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 might struggle trying to improve heuristic solutions. However, for a 32 node instance available in 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).