A highly efficient and accurate stochastic optimal control algorithm is proposed for active vehicle suspension systems. It uses the pseudo-excitation method for random vibration and the precise integration method (PIM) for the solution of the Riccati equations. Unlike conventional time-domain based methods, the present method computes the power spectral density (PSD) of each vehicle response directly from the road surface elevation PSD. Therefore root mean squares of the vehicle responses, including the frequency-weighted acceleration ones required by the ISO standard 2631–1(1997), can be obtained conveniently and accurately. The PIM used is proved to be numerically stable and highly precise, even if all eigenvalues of the corresponding Hamiltonian matrices are very close to the imaginary axis. Finally, numerical examples verify the effectiveness of the proposed method.