This paper presents a method of model reduction for interval systems using the Padé approximation to allow retention-dominant poles, where the denominator dm(s) of the reduced model is formed by retaining the dominant poles of the original interval system, and the numerator nm (s) of the reduced model is obtained by matching the first r interval time moments of the reduced model with that of the system. A numerical example illustrates the procedure.