Vibration-based piezoelectric energy harvesting has received significant attention recently due to its ability to convert vibration energy from the surrounding environment into electric energy to power wireless sensors. In this work, two typical types of vibration sources, that is, engine-induced floor vibration and train-induced bridge vibration, are considered. The array of piezoelectric energy harvesters (PEHs) is adopted to harvest the energy of vibration with different frequencies. The aim of this work is to give a stable system-level modeling method and present a method to select operating point for designing energy harvesters connected with nonlinear energy harvesting circuits under complex excitations. The system-level modeling method is based on Krylov’s subspace model order reduction algorithm. With the help of this method, the harvester can be modeled using conventional finite element software and co-simulated with nonlinear circuits at a system-level. The proposed method is validated using equivalent circuit method, equivalent impedance method and experiment results, and then employed to investigate the performances of harvesters under typical real vibrations. Finally, design suggestions for PEHs under different types of vibration sources are presented.