The computational efficiency of the Finite-Difference Time-Domain (FDTD) method can be significantly reduced by the presence of complex objects with fine features. Small geometrical details impose a fine mesh and a reduced time step, significantly increasing computational cost. Model order reduction has been proposed as a systematic way to generate compact models for complex objects, that one can then instantiate into a main FDTD mesh. However, the stability of FDTD with embedded reduced models remains an open problem. We propose a systematic method to generate reduced models for FDTD domains, and embed them into a main FDTD mesh with guaranteed stability. Models can be created for arbitrary domains containing inhomogeneous and lossy materials. The Courant-Friedrichs-Lewy (CFL) limit of the final scheme is provided by the theory, and can be extended with a simple perturbation of the coefficients of the reduced models. Numerical tests confirm the stability of the proposed method, and its potential to accelerate multiscale FDTD simulations.