Amplitude equations of Landau type, which describe the dynamics ofthe most amplified periodic disturbance waves in slightly supercritical flow systems, have been known to form reliable and sufficiently accurate low-dimensional models capable of predicting the asymptotic magnitude of saturated perturbations. However the derivation of similar models for estimating the threshold disturbance amplitude in subcritical systems faces multiple resonances which lead to the singularity of model coefficients. The observed resonances are traced back to the interaction between the mean flow distortion induced by the decaying fundamental disturbance harmonic and other decaying disturbance modes. Here we suggest a methodology of deriving a two-equation dynamical system of coupled cubic amplitude equations with non-singular coefficients which resolve the resonances and are capable of predicting the threshold amplitude for weakly nonlinear subcritical regimes. The suggested reduction procedure is based on the consistent use of an appropriate orthogonality condition which is different from a conventional solvability condition. As an example, a developed procedure is applied to a system of Navier-Stokes equations describing a subcritical plane Poiseuille flow. Predictions of the so-developed model are found to be in reasonable agreement with experimentally detected threshold amplitudes reported in literature.