We report simulations of exciton-exciton annihilation (EEA) in dye molecular aggregates with explicit laser pulse excitation in the framework of reduced density matrix theory. We use a supramolecular complex model Hamiltonian, in which each monomer of an aggregate is described by three electronic levels. An expectation value approximation is used to avoid the large amount of matrix elements for various exciton configuration distributions in molecular aggregates. Quantum fluctuations truncation is applied to derive the closed equations of motion for single- and two-site functions. Using these equations, we simulate the laser pulse driven EEA dynamics, and compare it with rate equations to highlight the importance of quantum coherence effects. We also discuss the effects of the aggregate configuration, the molecular transition dipoles and the duration of laser pulses on EEA. Furthermore, we analyze the relation between the coherence time of fusion and the annihilation process. Our study indicates that EEA process is more efficient for coherent excitation at the lowest energy level in J aggregates. Our simulation of the normalized population of the first excited state for various field amplitudes can qualitatively reproduce time-resolved experimental emission spectra. This confirms that the fundamental description of laser pulse induced EEA is correct within our model.