Mathematics of Love: Optimizing a Dining-Room Seating Arrangement for Weddings with Python

Dr. Hannah Fry‘s book, The Mathematics of Love [1], is one of those rare finds, smart, funny, and incredibly easy to read. I enjoyed it so much that I ended up giving copies to some of my closest friends. They loved it too, but a few of them hit a snag with Chapter 8, which, ironically, was the chapter they were most excited about.

This chapter goes into the tricky business of organizing wedding seating arrangements, using mathematical programming and optimization to figure out how to seat guests so everyone has the best time possible. Sounds cool, right? But here’s the catch, the book doesn’t actually walk you through how to set up or solve the problem, which left my friends feeling a bit lost.

Now, I know what you’re thinking: wedding seating? Really? But don’t be fooled , it’s a really difficult problem to solve, and the solution has way more applications than just weddings (although this one is already important). Think organizing tables on cruise ships [3], forming sports teams or workgroups, optimizing portfolios, and in my world as a business school professor, putting together effective teams for group projects (group conflicts and free riding have dramatically decreased since I started organizing groups with the methodology described in this article). So yeah, it’s a problem worth knowing how to solve, at least in my humble opinion.

In this article, I’m going to pick up where Dr. Fry left off and show you how to actually tackle this problem using Python and Gurobi. We’ll start by breaking down the problem, then jump into the details of it’s Mixed Integer Programming (MIP) formulation. After that, we’ll bring in a sample dataset and solve it step by step. And to wrap things up, I’ll talk about some cool extensions and advanced methods you can use to take this even further.

Now If you have never used Gurobi before, especially from Google Colab (The IDE that we will use in this article), 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, I am leaving the links down below. Another great source is William’s HP book about Mathematical Programming [2].

Linear programming: Theory and applications

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

Article index:

• The Problem
• MILP formulation of the RQMKP
• Installing the libraries and setting up the Colab environment
• Problem Data
• Solving the RQMKP
• Visualizing the RQMKP solution
• Conclusions
• References

The Problem

As Dr. Hanna Fry puts it, aside from choosing a terrible best man to give the speech, one of the biggest mistakes a couple can make at their wedding is seating enemies at the same table (In my home country not having enough tequeños would also be considered a mortal sin). Make that mistake, and you’re in for a night of icy glares, passive-aggressive comments, and, if alcohol is involved, maybe even a brawl. Getting the seating arrangement right can make a world of difference in the overall vibe of the big day.

Every optimization problem has an objective function, and in this case, our goal is to maximize the total happiness of the guests. Guest happiness depends on how they interact with each other. Let’s take an example: Pongo and Cruella. Cruella loves Pongo’s beautiful hair, so the interaction term from Cruella’s side is +1. But it’s not always mutual – Pongo, worried about Cruella’s obsession with his fur, feels uneasy, so his interaction term is -1. If they’re seated together, their happiness cancels out to zero. Now, imagine doing this for 100 guests instead of just two. How much harder does it get? If you try to brute force the Sh-t out of this, you’ll find that there are 100 factorial (100!) ways to arrange the seating – that’s an astronomical number. This is precisely why mathematical optimization is our best bet.

To make matters more complex, we have constraints to consider. It’s not just about maximizing happiness; we also need to respect the table capacity limits and keep couples together (they tend to be angry if they are separated). So, our problem is to find an assignment of guests to tables that maximizes total happiness while respecting these constraints.

Now, here’s a philosophical twist: Should we aim to maximize the total happiness of the entire venue, or should we try to maximize the happiness of the least happy table? Which approach is fairer to our guests? What if maximizing total happiness makes one table absolutely miserable? Don’t worry, this article will tackle both approaches, but I am leaving the question open as food for thought.

These kinds of problems are known in the literature as quadratic assignment problems. However, the specific problem we’re solving here can be addressed as a Restricted Quadratic Multiknapsack Problem (RQMKP).

You’ve probably heard of the Binary Knapsack problem – there are tons of articles on Medium about how to solve it, I’ll leave a link to a cool Knapsack article written by Maria Mouschoutzi, PhD below.

How Many Pokémon Fit?

But the RQMKP? That’s what the Knapsack problem evolves into after eating its vegetables, getting eight hours of sleep, and hitting the gym five times a week. This isn’t just a Knapsack; it’s "the" Knapsack problem – the Goku Ultra Instinct version of it. It’s been proven to be strongly NP-Hard [6], which makes it a real beast to solve.

Unlike the regular Knapsack problem, the RQMKP involves multiple ("M") knapsacks that need to be filled. The objective function includes both quadratic ("Q") and linear terms, which, in our case, represents happiness – a function that, I should remind you, we need to maximize. The "R" stands for restricted, meaning we must place each item into one of the available knapsacks without exceeding their individual capacities, but every item has to be placed somewhere.

In our wedding dining context, when we focus on maximizing the total happiness of the event, we’re dealing with the Max-Sum variant of the RQMKP. When we aim to maximize the happiness of the least happy table, we’re looking at the Max-Min variation (we are maximizing the minimum of the tables).

Now that we have a broad understanding of the problem, let’s discuss the formal mathematical programming formulations.

Mathematical formulation of the RQMKP

We’ll start out with the simpler version of the problem, the Max-Sum RQMKP. Once we’ve tackled that, we’ll move on to the more challenging Max-Min variation.

The Max-Sum RQMKP

A mathematical formulation of the Max-Sum RQMKP problem is presented in Figure 1 below.

At first glance, the formulation above might seem daunting, but it’s actually quite straightforward. It uses binary variables x_{i,k}, which equal 1 if guest i is seated at table k, and 0 otherwise. The objective function is simply the sum of the interaction terms c_{i,j} for all pairs of individuals {i,j} who are seated together at the same table.

1. Constraint (1) ensures that for each individual i, only one binary variable x_{i,k} equals 1, meaning each person can be seated at only one table.
2. Constraint (2) makes sure that the number of guests at table k does not exceed its maximum capacity B_k. Note that different tables might have different capacities, which is why there’s a specific parameter B_k for each table.
3. Constraint (3) ensures that all couples {i,j} in the set of couples E are seated together at the same table k.
4. Finally, Constraint (4) guarantees that the variables x_{i,k} remain binary.

The Max-Min RQMKP

A mathematical formulation of the Max-Min RQMKP problem is presented in Figure 2 below. This formulation is slightly more difficult, but nonetheless easy to grasp.

This formulation is very similar to the Max-Sum-RQMKP, but here, the objective function is a new variable ζ that we are trying to maximize. Constraint (1) forces ζ to be less than or equal to the total happiness of each individual table. This means that, in the best-case scenario, ζ represents the happiness of the least happy table.

The rest of the constraints (2–5) are exactly the same as in the Max-Sum case, so there’s no need to repeat them.

With the formulations in place, let’s proceed to set up the environment to solve these problems.

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 networkx as nx
import matplotlib.pyplot as plt
import gurobipy as gp
from gurobipy import GRB
from google.colab import userdata

We’re going to use the standard "mirepoix" (🥬🥕🧅) of almost any Python project: numpy, pandas, and matplotlib. Additionally, we’ll import networkx to visualize the relationships between our guests, and gurobipy to formulate and solve the problem. With all the libraries in place, the only thing left to do is initialize the Gurobi environment. You can easily do this with the code below. Please note that you’ll need your Gurobi API credentials to proceed. I’m storing mine in the "sec_rets" t_ab of Google Colab, but if you prefer, you can directly copy the API access values into the variables.

WLSACCESSID = userdata.get('gurobi_WLSACCESSID')
WLSSECRET = userdata.get('gurobi_WLSSECRET')
LICENSEID = int(userdata.get('gurobi_LICENSEID'))

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

Now that we’re familiar with the problem formulations and have the environment set up, we can move forward with solving the problem. But first, let’s take a look at the example data we’ll be working with.

Problem data

For this article, I’m going to use a small example with 21 guests for our wedding event. Each guest has previously expressed their preference to dine with another guest: a +1 if they want to share a table, a 0 if they are neutral, or a -1 if they want to avoid that person. The data for this example is shown in Table 1 below.

For this example, we’re going to consider only three couples: Pikachu and Meditite, Machoke and Tyrantrum, and Aipom and Pelipper. Remember, we need to seat these couples together.

To keep things simple, we’ll have just three tables, each with a maximum capacity of 7 guests. The data for this problem is available in my GitHub repo. Let’s go ahead and import it using the code snippet below.

url = 'https://raw.githubusercontent.com/ceche1212/math_love_RQMKP_medium/main/input_data.csv'
df = pd.read_csv(url,index_col=0)
idx_to_col = {i:v for i,v in enumerate(df.columns)}
col_to_idx = {v:i for i,v in enumerate(df.columns)}

data = df.to_numpy() # trnasform the data into a numpy array

The data is correctly imported, so let’s go ahead and generate a cool network graph to visualize the relationships between our guests.

G = nx.DiGraph()

nodes = df.columns.to_list()
G.add_nodes_from(nodes)

for i,n1 in enumerate(nodes):
  for j,n2 in enumerate(nodes):
    if i == j:
      continue
    else:
      c_i_j = data[i,j]

      if  c_i_j > 0:
        G.add_edge(n1,n2,weight=c_i_j,color = 'b')
      else:
        G.add_edge(n1,n2,weight=c_i_j,color = 'r')

def nudge(pos, x_shift, y_shift):
    return {n:(x + x_shift, y + y_shift) for n,(x,y) in pos.items()}

plt.figure(figsize=(16,16))
pos = nx.circular_layout(G)
pos_nodes = nudge(pos, 1, 1)
edges = G.edges()
colors = [G[u][v]['color'] for u,v in edges]
nx.draw(G, pos_nodes, edge_color=colors,node_color='grey', alpha=0.4, node_size=3600)
nx.draw_networkx_labels(G, pos=pos_nodes)
plt.title('Wedding Guests')
plt.show()

If everything goes according to plan, you should be looking at something similar to Figure 3 just below.

The graph above looks pretty (and complicated). Each red line corresponds to a negative interaction – for instance, very few people want to eat with Pikachu. I mean, he can be such a diva, and his conversations are just one note, so it’s understandable. In contrast, blue lines represent positive interactions. For instance despite his past work experiences with Team Rocket, and his intransigence regarding pollution and emissions, Koffing seems to be quite popular with the other guests.

To finalize the input data setup, we are also going to need the variables associated with table capacities and couples. The table capacities B_k is going to be a list with the capacity of each table. Couples are going to be stored as tuples inside a list E.

n_tables = 3
B_k = [7,7,7]
E = [(col_to_idx['Pikachu'],col_to_idx['Meditite']),
           (col_to_idx['Machoke'],col_to_idx['Tyrantrum']),
           (col_to_idx['Aipom'],col_to_idx['Pelipper'])]

I,K = range(len(col_to_idx)),range(n_tables)

Solving the RQMKP

Let’s start with the Max-Sum-RQMKP problem. Typically, to solve such models in commercial mathematical programming solvers, we would first need to linearize the quadratic terms x_{i,k}x_{j,k}. Normally this will require the creation of new binary variables w_{i,j,k} = x_{i,k}x_{j,k} that represent the quadratic terms of our formulation, the linearization also required the addition of three new set of constraints, for each i,j and k . This would make the problem substantially more difficult for the solver to solve. But luckily for us, Gurobi has a quadratic solver in place (That is very likely doing some linearization’s behind the scenes), so we can go ahead and directly create the model, as shown in Figure 1.

x_vars = [(i,k) for i in I for k in K] # create indices for variables

model = gp.Model('Max_Sum_RQMKP',env=env) #initialize model

X = model.addVars(x_vars,vtype=GRB.BINARY, name="X") # create variables

#objective function
model.setObjective( gp.quicksum( data[(i,j)]*X[i,k]*X[j,k] for i in I 
                                for j in I 
                                for k in K), GRB.MAXIMIZE ) 

# Constraint (1)
model.addConstrs( gp.quicksum( X[i,k] for k in K ) == 1 for i in I)

# Constraint (2)
model.addConstrs( gp.quicksum( X[i,k] for i in I ) <= B_k[k]  for k in K)

# Constraint (3)
for i,j in E:
  for k in K:
    model.addConstr(X[i,k] == X[j,k])

To solve the model, we just need to set up the MIP gap (stopping criteria based on the percentage distance from the optimal solution) and the time limit we’ll allow the solver to run. For this example, I want to run the solver to optimality if possible (so the MIP gap is equal to 0.00), and we’ll allow Gurobi to run for a maximum of 600 seconds.

tl = 600
mip_gap = 0.00

model.setParam('TimeLimit', tl)
model.setParam('MIPGap', mip_gap)
model.optimize()

After less than 1 second, the solver should stop and provide a final solution with a total happiness of 59. Now, let’s proceed to extract a detailed solution.

max_sum_tables = [ []  for k in K]
for x in X.keys():
  if X[x].x > 0.95:
    guest_idx = x[0]
    table = x[1]
    guest_name = idx_to_col[guest_idx]
    max_sum_tables[table].append(guest_idx)

The table organization (which guest is seated where) is now stored in the list max_sum_tables, now I wonder what is the happiness of each one of the tables? let’s find out.

for i,table in enumerate(max_sum_tables):

  Z = 0
  table_comp = [idx_to_col[x] for x in table]
  for n1 in table:
    for n2 in table:
      Z += data[n1,n2]
  print(f" Table {i+1} ({table_comp}), Hapiness  = {Z}")

The output generated by the code above is the following one:

>>>>

Table 1 Happiness  = 3
Table 2 Happiness  = 23
Table 3 Happiness  = 33

If you’re seated at Table 3, you’re in for a great time, but if you’re at Table 1, you’d better hope the food is good because the conversations might be pretty dull. This configuration maximizes the overall happiness of the event, but it does so at the expense of Table 1’s enjoyment, creating a significant imbalance.

Is it possible to find a more equitable arrangement that still maximizes happiness? The answer is yes, but it will likely come at the cost of the total overall happiness. We can achieve this by solving the Max-Min variation of the RQMKP.

x_vars = [(i,k) for i in I for k in K]

model = gp.Model('Max_Min_RQMKP',env=env) #initialize model

X = model.addVars(x_vars,vtype=GRB.BINARY, name="X") # create variables
Zeta = model.addVar(vtype=GRB.INTEGER, name="zeta") # create objective function variable

# Objective function
model.setObjective( Zeta, GRB.MAXIMIZE )

# Constraint (1)
model.addConstrs(gp.quicksum(X[i,k]*X[j,k]*data[(i,j)] for i in I for j in I) >= Zeta for k in K)

# Constraint (2)
model.addConstrs( gp.quicksum( X[i,k] for k in K ) == 1 for i in I)

# Constraint (3)
model.addConstrs( gp.quicksum( X[i,k] for i in I ) <= B_k[k]  for k in K)

# Constraint (4)
for i,j in E:
  for k in K:
    model.addConstr(X[i,k] == X[j,k])

Again we input the Mip-Gap and time limit, similar with the previous model and then we proceed to solve.

tl = 600
mip_gap = 0.00

model.setParam('TimeLimit', tl)
model.setParam('MIPGap', mip_gap)
model.optimize()

In less than 3 seconds (this variation is a little bit harder to solve), the solver will stop and give you a final solution of 14. But remember, this isn’t the overall happiness of the entire event – it’s just the happiness level of the saddest 😢 table. Now, let’s check out the total happiness of the event and see how each table is feeling individually.

max_min_tables = [ []  for k in K]
for x in X.keys():
  if X[x].x > 0.95:
    guest_idx = x[0]
    table = x[1]
    guest_name = idx_to_col[guest_idx]
    max_min_tables[table].append(guest_idx)

Total_Z = 0
for i,table in enumerate(max_min_tables):

  Z = 0
  table_comp = [idx_to_col[x] for x in table]
  for n1 in table:
    for n2 in table:
      Z += data[n1,n2]
  Total_Z += Z
  print(f" Table {i+1} ({table_comp}), Happiness  = {Z}")
print(f"Total Happiness = {Total_Z}")

The output generated by the code above is the following one:

Table 1 Happiness  = 14
Table 2 Happiness  = 18
Table 3 Happiness  = 14
Total Happiness = 46

The total happiness of this solution is 46, which is a 22% drop compared to the 59 happiness value from the Max-Sum solution. However, the Max-Min solution is significantly more equitable. Just look at Table 1 – their happiness jumped from a dismal 3 to a solid 14. In my opinion, this is a better outcome for everyone involved (assuming you actually like all your guests and want them to have an equally good time. If you’re not too fond of the folks at Table 1, feel free to let them stew in their misery with the Max-Sum configuration – but that begs the question, why invite them at all?).

Visualizing the RQMKP solution

Now that we have our solution, let’s create a table plan to visualize the dining arrangement before wrapping up this article. We’ll use the plotting capabilities of networkx once again. For this example, we’ll plot the solution from the Max-Min variation of the RQMKP, but the process is exactly the same for the Max-Sum problem.

# Create and populate the Table-Guest network graph
Tables_G = nx.Graph()

for i,table in enumerate(max_min_tables):
  table_comp = [idx_to_col[x] for x in table]
  table_name = f'Table {i+1}'
  Tables_G.add_node(table_name,type='table')

  for guest in table_comp:
    Tables_G.add_node(guest,type = 'guest')
    e = (table_name,guest)
    Tables_G.add_edge(*e)

# Find table nodes and guest nodes
table_nodes = [node for node, data in Tables_G.nodes(data=True) if data['type'] == 'table']
guest_nodes = [node for node, data in Tables_G.nodes(data=True) if data['type'] == 'guest']

# Adjust layout for the position of the tables and guests
pos = nx.circular_layout(Tables_G)

# For each table node, create a circular layout for its guests
for table in table_nodes:
    # Get the neighbors (guests) of the table
    guests = list(Tables_G.neighbors(table))

    # Calculate circular layout for guests around the table node
    num_guests = len(guests)
    angle_step = 2 * np.pi / num_guests
    radius = 0.5  # Adjust the radius to space guests further from the table

    for i, guest in enumerate(guests):
        angle = i * angle_step
        pos[guest] = pos[table] + np.array([radius * np.cos(angle), radius * np.sin(angle)])

# Draw the graph with differentiated node sizes
plt.figure(figsize=(10, 10))

# Draw the table nodes
nx.draw_networkx_nodes(Tables_G, pos, nodelist=table_nodes, node_color='lightgrey', node_size=3600, label='Table',)

# Draw the guest nodes
nx.draw_networkx_nodes(Tables_G, pos, nodelist=guest_nodes, node_color='lightblue', node_size=1000, label='Guest')

# Draw the edges and labels
nx.draw_networkx_edges(Tables_G, pos)
nx.draw_networkx_labels(Tables_G, pos, font_size=10)
plt.title('Table Configuration Max-Min')
plt.show()

The resulting plot is displayed in Figure 4 below.

Conclusions

In conclusion, this article has guided you through solving a dining-room seating arrangement problem by tackling it as a Restricted Quadratic Multi-Knapsack Problem (RQMKP). We explored two different formulations: the Max-Sum version and the Max-Min version.

Depending on the situation, one of these formulations may be more appropriate or even more ethical. For example, when forming student project groups, prioritizing overall happiness at the expense of a single miserable group doesn’t seem fair to me. That’s why, in such cases, the Max-Min approach is likely the better choice. However, in other contexts, the Max-Sum approach might be more suitable. Ultimately, the decision is yours – but now, you have the tools to formulate and solve both using Python.

Quadratic problems like these are challenging, and as the problem size increases, you might notice the solver struggling to find a solution. This is why such problems are often tackled using metaheuristics. I’ve included a link to a state-of-the-art metaheuristic approach for these quadratic assignment problems:

https://www.sciencedirect.com/science/article/abs/pii/S0377221720309814

Additionally, here’s a link to one of my research papers, where I explore a similar variant of the group formation problem using a hybrid method that combines quantum annealing and classical decomposition methods.

Restricted Min-Norm Semi-Assignment in Group Formation: Quantum and Binary Algorithmic Integration

You can find the entire notebook with the code developed for this article in the link to my GitHub repository just below.

math_love_RQMKP_medium/RQMKP_medium_com.ipynb at main · ceche1212/math_love_RQMKP_medium

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] Fry, Hannah. The mathematics of love: Patterns, proofs, and the search for the ultimate equation. Simon and Schuster, 2015.
• [2] Williams HP. Model Building in Mathematical Programming. fourth ed. Chichester: Wiley; 1999
• [3] Bergman, David. "An exact algorithm for the quadratic multiknapsack problem with an application to event seating." INFORMS Journal on Computing 31.3 (2019): 477–492.
• [4] Nasini, Stefano, and Luis Fernando Pérez Armas. "Restricted Min-Norm Semi-Assignment in Group Formation: Quantum and Binary Algorithmic Integration." Available at SSRN 4684889 (2024).
• [5] Silva, Allyson, Leandro C. Coelho, and Maryam Darvish. "Quadratic assignment problem variants: A survey and an effective parallel memetic iterated tabu search." European Journal of Operational Research 292.3 (2021): 1066–1084.
• [6] Cacchiani, Valentina, et al. "Knapsack problems – An overview of recent advances. Part II: Multiple, multidimensional, and quadratic knapsack problems." Computers & Operations Research 143 (2022): 105693.
• https://www.perfecttableplan.com/help/52/windows/html/seat_assignment.htm
• https://towardsdatascience.com/how-many-pok%C3%A9mon-fit-84f812c0387e