Suppose n teams play each other, and let Team Team Team denote some fixed ranking (where is some permutation of ). An upset occurs when a lower ranked team beats an upper ranked team. For each ranking, , let denote the total number of upsets. The minimum upset problem is to find an “efficient” construction of a ranking for which is as small as possible.

In general, let denote the number of times Team i beat team $j$ minus the number of times Team j beat Team i. We regard this matrix as the signed adjacency matrix of a digraph . Our goal is to find a Hamiltonian (undirected) path through the vertices of which goes the “wrong way” on as few edges as possible.

1. Construct the list of spanning trees of (regarded as an undirected graph).
2. Construct the sublist of Hamiltonian paths (from the spanning trees of maximum degree 2).
3. For each Hamiltonian path, compute the associated upset number: the total number of edges transversal in going the “right way” minus the total number going the “wrong way.”
4. Locate a Hamiltonian for which this upset number is as large as possible.

Use this sagemath/python code to compute such a Hamiltonian path.

def hamiltonian_paths(Gamma, signed_adjacency_matrix = []):
    """
    Returns a list of hamiltonian paths (spanning trees of 
    max degree <=2).

    EXAMPLES:
        sage: Gamma = graphs.GridGraph([3,3])
        sage: HP = hamiltonian_paths(Gamma)
        sage: len(HP)
        20
        sage: A = matrix(QQ,[
        [0 , -1 , 1  , -1 , -1 , -1 ],
        [1,   0 ,  -1,  1,  1,   -1  ],
        [-1 , 1 ,  0 ,  1 , 1  , -1  ],
        [1 , -1 , -1,  0 ,  -1 , -1  ],
        [1 , - 1 , - 1 , 1 , 0 , - 1  ],
        [1 ,  1  ,  1  , 1  , 1  , 0 ]
        ])
        sage: Gamma = Graph(A, format='weighted_adjacency_matrix')
        sage: HP = hamiltonian_paths(Gamma, signed_adjacency_matrix = A)
        sage: L = [sum(x[2]) for x in HP]; max(L)
        5
        sage: L.index(5)
        21
        sage: HP[21]                                 
        [Graph on 6 vertices,
         [0, 5, 2, 1, 3, 4],
         [-1, 1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1]]
        sage: L.count(5)
        1

    """
    ST = Gamma.spanning_trees()
    if signed_adjacency_matrix == []:
        HP = []
        for X in ST:
            L = X.degree_sequence()
            if max(L)<=2:
                #print L,ST.index(X), max(L)
                HP.append(X)
        return HP
    if signed_adjacency_matrix != []:
        A = signed_adjacency_matrix
        HP = []
        for X in ST:
            L = X.degree_sequence()
            if max(L)<=2:
                #VX = X.vertices()
                EX = X.edges()
		if EX[0][1] != EX[-1][1]:
                    ranking = X.shortest_path(EX[0][0],EX[-1][1])
		else:
		    ranking = X.shortest_path(EX[0][0],EX[-1][0])
		signature = [A[ranking[i]][ranking[j]] for i in range(len(ranking)-1) for j in range(i+1,len(ranking))]
                HP.append([X,ranking,signature])
        return HP

Wessell describes this method in a different way.

1. Construct a matrix, , with rows and columns indexed by the teams in some fixed order. The entry in the i-th row and the j-th column is defined by
2. Reorder the rows (and corresponding columns) to in a basic win-loss order: the teams that won the most games go at the
   top of , and those that lost the most at the bottom.
3. Randomly swap rows and their associated columns, each time checking if the
   number of upsets has gone down or not from the previous time. If it has gone down, we keep
   the swap that just happened, if not we switch the two rows and columns back and try again.

An implementaiton of this in Sagemath/python code is:

def minimum_upset_random(M,N=10):
    """
    EXAMPLES:
        sage: M = matrix(QQ,[
        [0 , 0 , 1  , 0 , 0 , 0 ],
        [1,   0 ,  0,  1,  1,   0  ],
        [0 , 1 ,  0 ,  1 , 1  , 0  ],
        [1 , 0 , 0,  0 ,  0 , 0  ],
        [1 , 0 , 0 , 1 , 0 , 0  ],
        [1 ,  1  ,  1  , 1  , 1  , 0 ]
        ])
        sage: minimum_upset_random(M)
        (
        [0 0 1 1 0 1]                    
        [1 0 0 1 0 1]                    
        [0 1 0 0 0 0]                    
        [0 0 1 0 0 0]                    
        [1 1 1 1 0 1]                    
        [0 0 1 1 0 0], [1, 2, 0, 3, 5, 4]
        )

    """
    n = len(M.rows())
    Sn = SymmetricGroup(n)
    M1 = M
    wins = sum([sum([M1[j][i] for i in range(j,6)]) for j in range(6)])
    g0 = Sn(1)
    for k in range(N):
        g = Sn.random_element()
        P = g.matrix()
        M0 = P*M1*P^(-1)
        if sum([sum([M0[j][i] for i in range(j,6)]) for j in range(6)])>wins:
            M1 = M0
            g0 = g*g0
    return M1,g0(range(n))