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 Piecewise Linear Separability problem (k-PAMF-PLS) where, given a set of m points { a 1 , ź , a m } ź R n and the corresponding observations { b 1 , ź , b m } ź R , we have to partition the domain R n into k piecewise linearly (or affinely) separable subdomains and to determine an affine submodel (i.e., an affine function) for each of them so as to minimize the total linear fitting error w.r.t. the observations bi.To solve k-PAMF-PLS to optimality, we propose a Mixed-Integer Linear Programming (MILP) formulation where symmetries are broken by separating shifted column inequalities. For medium-to-large scale instances, we develop a four-step heuristic involving, at each iteration, 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 both our exact and heuristic 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 the 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.We conclude the paper with an application to the identification of dynamical piecewise affine systems, for which we obtain promising results of comparable quality with those achieved with state-of-the-art methods from the literature (hinging hyperplane models and sum-of-norms regularization) on benchmark data sets. HighlightsWe address the fundamental piecewise affine model fitting problem.We propose an MILP formulation, strengthened with symmetry breaking constraints.We present a heuristic with point-reassignment and domain partition at each iteration.We report extensive results on real-world and structured random instances.We conclude with an application to the identification of hybrid dynamical systems.