The paper presents the methods for solving the homogenization problems for two-phase piezoelectric composites, which are implemented in ACELAN-COMPOS finite element package, developed by the authors. The considered composites consist of piezoelectric skeleton and elastic inclusions. The effective moduli method is used to determine the effective properties of the composite. In this method, the static electroelasticity problems with special boundary conditions are set for a representative volume of the composite. These problems are solved numerically using the finite element method. The developed algorithm for generating representative volumes for the composites with granular inclusions is described in details. The work of the algorithm is illustrated by an example of a composite made of PZT-4 piezoceramic with inclusions of α-corundum.