The traditional power system dominated by synchronous generators is gradually evolving into a modern power system featured by high-penetrated renewable energy. As a key technology for high-penetrated renewable energy, the grid-forming voltage source converter (GFM-VSC) has received increasing attention. However, the large-disturbance stability analysis of power systems with multiple GFM-VSCs is still a challenging problem due to various limitations of existing methods, including huge computational burden and difficulty in considering network losses. This paper is intended to address these issues from the perspective of reduced-order modeling and domain of attraction (DA) estimation. The innovations involve three aspects. First, the reduced-order modeling method for power systems with multiple GFM-VSCs is proposed using the standard dual-time-scale model in singular perturbation theory. Second, an expanding annular domain (EAD) algorithm is developed to estimate the DA with an entire boundary to analyze the large-disturbance stability of power systems. Third, the conditions of using the reduced-order modeling method based on singular perturbation theory have been clarified. The validity of the reduced-order modeling method is illustrated on a modified 39-bus system with 10 GFM-VSCs.