I would suggest you try out parameterizing your transition matrix via a product of spheres -- maybe not what you had in mind, but it has the merit of fitting into manopt fairly easily.
Specifically, let problem.M = obliquefactory(n, n); This defines a product of n spheres in R^n.
Then, each matrix X on that manifold is a matrix of size n x n with the property that each column has unit 2-norm. As a result, if you consider P = X.^2 (that is, entry-wise square the matrix), then that matrix P has nonnegative entries and its columns sum to 1.
If you would rather that it is the rows that sum to 1, you can call obliquefactory(n, n, true); (which transposes the matrices to fit the other common convention).
Now, you can specify your cost function problem.cost = @(X) .... and in the expression you would write X.^2 (or X.*X) instead of P.
My sense is that this might work quite well, but I'm curious to hear what you find.
finding the Riemannian gradient (which I believe is determined using the Euclidean gradient in this package)
Manopt allows you to define the Riemannian gradient directly, for example by setting problem.grad = ...; But it's true that most users (myself included) define the Euclidean gradient (or more generally the "gradient in the embedding space") by setting problem.egrad = ...; then counting on Manopt to do the appropriate transformation.