We all know that a POVM $\Pi$ is a set of operator such that $$ \Pi = \left\{\Pi_k \quad s.t.\quad \Pi_k \geq 0 \quad \& \quad \sum_{k=1}^n \Pi_k = \mathbb{I}\right\} $$ Let us assume that we are working with a qubit, so that $\mathbb{I} \to \mathbb{I}_2$ and $\Pi_k$ are $2\times2$ matrices. Let us also assume that $n=3$, so that we have three elements in the POVM and it can not be a projective measurement.
My question regards the following: how can I randomly generate such a POVM? Let us assume that we can write a generic positive operator such as $$ \Pi_k = r_k^0 \mathbb{I}_2 + r^i_k \sigma^i $$ where we used the Einstein notations and $\sigma^i$ are the Pauli matrices. Then the condition to be satisfied for $k=1,2$ for the vector $\{r^0_k,r^1_k,r^2_k,r^3_k\}$ is the following $$ r^0_k + \sqrt{{r_k^1}^2+{r_k^2}^2+{r_k^3}^2}>0 \quad \& \quad r^0_k - \sqrt{{r_k^1}^2+{r_k^2}^2+{r_k^3}^2}>0 $$ with the additional constraint that $$ \begin{cases} \sum_{k=1}^3 r^0_k = 1\\ \sum_{k=1}^3 r^j_k = 0 \quad \text{for } j=1,2,3 \end{cases} $$
So, the problem here is that I have a huge number of constraint that must be satisfied to randomly generate such a POVM and they are not feasible to be formulated in any informatical language (at least, as far as I now, the code on Mathematica is not working in small time i.e. 10 minutes...).
Is there any other way to formulate the problem in a more efficient way, or to reduce the random generation to more feasible things, such as unitary matrix and so on? I guess that someone has already tried it, but I was not able to find any results. I was thingking also to consider the possibility of randomly generating a PVM on a larger space and than tracing out some degrees of freedom and obtain a POVM, but I am not so sure is a right thing to do and if the sampling would be omogeneous or somehow there are some bias.