Integer Programming Models for Round Robin Tournaments

Round robin tournaments are omnipresent in sport competitions and beyond. We propose two new integer programming formulations for scheduling a round robin tournament, one of which we call the matching formulation. We analytically compare their linear relaxations with the linear relaxation of a well-known traditional formulation. We find that the matching formulation is stronger than the other formulations, while its LP relaxation is still being solvable in polynomial time. In addition, we provide an exponentially sized class of valid inequalities for the matching formulation. Complementing our theoretical assessment of the strength of the different formulations, we also experimentally show that the matching formulation is superior on a broad set of instances. Finally, we describe a branch-and-price algorithm for finding round robin tournaments that is based on the matching formulation.


Introduction
Integer programming continues to be a very popular way to obtain a schedule for a round robin tournament. The ability to straightforwardly model such a tournament, and next solve the resulting formulation using an integer programming solver, greatly facilitates practitioners. Moreover, it is usually possible to add all kinds of specific local constraints to the formulation that help addressing particular challenges. We substantiate this claim of the widespread use of integer programming by mentioning some of the works that use integer programming to arrive at a schedule for a round robin tournament. Indeed, from the literature, it is clear that for national football leagues (which are predominantly organized according to a so-called double round robin format), integer programming-based techniques are used extensively to find schedules. Without claiming to be exhaustive we mention Alarcón et al. [2], Della Croce and Oliveri [9], Durán et al [11,12,10], Goossens and Spieksma [17], Rasmussen [27], Recalde et al. [29], Ribeiro and Urrutia [30]. Other sport competitions that are organized in a round robin fashion (or a format close to a round robin) have also received ample attention: we mention Cocchi et al. [8] and Raknes and Pettersen [26] who use integer programming for scheduling volleyball leagues, Fleurent and Ferland [16] who use integer programming for scheduling a hockey league, Kim [21] and Bouzarth et al. [4] for baseball leagues, Kostuk and Willoughby [23] for Canadian football, Nemhauser and Trick [25] and Westphal [36] for basketball leagues. Further, there has been work on studying properties of the traditional formulation, among others, by Trick [33] and Briskorn and Drexl [5]. Wellknown surveys are given by Rasmussen and Trick [28], Kendall et al. [20] and Goossens and Spieksma [18]; we also refer to Knust [22], who maintains an elaborate classification of literature on sports scheduling. More recently, the international timetabling competition [35] featured a round robin sports timetabling problem, and most of the submissions for this competition used integer programming in some way to obtain a good schedule.
All this shows that integer programming is one of the most preferred ways to find schedules for competitions organized via a round robin format.
In this paper, we aim to take a fresh look at the problem of finding an optimal schedule for round robin tournaments using integer programming techniques. Depending upon how often a pair of teams is required to meet, different variations of a round robin tournament arise: in case each pair of teams meets once, the resulting format is called a Single Round Robin, in case each pair of teams is required to meet twice, we refer to the resulting variation as a Double Round Robin. These formats are the ones that occur most in practice; in general we speak of a k-Round Robin to describe the situation where each pair of teams is required to meet k times.
We have organized the paper as follows. In Section 2, we precisely define the problem corresponding to the Single Round Robin tournament, and we present three integer programming formulations for it. We call them the traditional formulation (Section 2.1), the matching formulation (Section 2.2), and the permutation formulation (Section 2.3); the latter two formulations are, to the best of our knowledge, new. We show that their linear relaxations can be solved in polynomial time. We prove in Section 3 that the matching formulation is stronger than the other formulations. In Section 4 we provide a class of valid inequalities for the matching formulation. We show in Section 5 how our results extend to the k-Round Robin tournament. In Section 6, we generate instances of our problem with two goals in mind: (i) to experimentally assess the quality of the bounds found by our models (Section 6.2), and (ii) to report on the performance of a branch-and-price algorithm (Section 6.3). We conclude in Section 7.

Problem definition and formulations
In this section, we provide a formal definition of our problem and introduce the necessary terminology and notation. We start by describing the so-called Single Round Robin (SRR) tournament, where every pair of teams has to meet exactly once, and we return to the general version of the problem, where every pair of teams has to meet k times (k ≥ 1), in Section 5.
Throughout the entire paper, we assume that n is an even integer that denotes the number of teams; for reasons of convenience we assume n ≥ 4. We denote the set of all teams by T . A match is a set consisting of two distinct teams and the set of all matches is denoted by M, in formulae, M = {m = {i, j} : i, j ∈ T, i = j}. We denote, for each i ∈ T , by M i = {{i, j} : j ∈ T \ {i}} the set of matches played by team i. As we assume in this section that every pair of teams meets once, and as n is even, the matches can be organized in n − 1 rounds, which we denote by R; hence, we deal in this section with a compact single round robin tournament.
Prepared with this terminology and notation, we are able to provide a formal definition of the SRR problem. Problem 1. (SRR) Given an even number n ≥ 4 of teams with corresponding matches M, a set of n − 1 rounds R, as well as an integral cost c m,r for every match m ∈ M and round r ∈ R, the single round robin (SRR) problem is to find an assignment A ⊆ M × R of matches to rounds that minimizes the cost (m,r)∈A c m,r such that every team plays a single match per round and each match is played in some round.
Since the SRR problem is NP-hard (see Easton [13], Briskorn et al. [6], and Van Bulck and Goossens [34]), there does not exist a polynomial time algorithm to find an optimal assignment unless P = NP. For this reason, several researchers have investigated integer programming (IP) techniques for finding an optimal assignment of matches to rounds. We follow this line of research and discuss three different IP formulations for the SRR problem: a traditional formulation with polynomially many variables and constraints (Section 2.1) as well as two formulations that involve exponentially many variables (Sections 2.2 and 2.3). To the best of our knowledge, the latter models have not been discussed in the literature before.

The traditional formulation
The traditional formulation of the SRR problem has been discussed, among others, by Trick [33] and Briskorn and Drexl [5]. To model an assignment of matches to rounds, this formulation introduces, for every match m ∈ M and round r ∈ R, a binary decision variable x m,r to model whether match m is played at round r (x m,r = 1) or not (x m,r = 0). With these variables, problem SRR can be modeled as: Constraints (T2) ensure that each pair of teams meets once, and Constraints (T3) imply that each team plays in each round. This model has O(n 2 ) constraints and O(n 3 ) variables. Note that Constraints (T4) can be replaced by x m,r ∈ Z + as the upper bound x m,r ≤ 1 is implicitly imposed via Constraints (T2) and non-negativity of variables. The linear programming relaxation of (T) arises when we replace (T4) by x m,r ≥ 0; given an instance I of SRR, we denote the resulting value by v LP tra (I).

The matching formulation
Consider the complete graph that results when associating a node to each team, say K n = (T, M). Clearly, a single round of a feasible schedule can be seen as a perfect matching in this graph. This observation allows us to build a matching based formulation by introducing a binary variable for every perfect matching in K n ; we denote the set of all perfect matchings in K n by M.
We employ a binary variable y M,r for each perfect matching M ∈ M and round r ∈ R. If y M,r = 1, the model prescribes that matching M is used for the schedule of round r, whereas y M,r = 0 encodes that a different schedule is used. To be able to represent the cost of round r ∈ R, the total cost of all matches in M is denoted by d M,r := m∈M c m,r , which leads to the model Constraints (M2) ensure that a matching is selected in each round, while Constraints (M3) enforce that each pair of teams meets in some round. Similarly to the traditional formulation, we can replace (M4) by y M,r ∈ Z + . In this way, the linear programming relaxation of (M) arises when replacing (M4) by y M,r ≥ 0; given an instance I of SRR, the resulting value is denoted by v LP mat (I). Notice that this formulation uses an exponential number of variables, as the number of matchings grows exponentially in n. Thus, a relevant question is whether we can find v LP mat in polynomial time. The following observation shows that it can be answered affirmatively. Proof. Due to the celebrated result by Grötschel et al. [19], it is sufficient to show that the separation problem for the constraints of the dual of the linear relaxation of Model (M) can be solved in polynomial time. To avoid an exponential number of variables in the dual, we replace Constraint (M4) by y M,r ∈ Z + as explained above. Then, by introducing dual variables α r , r ∈ R, corresponding to Constraints (M2) and β m , m ∈ M, corresponding to Constraints (M3), the constraints of the dual of the LP relaxation of Model (M) are: Given values for the dual variables, say (ᾱ,β), the separation problem is to decide whether it satisfies all dual constraints. For fixed r ∈ R, we show that this problem can be solved in polynomial time. Thus, the assertion follows as there are only O(n) rounds. Indeed, if r ∈ R is fixed, the problem reduces to check whether there exists a matching M ∈ M such thatᾱ The latter inequality asks whether there exists a perfect matching of teams with weight greater than −ᾱ r , where an edge m between two teams is assigned weight (β m − c m,r ). This problem can be solved in polynomial time by Edmonds' blossom algorithm [14,15], which concludes the proof.

The permutation formulation
Instead of fixing the schedule of a round, the permutation formulation fixes, for a given team, the order of the teams against which the given team plays its successive matches. That is, it introduces a variable for each team i and each permutation of T \ {i}. We denote the set of all such permutations by Π −i . Moreover, for a team j ∈ T and round r ∈ R, denote the set of all permutations where j occurs at position r in the permutation by Π −i j,r . Permutations from Π −i j,r thus encode that team i plays against team j on round r. For a permutation π ∈ Π −i and round r ∈ R, we refer to the opponent of team i at round r as π r ∈ T \ {i}. The cost of a schedule encoded via permutations Π −i for a team i ∈ T is then given by e i,π := r∈R c {i,πr},r . Using binary variables z i,π , where i ∈ T and π ∈ Π −i , that encode whether i plays against its opponents in order π (z i,π = 1) or not (z i,π = 0), the permutation formulation is Constraints (P2) ensure that a permutation is selected for each team, while Constraints (P3) enforce that, given a round and a pair of teams, these teams meet in that round, or they do not meet in that round. Due to rescaling the objective by 1 2 , we find the cost of an optimal SRR schedule. Moreover, we can again replace Constraint (P4) by z i,π ∈ Z + . The linear programming relaxation of (P) then arises when replacing Constraints (P4) by z i,π ≥ 0; given an instance I of SRR, we denote the resulting value by v LP per (I). Since this model has n! variables, we again investigate whether its LP relaxation can be solved efficiently. Proof. As in the proof of Lemma 2.1, it is sufficient to show that the separation problem corresponding to the constraints of the dual of the relaxation of the permutation formulation can be solved in polynomial time. Again, to avoid exponentially many variables in the dual, we replace (P4) by z i,π ∈ Z + . We introduce dual variables α i for each constraint of type (P2) and β {i,j},r for each constraint of type (P3). To normalize Constraint (P3), we assume it to be given by π∈Π −i j,r z i,π − π∈Π −j i,r z j,π = 0 with i < j. Then, the dual constraints are given by If i ∈ T is fixed, the separation problem for dual values (ᾱ,β) is to decide whether there exists a permutation π ∈ Π −i such that To answer this question, it is sufficient to find a permutation maximizing the left-hand side expression. Such a permutation can be found by computing a maximum weight perfect matching in the complete bipartite graph with node bipartition (T \ {i}) ∪ R and edge weights defined for each j ∈ T \ {i} and r ∈ R by Since this problem can be solved in polynomial time, the assertion follows by solving this problem for each of the n teams.

Comparing the strength of the different formulations
In the previous section, we have introduced three different models for finding an optimal schedule for problem SRR. While the traditional formulation contains both polynomially many variables and constraints, the matching and permutation formulation make use of an exponential number of variables. The aim of this section is to investigate whether the increase in the number of variables in comparison with the traditional formulation leads to a stronger formulation. We measure the strength of a formulation based on the value of its LP relaxation, where a higher value of the LP relaxation indicates a stronger formulation as the LP relaxation's value is closer to the optimum value of the integer program, as encapsulated by the following definitions.
Definition 3.1. Let f and g be mixed-integer programming formulations of the SRR problem and denote by v LP f (I) and v LP g (I) the value of the respective LP relaxations for an instance I of SRR. • We say that f and g are relaxation-equivalent if, for each instance I of problem SRR, the value of the linear programming relaxations are equal, i.e., v LP f (I) = v LP g (I).
• We say that f is stronger than . We now proceed by formally comparing the strength of the formulations from Section 2 using the terminology of these definitions. We state our results using three lemmata, and summarize all our results in Theorem 3.5.
First, we show that the traditional and permutation formulation have equivalent LP relaxations.

Lemma 3.2. The permutation formulation (P) is relaxation-equivalent to the traditional formulation (T).
Proof. To prove this lemma, we show that there is a one-to-one correspondence of feasible solutions of the traditional formulation's and the permutation formulation's LP relaxations that preserves the objective value. First, we construct a solution of the LP relaxation of the traditional formulation from a solution z of the LP relaxation of the permutation formulation. To this end, define for each {i, j} ∈ M and r ∈ R a solution x ∈ R M×R via x {i,j},r = Note that x is non-negative as all z-variables are non-negative. Moreover, it is well-defined as π∈Π −i j,r z i,π = π∈Π −j i,r z j,π due to (P3). Finally, all constraints of type (T2) and (T3) are satisfied since for each m ∈ M : for each i ∈ T, r ∈ R : j∈T \{i} We conclude the proof by constructing a feasible solution for the LP relaxation of the permutation formulation from a feasible solution x of the traditional formulation's LP relaxation. Let x be such a solution and let i ∈ T . Consider the matrix X i ∈ R (T \{i})×R with entries X i j,r = x {i,j},r . Due to all constraints of the traditional formulation's LP relaxation, X i is a doubly stochastic matrix and is thus contained in the Birkhoff polytope, see [37]. Consequently, X i can be written as a convex combination of all permutation matrices. That is, if P i,π is the permutation matrix associated with π ∈ Π −i , there exist multipliers λ i π ≥ 0, π ∈ Π −i , such that X i = π∈Π −i λ i π P i,π and π∈Π −i λ i π = 1. Based on these multipliers, we define a solution z of the permutation formulation via z i,π = λ i π . To conclude the proof, we need to show that this solution z is feasible for the permutation formulation's LP relaxation and has the same objective value as x. Observe that z is non-negative since all λ's are non-negative. Constraints (P2) and (P3) are satisfied as for each i ∈ T : since π∈Π −i λ i π = 1 and x {i,j},r is a convex combination of permutation matrices that assign team j (or i) to round r, respectively. Consequently, z is feasible for the permutation formulation's LP relaxation. Finally, both x and z have the same objective value because which proves that both formulations are relaxation-equivalent.
Next, we turn our focus to the matching formulation and compare it with the traditional formulation (and thus, by the previous lemma, also with the permutation formulation).

Lemma 3.3. For each n ≥ 6, the matching formulation (M) is stronger than the traditional formulation (T).
Proof. First, we show that we can transform any feasible solution of the matching formulation's LP relaxation to a feasible solution of the traditional formulation's LP relaxations. Afterwards, to show that the matching formulation is stronger than the traditional formulation, we show that, for any even n ≥ 6, there exists an instance of SRR for which the LP relaxation of the matching formulation has a strictly larger value than the traditional formulation's LP relaxation.
Let y be a feasible solution of the matching formulation's LP relaxation. We construct a solution x for the traditional formulation by setting x m,r = M ∈M : m∈M y m,r . Since y is non-negative, also x is non-negative. Moreover, Conditions (T2)  that is, the traditional formulation cannot be stronger than the matching formulation.
To prove that the matching formulation dominates the traditional formulation for n ≥ 6 even, we distinguish three cases. In the first case, assume n ≥ 10. Consider the pairs of teams given by Interpreting P as the edges of an undirected graph, P defines three connected components consisting of two 3-cycles and an even cycle. We construct an instance of the SRR problem by specifying the cost function c ∈ R M×R via It is easy to verify that x ∈ R M×R given by if m / ∈ P and r ∈ {1, 2}, 1 n−3 , otherwise, is feasible for the LP relaxation of the traditional formulation and has objective value 0. Hence, x is optimal.
Solving the LP relaxation of the matching formulation for this instance, however, results in an objective value that is at least 2. Indeed, each perfect matching M ∈ M contains at least one match m ∈ M with m ∈ {{i, j} : (i, j) ∈ {1, 2, 3} × {4, 5, . . . , n}}. Since c m,1 = c m,2 = 1 for such a match, it follows that in both rounds 1 and 2, matchings are selected with total weight at least 1 due to (M3),leading to a solution with total cost at least 2.
In the second case, we consider n = 6. To prove the statement, we use the same construction as before, however, we do not require the even cycle anymore. That is, P defines two 3-cycles and the argumentation remains the same as before.
In the last case n = 8, we consider the set of pairs If we interpret P as edges of an undirected graph, the corresponding graph has two connected components being a 3-cycle and a 5-cycle, respectively. We choose the cost-coefficients c ∈ R M×R to be Simple calculations show that an optimal solution of the traditional formulation's LP relaxation is x ∈ R M×R with ∈ E and r ∈ {1, 2}, 1 5 , otherwise, which has objective value 0, whereas the matching formulation's LP relaxation has value 2.
The previous two lemmata completely characterize the relative strength of the three different formulations except for n = 4. The status of this missing case is settled next. Proof. Observe that exactly the same arguments as in the proof of Lemma 3.3 can be used to show that the matching formulation is at least as strong as the traditional formulation if n = 4. Hence, it remains to show that every solution x of the traditional formulation's LP relaxation can be turned into a solution of the LP relaxation of the matching formulation if n = 4. Let x be such a solution. Let M = {i, j}, {k, l} ∈ M. Since x satisfies Equations (T3), summing the equations for i, j and subtracting the equations for k, l yields 2x {i,j},r − 2x {k,l},r = 0. That is, the x-variables for the two matches of a matching within the same round have the same value. Consequently, the solution y ∈ R M×R given by y M,r = x {i,j},r is well-defined and it is immediate to check that y is feasible for the matching formulation's LP relaxation and has the same objective value as x.
Summarizing the previous results of this section, we can provide a complete comparison of the strength of the traditional, matching, and permutation formulation. Besides verifying that all three models are equivalent for n = 4, we can also show that the matching formulation's integer hull is already completely characterized by (M2), (M3), as well as non-negativity inequalities for all variables. The non-trivial constraints from Formulation (M) are Equations (M2) and (M3), which yield system Note that the last three equations are redundant and can be removed. The constraint matrix of the remaining equations is the node-edge incidence matrix of a bipartite graph and hence totally unimodular, which concludes the proof.
Thus, for n = 4, simply solving the LP-relaxation of the matching formulation by the simplex method, suffices to find an optimum integral solution.

Strengthening the formulations
In this section, we continue our investigations of the structure of the formulations. In Section 4.1, we derive an exponentially sized class of valid inequalities for the matching formulation. Also, we show in Section 4.2 that adding the so-called odd-cut inequalities to the traditional formulation yields a formulation that is relaxation-equivalent to the matching formulation.

Strengthening the matching formulation
Observe that Theorem 3.5 does not rule out the possibility that, for n ≥ 6, every vertex of the matching formulation is integral. That, however, is not the case already for n = 6 as we will show next. To this end, we first provide a fractional point y that is contained in the LP relaxation of the matching formulation for n = 6. Afterwards, we derive a class of valid inequalities for the integer hull matching formulation, and finally, we provide one such inequality that is violated by y .  Figure 1, we depict a fractional solution of the matching formulation's LP relaxation. For each round r ∈ R, we provide two perfect matchings between the teams T , the blue and green (dashed) matching M , whose corresponding variables y M,r have value 1 2 in the corresponding solution; all remaining variables have value 0. It is easy to verify that this fractional solution is indeed feasible for the LP relaxation of (M). To describe our class of valid inequalities, consider the following lemma.
is a valid inequality for (M). In particular, it is a Chvátal-Gomory cut derived from the LP relaxation of (M).
Proof. It is sufficient to prove that (4)

2.
Note that Inequality (4) is a so-called {0, 1 2 }-cut [7] as all multipliers used in the derivation are 1 2 (and 0 for inequalities/equations that have not been used). By taking more equations in the generation of a valid inequality into account, we can generalize (4) to an exponentially large class of inequalities.
is a valid inequality for (M). In particular, it is a Chvátal-Gomory cut derived from the LP relaxation of (M). Since all y-variables are non-negative, we derive the inequality and by integrality of the y-variables, we can round up the right-hand side, which leads to the desired inequality.
While Inequalities (4) can trivially be separated in polynomial time, an efficient separation algorithm for (5) is not immediate. We leave the complexity status of separating (5) open for future research.

Strengthening the traditional formulation
Revisiting the proof of Lemma 3.3, it becomes clear that it is possible to assign, for a fixed round, each edge (match) of an odd cycle in K n a weight of 1 2 . That is, the traditional formulation can assign an odd cycle of length k a weight of k 2 . Such a solution, however, cannot be written as a convex combination of integer feasible solutions, because each such solution defines a perfect matching on the matches of a fixed round, i.e., the total weight of an odd cycle can be at most k−1 2 . To strengthen the traditional formulation, one can thus add facet defining inequalities for the perfect matching polytope P M to Model (T), which results in the additional inequalities i∈U j∈T \U which correspond to the odd-cut inequalities for the matching polytope and can be separated in polynomial time. Proof. We use the same proof strategy as for Lemma 3.3. Therefore, consider again the solution x ∈ R M×R given by x m,r = M ∈M : m∈M y M,r for a solution y of the matching formulation's LP relaxation. Due to the proof of Lemma 3.3, it is sufficient to show that x satisfies (6) to prove that the matching formulation is at least as strong as the enhanced traditional formulation.
Let U ⊆ T have odd cardinality. Since every M ∈ M is a perfect matching, there is at least one team i ∈ U that does not play against another team in U since U is odd. Hence, for each M ∈ M, there is a match {i, j} ∈ M with i ∈ U and j / ∈ U . Then, i∈U j∈T \U Consequently, the matching formulation is at least as strong as the enhanced traditional formulation.
To prove that the enhanced traditional formulation is not weaker than the matching formulation, we use a strategy similar to the one pursued in the proof of Lemma 3.2. Since the enhanced traditional formulation contains, per round r, all facet defining inequalities as well as equations for the perfect matching polytope P M , each vector X r ∈ R M given by X r {i,j} = x {i,j},r is contained in P M . Hence, there exist non-negative multipliers λ which concludes the proof.

Remark 4.5.
Since the traditional and permutation formulation are equivalent, one might wonder whether also the permutation formulation can be enhanced by odd-cut inequalities. Indeed, using the transformation x {i,j},r = π∈Π −i j,r z i,π as in the proof of Lemma 3.2, one can show that the corresponding version of odd-cut inequalities is given by and that the enhanced traditional and permutation formulation are equivalent.

An extension: k-round robin tournaments
In this section, we generalize the models for single round robin tournaments to k-round robin tournaments, where each pair of teams is required to meet exactly k times, for k ≥ 1. As a consequence, the total number of matches that need to be scheduled becomes 1 2 kn(n − 1), and we set R := {1, 2, . . . , k(n − 1)}. Problem 2 (kRR). Let n ≥ 4 be an even integer and let k ≥ 1 be integral. Given n teams with corresponding matches M, a set of 1 2 kn(n − 1) rounds R (k ≥ 1), as well as an integral cost c m,r for every m ∈ M and round r ∈ R, the k-round robin (kRR) problem is to find an assignment A ⊆ M × R of matches to rounds such that (i) every team plays a single match per round, (ii) each match is played in k, pairwise distinct, rounds, while total cost (m,r)∈A c m,r is minimized.
Problem 1 (SRR) is a special case of kRR as it arises when k = 1. Another very prominent special case arises when k = 2, the so-called Double Round Robin tournament, denoted hereafter by DRR.
In principle, it is easy to generalize the models from Section 2 to account for meeting k times instead of once. Indeed, by replacing the right-hand side of constraints (T2), or the right-hand side of constraints (M3) by k, or by redefining Π −i and Π −i j,r to ordered lists for team i that features every opponent j exactly k times, the resulting formulations for kRR directly arise. In fact, we claim that it is not difficult to verify that the results concerning the polynomial solvability of the linear relaxations (Lemmata 2.1 and 2.2), as well as the strength of the relaxations (Theorem 3.5) hold for the kRR for each k ≥ 1.
However, in practice, a number of additional properties become relevant when considering k-round robin tournaments: phased tournaments and tournaments where playing home or away matters. We now discuss these properties, and their consequences for the formulations, in more detail.
Phased (PH) The tournament is split into k parts such that each pair of teams meets once in each part. Here a part of the tournament refers to n − 1 consecutive rounds, starting at round (n − 1) + 1, for ∈ {0, . . . , k − 1}. Moreover, we use R := { (n − 1) + 1, . . . , ( + 1)(n − 1)} to denote the rounds in part ∈ {0, . . . , k − 1}, and R := k−1 =0 R . Without the presence of any additional constraints, a phased tournament can be trivially decomposed in multiple single-round robin tournaments: one for each set of rounds R .
Home-away (HA) Each team has a home venue, implying that to specify a schedule it is no longer sufficient to specify the matches in each round; instead, one also has to specify, for each match, which teams plays home, and which team plays away. We denote this by redefining a match between teams i, j ∈ T (i = j) where i is the home-playing team, by an ordered pair (i, j) (in contrast to an unordered pair {i, j}).
Let k be a positive integer. For n teams and a k-round robin setting, denote T := {1, . . . , n} and R := {1, . . . , k(n − 1)}. Let M := {(i, j) : i, j ∈ T, i = j} be the set of ordered matches. The assignment of match (i, j) ∈ M to round r ∈ R comes at a cost c (i,j),r , and in contrast to the SRR case, c (i,j),r and c (j,i),r can be different. We proceed by describing the phased k-round robin problem with home-away patterns (kRR-PH-HA): Problem 3 (kRR-PH-HA). Given an even number n ≥ 4 of teams with corresponding matches M, a set of 1 2 kn(n − 1) rounds R (k ≥ 1), as well as an integral cost c m,r for every m ∈ M and round r ∈ R, the kRR-PH-HA problem is to find an assignment A ⊆ M × R of matches to rounds such that (i) every team plays a single match per round, (ii) each pair of teams meets once in part R , ∈ {1, . . . , k}, and (iii) each ordered match is played k 2 or k 2 times so that each pair of teams meets in total k times, while total cost (m,r)∈A c m,r is minimized.
We will show how the three formulations of Sections 2.1-2.3 can be adapted to deal with these properties.
Extending the traditional formulation for k-RR tournaments For reasons of convenience, we assume k is even; this implies that for each pair of distinct teams i, j ∈ T match (i, j) and Finally, we discuss our experience with solving instances of the SRR problem using the matching formulation, i.e., not only solving the LP relaxation but also the corresponding integer program. To this end, we have implemented a branch-and-price algorithm whose details are described in Section 6.3.

Test set
We have generated 1000 instances 1 of the SRR problem to evaluate the quality of the LP relaxations of the different models. Our test set comprises of instances of different sizes and thus different levels of difficulty, which are parameterized by a tuple (n, ρ) and have cost coefficients attaining values 0 or 1. Parameter n encodes the number of teams and has range n ∈ {6, 12, 18, 24}; parameter ρ controls the number of 1-entries in the objective. More precisely, we pick a set of match-round pairs M × R of size ρ · |M × R| uniformly at random, denoted by S ⊆ M × R, where ρ ∈ {0.5, 0.6, 0.7, 0.8, 0.9}. The generated instance consists of n teams and has cost coefficients c m,r = 1 if (m, r) ∈ S and c m,r = 0 otherwise. For each combination of n ∈ {6, 12, 18, 24} and ρ ∈ {0.5, 0.6, 0.7, 0.8, 0.9} we have generated 50 instances.

A computational comparison of the linear relaxations
In this section, we provide a computational comparison between the LP relaxation values of the traditional formulation (T) (and thus by Lemma 3.2 also of the permutation formulation) and LP relaxation values of the matching formulation (M), and compare these to the actual optimal (or best found) integral solutions. Before we discuss our numerical results, we provide details about our implementation as well as on how we find optimal integral solutions first.

Implementation details
To find the LP relaxation values, we implemented both formulations in Python 3 using the PySCIPOpt 4.1.0 package [24] for SCIP 8.0.0 [3], with CPLEX 20.1.0.0 as LP solver. The traditional formulation is implemented as a compact model. For the matching formulation, we use a column generation procedure that receives a subset of all variables, solves the corresponding LP relaxation restricted to these variables, and adds further variables until it can prove that an optimal LP solution has been found. To identify whether new variables need to be added, we solve the so-called pricing problem, which corresponds to separating a corresponding solution of the dual problem. The separation problem can be solved by finding a maximum weight perfect matching as detailed in the proof of Lemma 2.1. We start with the empty set of variables, which means that the primal problem is initially infeasible. Analogously to the proof of Lemma 2.1, we resolve infeasibility by adding variables to the problem that are associated with a dual constraint that violate a dual unbounded ray of this infeasible problem. The column generation procedure has been embedded in a so-called pricer plug-in of SCIP, which adds newly generated variables to the matching formulation. The maximal weight perfect matchings are computed using NetworkX 2.5.1, which provides an implementation of Edmonds' blossom algorithm.
Finding optimal integer solutions To obtain the optimal integer solution value of as many instances as possible, we have used two different solvers to solve the integer program of Model (T). On the one hand, we have used SCIP as described in the above setup. On the other hand, we have modeled (T) using Gurobi 9.1.2 via its Python 3 interface. For each instance and solver, we have imposed a time limit of 48 h to find an optimal integer solution. Using SCIP, we managed to solve 852 of the 1000 instances to optimality. With Gurobi, we were able to solve 866 of the 1000 instances to optimality. There were 45 instances where SCIP found a better primal objective value, and 79 instances where Gurobi found a better primal objective value. All experiments have been run on a compute cluster with identical machines, using one (resp. two) thread(s) on Xeon Platinum 8260 processors, with 10.7 GB (resp. 21.4 GB) memory, respectively for SCIP and Gurobi. Numerical results Table 1 shows the aggregated computational results of our experiments. For each number of teams n and ratio ρ, we provide the average of the objective values of the relaxation of the traditional formulation (column "v LP tra "), the average of the objective values of the relaxation of the matching formulation (column "v LP mat "), and the average optimum value (column "v IP "). Notice that for n = 24 we have not been able to solve all instances to optimality; in this case, we use the value of the best known solution instead of the (unknown) optimum for that instance in the v IP column. Recall that each value is an average over 50 instances. The number of optimally solved instances (resp. instances not terminating within the time limit) are shown in column "O" (resp. "T").
To be able to assess the strength of the matching formulation compared to the traditional formulation, we focus, in the right side of the table, on those instances for which v LP tra < v IP ; their number (out of 50) is given in the column labeled "#". From this column, we see that the fraction ρ that leads to instances with a gap between v LP tra and v IP slowly increases with n. Indeed, for n = 6, most instances do not have a gap, for n = 12, almost all instances with ρ ∈ {0.6, 0.7, 0.8} have a gap, and for n = 18, almost all instances with ρ ∈ {0.7, 0.8, 0.9} have a gap.
We use the notion of the relative gap that is closed by the matching formulation relative to the traditional formulation, given by A value of zero for rgap(I) implies that the relaxation values of the traditional formulation and the matching formulation are equal, while a value of one (i.e., 100%) implies that the relaxation of the matching formulation is equal to the true objective of the optimal integral solution. The column "average" gives the average rgap, whereas column "maximal" shows the maximum relative closed gap for an instance of this sub test set. For n = 6, there are few instances with a gap. However, for those instances for which there is a gap, it is clear that a sizable part of that gap is closed by the relaxation of the matching formulation. For larger values of n, many instances have a gap. We observe that a significant percentage of the gap is closed by the relaxation of the matching formulation. If n is getting larger, however, both the value of the average gap closed as well as the value of the maximal gap closed decrease. We conclude that for small values of n, and thus for many realistic applications, the matching formulation provides a much better relaxation value than the traditional formulation.

A branch-and-price algorithm
Since the matching formulation can dominate the traditional formulation, a natural question is whether the stronger formulation also allows to solve the SRR problem faster than the traditional formulation. For this reason, we have implemented a branch-and-price algorithm (in the computational setup as described above) to compute optimal integral solutions of the matching formulation. That is, we use a branch-and-bound algorithm to solve the matching formulation, where each LP relaxation is solved using a column generation procedure.

Implementation details
In classical branch-and-bound algorithms, the most common way to implement the branching scheme is to select a variable x i whose value x i in the current LP solution is non-integral and to generate two subproblems by additionally enforcing either x i ≤ x i or x i ≥ x . In principle, this strategy is also feasible for the matching formulation, where the subproblems correspond to forbidding a schedule M ∈ M for a round r ∈ R or fixing the schedule in round r to be M . This branching scheme, however, leads to a very unbalanced branch-andbound tree as the former subproblem only rules out a very specific schedule, while the latter one fixes the matches of an entire round. Another difficulty of the classical scheme is that it might affect the structure of the pricing problem in the newly generated subproblems. Ideally, the pricing problem should not change such that the same algorithm can be used for adding new variables to the problem. We will address both issues next.
To obtain a more balanced branch-and-bound tree, we have implemented a custom branching rule following the Ryan-Foster branching scheme [31]: Our scheme selects a match {i, j} ∈ M at a round r ∈ R and creates two children. In the left child, we forbid that {i, j} is played in round r, and in the right child, we enforce that {i, j} is played in round r. Note that for all matchings M ∈ M this branching decision fixes all variables y M,r to zero if {i, j} ∈ M for the left child, and {i, j} / ∈ M for the right child. Using this branching strategy, the structure of the pricing problem at each subproblem remains a matching problem. At the root node of the branch-and-bound tree, we need to solve a maximum weight perfect matching problem in a weighted version of K n as described above. At other nodes of the branch-and-bound tree, we have added branching decisions that enforce that two teams i and j either do meet or do not meet in a round r ∈ R. These decisions can easily be incorporated by deleting edges from K n . When generating variables for round r, we remove edge {i, j} from K n if i and j shall not meet in this round; if the match {i, j} shall take place, then we remove all edges incident with i and j except for {i, j}. Consequently, our branching strategy allows to solve the LP relaxations of all subproblems in polynomial time.
Since our Python implementation of the traditional and matching formulation took too much time to be used in a branching scheme, we decided to implement our branch-and-price algorithm as a plug-in using the C-API of SCIP. The pricer plug-in is analogous, and maximal weight perfect matchings are now computed using the LEMON 1.3.1 graph library. To ensure that the branching decisions are taken into account, we also implemented a constraint that fixes y M,r to zero if the matching M violates the branching decisions for round r, and added a plug-in that implements the branching decisions.
The branching rule sketched above admits some degrees of freedom in selecting the match {i, j} and round r. In our implementation, we decided to mimic two well-known branching rules: most infeasible branching and strong branching on a selection of variables, see Achterberg et al. [1] for an overview on branching rules. Most infeasible branching branches on a binary variable with fractional value in an LP solution that is closest to 0.5, and strong branching branches on the variable that yields the largest dual bound improvement based on some metric. Since strong branching requires significant computational effort, it is common to make a limited branching candidate selection and apply strong branching on those.
The pseudocode for our branching rule is given in Algorithm 1. We start by computing the fractional match-on-round assignment values induced by the y-variables. Then, we make a selection of potentially good branching candidates (m, r) ∈ M × R, that is based on the score m,r metric shown in Line 8: this score prioritizes match-on-round assignments for which (i) the assignment value is close to 0.5, (ii) the cost coefficients are large, and (iii) the assignment values are relatively high. Using this score, we hope to resolve fractionality soon (by (i)). By (ii), we want to enforce a significant change of the objective value in the child that forbids match m, whereas the child enforcing that m is played selects a match that is most likely played due to (iii). Experiments show that full strong branching leads to a smaller number of nodes to solve the problem, but this turns out to be very costly computationally. Therefore, we only apply strong branching for branch-andbound tree nodes close to the root, and only evaluate a subset of branching candidates that have the highest score m,r -metric. This is our candidate pre-selection. The higher the depth of the considered branch-and-bound tree node, the smaller the number of candidates considered. Of those branching candidates, we pick the candidate that maximizes score m,r as defined in Line 18. The goal of this score is to choose the candidate where the objective values of the hypothetical children are different from the current node's objective. By considering the product of this difference, we prioritize if the objective of both hypothetical children have some difference with the current node objective. If the number of candidates is only one, no strong branching is applied and that candidate match-on-round assignment is chosen for branching. All experiments have been run on a Linux cluster with Intel Xeon E5 3.5 GHz quad core processors and 32 GB memory. The code was executed using a single thread and the time limit for all computations was 2 h per instance. Table 2 summarizes our results for our instances for n ∈ {6, 12, 18}. We distinguish the instances by their parameters n and ρ, and we report on the number of instances that could be solved (resp. could not be solved) within the time limit in column "O" (resp. "T"). Moreover, we report on the the minimum, mean, and maximum running time per parameterization. The mean of all running times t i is reported in shifted geometric mean 50 i=1 (t i + s) 1 50 − s using a shift of 10 s to reduce the impact of instances with very small running times.

Numerical results
We observe that instances with 6 and 12 teams can be solved very efficiently within fractions of seconds in the former and within seconds in the latter case. Instances with 18 teams are more challenging, in particular, if the ratio ρ ∈ {0.7, 0.8}. In this case, only 3 and 22 instances could be solved, respectively, but note that not all instances are equally difficult. For instance, for n = 18 and ρ = 0.8, there exists an instance that can be solved within roughly five minutes, whereas the mean running time is more than an hour. To fully benefit from the strong LP relaxation of the matching formulation, it might be the case that additional algorithmic enhancements can further improve the performance of the branch-and-price algorithm.

Conclusion
The use of integer programming for finding schedules of round robin tournaments is widespread. We have introduced and analyzed two new formulations for this problem, one of which (the matching formulation) is stronger than the other formulations. We have proposed a class of valid inequalities for the matching formulation, which may be of use when developing cutting-plane based techniques for this problem. By randomly generating instances, we studied the strength of the formulations, and we implemented a branch-and-price algorithm based on the matching formulation to see its efficiency. Although this algorithm is able to solve small-scale instances rather efficiently, solving large instances of the SRR efficiently remains a challenge.
Possible directions of future research are thus to further strengthen our integer programming formulations/techniques. On the one hand, once can investigate additional cutting planes to strengthen both the traditional and matching formulation. For the matching formulation, cutting planes will in particular affect the pricing problem and thus might change its structure. Thus, the trade-off between the strength of cutting planes and the difficulty of solving the pricing problem that needs to be investigated. On the other hand, one can enhance our branch-and-price algorithm in several directions, e.g., the development of more sophisticated branching rules or heuristics for producing good schedules.