We are concerned with the overall elastic energy in martensitic polycrystals. These are polycrystals whose constituent crystallites can undergo shape-deforming phase transitions as a result of changes in their stress or temperature. We approach the problem of calculation of the nonlinear overall energy via a statistical optimization method which involves solution of a sequence of linear elasticity problems. As a case study we consider simulations on a two-dimensional model in which circular randomly-oriented crystallites are arranged in a square pattern within an elastic matrix. The performance of our present code suggests that this approach can be used to compute the overall energies in realistic three-dimensional polycrystals containing grains of arbitrary shape. In addition to numerical results we present upper bounds on the overall energy. Some of these bounds apply to the square array mentioned above. Others apply to polycrystals containing circular, randomly-oriented crystallites with sizes ranging to infinitesimal, and no intergrain matrix. The square-array bounds are consistent with our numerical results. In some regimes they approximate them closely, thus providing an insight on the convergence of the numerical method. On the other hand, in the case of the random array the bounds carry substantial practical significance, since in this case the energy contains no artificial contributions from an elastic matrix. In all the cases we have considered our bounds compare favorably with those obtained under the well-known Taylor hypothesis; they show that, as far as polycrystalline martensite is concerned, calculations of the elastic energy based on the Taylor assumption may lead to substantial overestimates of this quantity.