In this paper, an iterative method is constructed to solve the following constrained linear matrix equations system: [A(1)(X), A(2)(X), ..., A(r)(X)] = [E-1, E-2, ..., E-r], X 2 S = {X|X = U(X)}, where A(i) is a linear operator from C-mxn onto C-pixqi, E-i is an element of C-pixqi, i = 1, 2, ..., r, and U is a linear self-conjugate involution operator. When the above constrained matrix equations system is consistent, for any initial matrix X0 2 S, a solution can be obtained by the proposed iterative method infinite iteration steps in the absence of roundoff errors, and the least Frobenius norm solution can be derived when a special kind of initial matrix is chosen. Furthermore, the optimal approximation solution to a given matrix can be derived. Several numerical examples are given to show the efficiency of the presented iterative method.