Summary: We consider a multiphase material which consists of a homogeneous elastic-plastic matrix containing a homogeneous statistically uniform random set of ellipsoidal elastic-plastic inclusions. The elastic properties of the matrix and the inclusions are the same, but the so-called “stress-free strains”, i.e. the strain contributions due to temperature loading, phase transformations, and the plastic strains, fluctuate. We employ a general theory of yielding for arbitrary loading (by stress and by temperature). The realization of an incremental plasticity scheme is based on averaging over each component of the nonlinear yield criterion. Usually, averaged stresses are used inside each component for this purpose. In distinction from this usual practice, we apply here physically consistent assumptions on the dependence of these functions on the components of the second stress moments. The application of the proposed theory to the prediction of thermomechanical deformation behavior of a model material is presented as an example.