In this work, we introduce a hybrid computational scheme for modeling multistate nonadiabatic dynamics of molecules that treats population transfer among electronically excited states with a linear-response approach and S1 -> S0 dynamics with a state-averaged method. Specifically, a time-dependent density functional theory (TDDFT)-based method and an ensemble density functional theory-based method are employed to model the gas-phase photodynamics of trans- and cis-azobenzene. The proposed computational protocol accurately reproduces the (approximately) 2-fold difference in the quantum yield of photoisomerization for the trans-azobenzene excited in the n pi* and pi pi* states, in agreement with other theoretical simulations and experimental data. These results demonstrate that the proposed hybrid computational scheme provides a reliable description of multistate photochemical processes in molecular systems.