We propose and implement a computational procedure to establish data-driven surrogate
constitutive models for heterogeneous materials. We study the multiaxial response of non-linear
n-phase composites via Finite Element (FE) simulations and computational homogenisation.
Pseudo-random, multiaxial, non-proportional histories of macroscopic strain are imposed on
volume elements of n-phase composites, subject to periodic boundary conditions, and the
corresponding histories of macroscopic stresses and plastically dissipated energy are recorded.
The recorded data is used to train surrogate, phenomenological constitutive models based on
neural networks (NNs), and the accuracy of these models is assessed and discussed. We analyse
heterogeneous composites with hyperelastic, viscoelastic or elastic-plastic local constitutive
descriptions. In each of these three cases, we propose and assess optimal choices of inputs and
outputs for the surrogate models and strategies for their training. We find that the proposed
computational procedure can capture accurately and effectively the response of non-linear nphase composites subject to arbitrary mechanical loading.