Model reduction methods for bilinear control systems are compared by means of practical examples of Liouville-von Neumann and Fokker-Planck type. Methods based on balancing generalized system Gramians and on minimizing an H2-type cost functional are considered. The focus is on the numerical implementation and a thorough comparison of the methods. Structure and stability preservation are investigated, and the competitiveness of the approaches is shown for practically relevant, large-scale examples.