  We propose an efficient numerical method to compute configuration averages of
observables in disordered open quantum systems whose dynamics can be unraveled
via stochastic trajectories. We prove that the optimal sampling of trajectories
and disorder configurations is simply achieved by considering one random
disorder configuration for each individual trajectory. As a first application,
we exploit the present method to the study the role of disorder on the physics
of the driven-dissipative Bose-Hubbard model in two different regimes: (i) for
strong interactions, we explore the dissipative physics of fermionized bosons
in disordered one-dimensional chains; (ii) for weak interactions, we
investigate the role of on-site inhomogeneities on a first-order dissipative
phase transition in a two-dimensional square lattice.
