We study the constrained system of linear equations

$Ax=b$, $x∈\mathcal{R}(A^k)$

for  $A ∈ \mathbb{C}^{n×n}$  and  $b ∈\mathbb{C}^n$, $k=Ind(A)$. When the system is consistent, it is well known that it has a unique $A^Db$. If the system is inconsistent, then we seek for the least squares solution of the problem and consider
$$\min _{x \in \mathcal{R}\left(A^{k}\right)}\|b-A x\|{_2,}$$

where $\|\cdot \|_2$ is  the 2-norm. For the inconsistent system with a matrix $A$ of index one, it was proved recently that the solution is $A^⊕b$ using the core inverse $A^⊕$ of $A$. For matrices of an arbitrary index and an arbitrary $b$, we show that the solution of the constrained system can be expressed as $A^⊕b$ where $A^⊕$ is the core-EP inverse of $A$. We establish two Cramer's rules for the inconsistent constrained least squares solution and develop several explicit expressions  for the core-EP inverse of matrices of an arbitrary index. Using these expressions, two Cramer's rules and one Gaussian elimination method for computing the core-EP inverse of matrices of an arbitrary index are proposed in this paper. We also consider the $W$-weighted core-EP inverse of a rectangular matrix and apply the weighted core-EP inverse to a more general constrained system of linear equations.