Let $A\in\{0,1\}^{(n,n)}$ be a $n$ by $n$ boolean matrix (in particular think of an adjacency matrix of a graph), and consider the following optimization problem:

$$\begin{align*}&&\max_{P\in\{0,1\}^{(n,n)}}& \sum_{i =1}^n \sum_{j=1}^n (P^{t}AP)_{i,j} \\ &&\text{s.t}.\hspace{5mm}& (\vec{1})^{t}P=(\vec{1})^{t}\\ &&& \hspace{4.4mm}P\hspace{0.7mm}\vec{1} = \vec{1}\end{align*}$$

This optimization problem is asking how to reorder the rows/columns of a matrix through a permutation $P\in S_{n}$ (with its representation as a matrix) in such a way that the sum of the upper triangular is maximal, this is known to be an NP-Hard problem for adjacency matrices of general digraphs, but from my searches it seems that it is not yet known for tournaments (i.e. orientations of the complete graph).

This can also be implemented as a linear problem noting that $(P^{t}AP)_{i,j} = A_{\pi(i),\pi(j)}=\sum\limits_{k=1}^{n}\sum\limits_{l=1}^{n} A_{k,l}\min\{P_{k,i},P_{l,j}\}$, where $\pi$ is the permutation (as a function $\pi:[n]\to[n])$ given by the matrix $P$.

I've implemented the linearization on PySCIPOpt (Github link here), but it's really slow even for small matrices. I guess it's an inherent problem of optimizing over permutations, is there any "quick" way to solve (or optimize the solving algorithm) of this kind of problems?.