This paper deals with developing a robust iterative algorithm to find the least-squares (P, Q)-orthogonal symmetric and skew-symmetric solution sets of the generalized coupled matrix equations. To this end, first, some properties of these type of matrices are established. Furthermore, an approach is offered to determine the optimal approximate (P, Q)-orthogonal (skew-)symmetric solution pair corresponding to a given arbitrary matrix pair. Some numerical experiments are reported to confirm the validity of the theoretical results and to illustrate the effectiveness of the proposed algorithm.