A highly applicable analytical approach is proposed for the dynamic analysis of composite plate structures with complex variation of geometric characteristics (shape/boundary and thickness) and material characteristics (Young's modulus, shear modulus and density) resting on Winkler-Pasternak foundation with general elastic boundary conditions. Boundary conditions of the plates are treated by the penalty method, so that the admissible functions can use any complete set of orthogonal polynomials, which is obtained by Gram-Schmidt orthogonalization process. The interval where the in-plane coordinate variables of the plate domain are located can be selected as the orthogonalization intervals, so as to avoid the normalization procedure of coordinates and simplify the operation process. In order to calculate the energy integral of the arbitrarily shaped plate, the plate domain segmentation strategy is introduced. To verify the effectiveness and accuracy of the proposed method, convergence study and numerical verifications for different orthogonalization intervals, weight functions, penalty spring stiffness and the truncation number of orthogonal polynomials are carried out. Taking regular hexagonal plates as examples, the effects of the parameters of thickness, nonhomogeneity and foundations on the vibration characteristics of the plates are studied. These results can be used as reference data for future research.