This paper presents complete nonlinear electromechanical models for energy harvesting devices consisting of multiple piezoelectric bimorphs (PBs) connected in parallel and series, for the first time. The proposed model is verified against available experimental results for a specific case. The piezoelectric and beam constitutive equations and different circuit equations are utilized to derive the complete nonlinear models for series and parallel connections of the PBs as well as those of piezoelectric layers in each bimorph, i.e., four nonlinear models in total. A multi-modal Galerkin approach is used to discretize these nonlinear electromechanical models. The resultant high-dimensional set of equations is solved utilizing a highly optimized and efficient numerical continuation code. Examining the system behavior shows that the optimum load resistance for an energy harvester array of 4 PBs connected in parallel is almost 4% of that for the case with PBs connected in series. It is shown an energy harvesting array of 8 PBs could reach a bandwidth of 14 Hz in low frequency range, i.e., 20–34 Hz. Compared with an energy harvester with 1 PB, it is shown that the bandwidth can be increased by more than 300% using 4 PBs and by more than 500% using 8 PBs. Additionally, the drawbacks of a multi-PB energy harvesting device are identified and design enhancements are proposed to improve the efficiency of the device.