What we want is to generate a bistochastic matrix according to the Haar measure, which is the unique distribution which is invariant to multiplication by bistochastic matrices from both sides.
The standard algorithm is to take some iid matrix (each entry is chosen iid from some distribution over non-negative numbers) and then repeatedly make it row-stochastic and column-stochastic - this is like projecting the matrix to the linear subspace of stochastic matrices. The process converges pretty fast, and we get some random bistochastic matrix, even though not according to the wanted distribution.
A recent paper discusses these issues in the context of estimating the measure of all bistochastic matrices (an open problem).