3 May 2025
'The Muffin Problem' is a recreational maths puzzle that is simple to state but hard to solve in general. For example:
You have 11 muffins and 6 students. You want to divide the muffins equally between the students, but no student wants a small piece. How should you slice the muffins to maximize the size of the smallest piece?
We're one muffin short of giving each student two whole muffins. An obvious solution is to give each student one whole muffin. Then slice 1/6 from each of the five remaining muffins, giving 5/6 muffin to five students and 5 x 1/6 muffin to the sixth student. Therefore, each student gets 1 + 5/6 muffins. In this solution, the smallest slice is 1/6 muffin. But is that the maximum size we can make the smallest slice?
Spoiler alert: it isn't. Have a go at working out a way to slice the muffins that maximizes the size of the smallest slice while giving each student an equal share of the total. Note that we assume each muffin is circular and they are sliced from the edge to the centre. In this small example, finding a good solution is straightforward. Finding an optimal solution is more difficult. In many cases, proving that a solution is optimal is very difficult.
In this article, we solve The Muffin Problem using a Mixed Integer Linear Program (MILP) for a range of combinations for the number of muffins and the number of students. A MILP is not the most efficient method for solving this problem, but it illustrates a way to approach many types of recreational maths puzzles. Along the way, we explore using fractions in a Python program, speed up the process by solving optimization model instances in parallel, include symmetry-breaking constraints and objective function bounds, and we eat a few muffins – purely for research purposes, of course.
Download the models
The model described in this article is built in Python, using the Pyomo library.
The files are available on GitHub.
Situation
Given a number of students and a number of muffins, we want to divide the muffins equally between the students. Our goal is to maximize the size of the smallest slice of muffin. Each muffin is circular and the slices are cut from the edge to the centre of a muffin.
We want to solve The Muffin Problem for a range of values for the number of muffins and the number of students – say, up to two dozen muffins and two dozen students.
History of The Muffin Problem
The Muffin Problem was proposed by Alan Frank in 2009 as a recreational maths puzzle. Since then, it has been the subject of numerous discussions, presentations, papers, and even, as shown in Figure 1, a 228 page book Mathematical muffin morsels: Nobody wants a small piece.
Solving The Muffin Problem is straightforward for a small number of muffins and students. But as the number of muffins and/or students increases, the puzzle soon becomes hard to solve. Good solutions can be found readily, but proving the optimality of a solution is surprisingly difficult.
In the decade following creation of The Muffin Problem, various people proposed partial solution techniques and optimal solutions for specific cases. For a selection of those solutions and related discussions, see The muffin website by William Gasarch (one of the muffin book's authors).
A decade after the puzzle's creation, in 2019, Richard E. Chatwin published an academic paper describing a complete set of solution techniques for The Muffin Problem: An optimal solution for The Muffin Problem. The paper has 89 pages of detailed theorems, proofs, and algorithms. What started as a simple recreational maths puzzle turned into a complex combinatorial optimization problem.
On GitHub there is a collection of algorithms for solving The Muffin Problem using the Julia programming language: The Muffin Package. The data folder includes solutions for many cases – we use those solutions to verify our model's results. Note that the solutions exclude trivial cases and results that can be derived from other cases.
Example of 11 muffins and 6 students
Simple solution
Returning to the example question we posed at the start of this article, suppose we have 11 muffins and 6 students. How can we divide the muffins equally between the students?
Figure 2 shows a simple solution. We start by giving each student one whole muffin. Then we slice 1/6 from each remaining muffin. The five muffins of size 5/6 go to the first five students. The last student gets the five 1/6 slices. Therefore, each student gets a total of 1 + 5/6 muffins. There are three different sizes (6 x 1, 5 x 5/6, and 5 x 1/6), with the smallest size being 1/6 (0.1667) of a muffin. This is not a good solution, because the students don't like receiving small pieces of muffin.
Optimal solution
In contrast to the simple solution, the optimal solution for 11 muffins and 6 students is more complex, as shown in Figure 3. Each student still gets 1 + 5/6 muffins. but the smallest size is 7/18 (0.3889) of a muffin – more than twice the size of the simple solution. In total, there are six different sizes (1 x 1, 7 x 7/18, 2 x 8/19, 2 x 9/18, 2 x 10/18, and 7 x 11/18).
Model formulation
Since The Muffin Problem is a combinatorial optimization problem, an obvious way to solve it is using a Mixed Integer Linear Program (MILP). Erwin Kalvelagen's article The Muffin Problem includes a MILP formulation, which we replicate. In addition, we include some bounds defined in the literature and a pair of symmetry-breaking constraints intended to make the model easier to solve.
Figure 4 shows the formulation for our model. That is:
- Equation (1). Maximize the size of the smallest slice of muffin.
- Equation (2). Each student gets an equal share of the total muffins.
- Equation (3). The slices of a muffin sum to a whole muffin.
- Equation (4) and (5). Linearized implications: If \(\delta\) = 0, then \(Size\) must be 0; and If \(\delta\) = 1, then \(Smallest\ \leq Size\).
- Equation (6) and (7). Symmetry-breaking constraints. We tried a few other variations, but none worked better than these constraints.
- Equation (8) to (10). Known bounds on the size of the smallest slice.
Implementation
Pyomo model
We implement the model in Python, using the Pyomo library. For a given number of muffins and number of students, the objective function and core constraints (excluding symmetry-breaking and bounds) are shown in Figure 5.
model.objective = Objective(expr=model.smallest, sense=maximize)
model.constraints = ConstraintList()
for s in model.students:
model.constraints.add(sum(model.size[m, s] for m in model.muffins) == num_muffins / num_students)
for m in model.muffins:
model.constraints.add(sum(model.size[m, s] for s in model.students) == 1)
for m in model.muffins:
for s in model.students:
model.constraints.add(model.size[m, s] ≤= model.delta[m, s])
model.constraints.add(model.smallest ≤= model.size[m, s] + (1 - model.delta[m, s]))In addition, Figure 6 shows the symmetry-breaking constraints and some simple bounds defined in the literature.
for s in range(1, num_students):
model.constraints.add(model.size[1, s] <= model.size[1, s + 1]) # Symmetry-breaking for students
for m in range(1, num_muffins):
model.constraints.add(model.size[m, 1] <= model.size[m + 1, 1]) # Symmetry-breaking for muffins
model.bounds = ConstraintList() # Known bounds from the literature
model.bounds.add(model.smallest >= 1/num_students)
if not (num_muffins % num_students == 0):
model.bounds.add(model.smallest <= 0.5)
if num_muffins > num_students:
model.bounds.add(model.smallest >= 1/3)Run cases in parallel using the multiprocessing library
We want to solve many cases, representing all combinations for the number of muffins and the number of students in given ranges. We're using the HiGHS solver, which is mostly single threaded so, to take advantage of the multiple CPU cores in our PC, we run multiple cases in parallel. We use, as a template, the multiprocessing library structure described in our previous article 10 times faster, running cases in parallel.
The main code for running the parallel cases is shown in Figure 7. That is, we create a multiprocessing Manager and some shared dictionaries for collating results. We then create all combinations of 1..num_muffins and 1..num_students in the generate_cases function. Each case is independent, so they are run in parallel in the tasks function. Each task creates a Pyomo model with the given number of muffins and students, solves that model, and prints a one-line result. A summary of collated results is printed after all cases have completed.
time_start = datetime.now()
manager = mp.Manager()
results_dict, status_dict, time_dict = initialize_dicts(manager, UB_MUFFINS, UB_STUDENTS)
cases = generate_cases(UB_MUFFINS, UB_STUDENTS)
lock = manager.Lock()
pool = mp.Pool(processes=PROCESSES_TO_USE) # Create a pool with the user-input number of processes
pool.map(func=setup, iterable=[(UB_MUFFINS, UB_STUDENTS)]) # Do one-off setup
try:
results = pool.map_async(tasks, [(i, case[0], case[1], lock, results_dict, status_dict, time_dict) for i, case in enumerate(cases)])
pool.close()
pool.join()
except KeyboardInterrupt:
pool.terminate()
pool.join()
print_summary(cases, UB_MUFFINS, UB_STUDENTS, results_dict, status_dict, time_start)Model results
One line for each case
Each case prints a one-line result showing, amongst other things, the size of the smallest muffin slice (the Objective and Fraction columns) and the specific sizes that the muffins are cut into. A portion of that output is shown in Figure 8.
In this run of the model we have a time limit of 300 seconds. All the cases have a solution, though some of the larger cases have not proven optimality within that time.
Note that the case numbers are in the order they completed. As is often the way with running cases in parallel, the completion order is mixed up, depending on when each case started, the speed of the CPU core that it is assigned to, and how long that case takes to solve.
Case Muffins Students Time Objective Fraction Gap Status Per student Sizes
---------------------------------------------------------------------------------------------------------------------------------------------------------------
[...]
38 12 4 0.0 1.000000 1 0.00% optimal 3 12x1
66 12 11 19.9 0.363636 4/11 0.00% optimal 1 1/11 6x4/11, 6x5/11, 6x6/11, 6x7/11
45 12 5 300.4 0.400001 2/5 20.00% maxTimeLimit 2 2/5 8x2/5, 8x3/5, 4x1
44 11 5 300.4 0.433333 13/30 3.85% maxTimeLimit 2 1/5 8x13/30, 2x7/15, 2x1/2, 2x8/15, 8x17/30
55 11 7 300.4 0.392857 11/28 23.07% maxTimeLimit 1 4/7 4x11/28, 4x13/28, 6x1/2, 4x15/28, 4x17/28
59 11 8 300.5 0.333333 1/3 37.50% maxTimeLimit 1 3/8 8x1/3, 2x3/8, 2x5/12, 2x11/24, 2x1/2, 2x13/24, 2x7/12, 2x5/8, 2x2/3
60 12 8 0.0 0.500000 1/2 0.00% optimal 1 1/2 10x1/2, 7x1
61 10 9 300.4 0.333333 1/3 33.33% maxTimeLimit 1 1/9 4x1/3, 6x4/9, 6x5/9, 4x2/3
63 12 9 300.4 0.333333 1/3 33.33% maxTimeLimit 1 1/3 7x1/3, 4x1/2, 7x2/3, 3x1
56 12 7 2.2 0.428571 3/7 0.00% optimal 1 5/7 12x3/7, 12x4/7
46 7 6 3.0 0.333333 1/3 0.00% optimal 1 1/6 3x1/3, 2x5/12, 4x1/2, 2x7/12, 3x2/3
50 11 6 300.4 0.388889 7/18 7.14% maxTimeLimit 1 5/6 7x7/18, 3x4/9, 2x1/2, 3x5/9, 7x11/18
54 10 7 300.5 0.333333 1/3 42.86% maxTimeLimit 1 3/7 5x1/3, 1x8/21, 3x3/7, 2x10/21, 2x11/21, 3x4/7, 1x13/21, 2x2/3, 1x1
62 11 9 300.4 0.361111 13/36 12.82% maxTimeLimit 1 2/9 4x13/36, 2x7/18, 4x5/12, 2x1/2, 4x7/12, 2x11/18, 4x23/36
64 11 10 300.4 0.350000 7/20 26.74% maxTimeLimit 1 1/10 4x7/20, 2x2/5, 4x9/20, 2x1/2, 4x11/20, 2x3/5, 4x13/20Summary tables of all cases
After all cases are complete, the model prints a summary of the results. As shown in Figure 9 for up to 12 muffins and 12 students, the optimal solution can always be expressed as a rational fraction. To help us work with the fractions in the output, the model code makes extensive use of Python's fractions library.
Smallest piece for each case
============================
Students
Muffins 1 2 3 4 5 6 7 8 9 10 11 12
--------------------------------------------------------------------------------------------
1 0 0 0 0 0 0 0 0 0 0 0 0
2 1 0 0 0 0 0 0 0 0 0 0 0
3 1 1/2 0 0 0 0 0 0 0 0 0 0
4 1 1 1/3 0 0 0 0 0 0 0 0 0
5 1 1/2 5/12 3/8 0 0 0 0 0 0 0 0
6 1 1 1 1/2 2/5 0 0 0 0 0 0 0
7 1 1/2 5/12 5/12 1/3 1/3 0 0 0 0 0 0
8 1 1 4/9 1 2/5 1/3 5/14 0 0 0 0 0
9 1 1/2 1 7/16 2/5 1/2 5/14 3/8 0 0 0 0
10 1 1 4/9 1/2 1 5/12 1/3 3/8 1/3 0 0 0
11 1 1/2 11/24 9/20 13/30 7/18 11/28 1/3 13/36 7/20 0 0
12 1 1 1 1 2/5 1 3/7 1/2 1/3 2/5 4/11 0
Status for each case
====================
* = optimal, - = not optimal, . = no solution
Students
Muffins 1 2 3 4 5 6 7 8 9 10 11 12
--------------------------------------------------------------------------------------------
1 . . . . . . . . . . . .
2 * . . . . . . . . . . .
3 * * . . . . . . . . . .
4 * * * . . . . . . . . .
5 * * * * . . . . . . . .
6 * * * * * . . . . . . .
7 * * * * * * . . . . . .
8 * * * * * * * . . . . .
9 * * * * * * * * . . . .
10 * * * * * * - * - . . .
11 * * * * - - - - - - . .
12 * * * * - * * * - * * .
Elapsed time: 603.0 seconds, total time 3121.8 secondsMuffin duality theorem
A feature of our model is that it evaluates only the lower triangle of the \(m \times s\) matrix. The solutions for the main diagonal are all obviously 1, as we don't need to slice any muffins if the number of muffins and the number of students is the same.
Less obviously, the literature describes what is called the "duality theorem" for The Muffin Problem. Given that \(f(m,s)\) is defined as the maximum size of the smallest piece of muffin for \(m\) muffins and \(s\) students, the duality theorem states that the solution for \(f(s,m)\) is inversely symmetrical, as shown in Figure 10.
\(f(s,m) = \frac{s}{m} \times f(m,s) \text{ for any }s > m\)
For example, given that \(f(11,6) = \frac{7}{18}\), we can calculate \(f(6,11) = \frac{6}{11} \times \frac{7}{18} = \frac{18}{33} \times \frac{7}{18} = \frac{7}{33} = 0.2121\).
Results for up to 24 muffins and 24 students
The result of running the model for up to 24 muffins and 24 students, then filling in the main diagonal with 1 and upper triangle using the duality theorem, gives the matrix shown in Figure 11. Select or hover over a cell to see the optimal objective function value for that combination of muffins and students.
There are obvious patterns in the matrix, especially for solutions with 1, 1/2, 1/3, 1/4, ... muffins as the smallest piece. But there are also odd solutions, like 53/130 for 23 muffins and 13 students. Using the results in the data folder of The Muffin Package, we've verified that all these solutions are optimal, even though many of the larger cases have not been proven optimal by our model within the 300 second time limit.
Performance of our muffin model
For up to 9 students and 9 muffins, HiGHS solves all 36 cases in 15 seconds. However, as the number of muffins and students gets larger, it becomes harder for our MILP model to prove optimality. The solver typically finds an optimal solution quickly, but within a time limit of 300 seconds some of the larger cases are not proven optimal. We only know the solutions from our model are optimal because we compare them with known optimal solutions.
In some cases, letting the HiGHS solver run for an hour or more leads to no further progress. We also tried some of the unproven cases with the CPLEX and Gurobi solvers – they both find good (actually optimal) solutions quickly and then converge very slowly, failing to prove optimality within several hours.
So, our model finds optimal solutions quickly and then struggles to prove optimality for many cases. The symmetry-breaking constraints and bounds help somewhat. Figure 12 shows the run times and number of proven optimal solutions for solving all 276 cases in the lower triangle of the matrix for up to 24 muffins and 24 students. Running 16 parallel processes on a CPU with 20 cores / 28 threads takes about 1 hour elapsed time, when each case has a maximum time limit of 300 seconds, with no significant variation with or without the bound and symmetry-breaking constraints.
But including the bounds and symmetry-breaking constraints increases the number of cases that are proven optimal within the time limit, with the bounds and symmetry-breaking constraints both contributing. In addition, including the bounds and symmetry-breaking constraints reduces the total CPU time to solve the cases, by 9% from 14.4 hours to 13.1 hours. Running the cases in parallel is 13 or 14 times faster, which is a good improvement over running them serially.
Conclusion
In this article, we use a Mixed Integer Linear Program (MILP) to solve a recreational maths puzzle: The Muffin Problem.
Our goal is to divide a number of muffins equally between a number of students, maximizing the size of the smallest piece that any student receives. Some cases are trivial, while some others are easy. But once we have more than a few muffins and/or students, it becomes difficult for our MILP to prove optimality for many of the cases – even though it does find optimal solutions in all of the cases we looked at.
After a decade of research, efficient methods have been developed to solve any Muffin Problem case – methods that are generally faster than a MILP. But that isn't the point. We just wanted to apply a MILP model to this situation to illustrate how mathematical programming can be a useful tool for solving a variety of recreactional maths puzzles. Along the way, we play with some parallel computing to substantially reduce the run time, and throw in some symmetry-breaking constraints and objective function bounds to help the solver.
Anyway, who needs an excuse to think about, and eat, some muffins.
If you would like to know more about this model, or you want help with your own models, then please contact us.