This paper investigates the nonlinear electromechanical behaviour of a constrained bimorph piezoelectric energy harvester utilising a coupled multi-modal fully nonlinear model, for the first time. The electromechanical system is modelled via using the nonlinear beam theory of Euler-Bernoulli, while assuming an inextensible centreline, and the piezoelectric constitutive equations. The motion constraints are modelled as nonlinear springs, consisting of linear and cubic terms. The coupled electromechanical model, including the electrical circuit equation, is obtained for both parallel and series connections of the piezoelectric layers. The Galerkin scheme is used to discretise the partial differential-type model while retaining a large number of modes to guarantee converged results. Next, the electromechanical model consisting of geometric, inertial, impact, and piezoelectric-type nonlinearities is solved using a well-optimised continuation code. The proposed model and numerical method are verified through comparison to experimental data in the literature. It is shown that adding a tip mass of 50% of the bimorph cantilever mass increase the power output by 263%. Additionally, it is shown that for a power output threshold of 15 mW, the addition of motion constraints results in a 66.8% resonance bandwidth (normalised with respect to first natural frequency) which is almost 10 times the bandwidth of the unconstrained system.