In this work, we developed a higher-order asymptotic homogenization method (AHM) for efficient multiscale analysis of piezoelectric composite structures. Instead of solving the composite structural problems directly with all heterogeneities considered, AHM first solves the microscale problem at each order of expansion to obtain the corresponding microfluctuation functions and effective properties; the latter is then passed to the macroscale to solve coupling governing differential equations, with the macroscale physical quantities passed back to the microscale and combined with microfluctuation functions to recover the local response. Governing differential equations are formulated and solved in Cartesian coordinates, which is much beneficial for understanding the multiscale-multiphysics behavior of flat piezocomposite structures such as Micro Fiber Composite transducers. Two numerical examples are solved using AHM, i.e., a 1D piezocomposite rod and a 2D piezocomposite plate, and the results are compared with direct numerical solution (DNS) to validate the accuracy and efficiency of the proposed method. The developed multiscale analysis algorithm largely saves the computational cost compared with DNS when a large number of inclusions are included in piezocomposites, thus making it an efficient tool for solving multiphysics composite structural problems with piezoelectricity.