Fitting piecewise affine models to data points is a pervasive task in many scientific disciplines. In this work, we address the k- Piecewise Affine Model Fitting with Pairwise Linear Separability problem (k-PAMF-PLS) where, given a set of real points and the corresponding observations, we have to partition the real domain into k pairwise linearly separable subdomains and to determine an affine submodel (function) for each of them so as to minimize the total linear fitting error w.r.t. the observations. To solve k-PAMF-PLS to optimality, we propose a mixed-integer linear programming (MILP) formulation where symmetries are broken by separating the so-called shifted column inequalities. For medium-to-large scale instances, we develop a four-step heuristic involving, among others, a point reassignment step based on the identification of critical points and a domain partition step based on multicategory linear classification. Differently from traditional approaches proposed in the literature for similar fitting problems, in our methods the domain partitioning and submodel fitting aspects are taken into account simultaneously. Computational experiments on real-world and structured randomly generated instances show that, with our MILP formulation with symmetry breaking constraints, we can solve to proven optimality many small-size instances. Our four-step heuristic turns out to provide close-to-optimal solutions for small-size instances, while allowing to tackle instances of much larger size. The experiments also show that the combined impact of the main features of our heuristic is quite substantial when compared to standard variants not including them.