The equilibrium concentration of polyvacancy in crystal is formulated by modifying the grand canonical ensemble. The generalized equation can be applied to obtain the concentration of vacancy in any order and shape in crystal. In this work, a Widom-like particle insertion method is adopted implemented with molecular dynamics simulation to obtain individual formation free energies for constituent vacancies in polyvacancy. As a case study, equilibrium concentrations and formation free energies were obtained of four and nine trivacancies identified in face-centred-cubic (FCC) and hexagonal-close-packed (HCP) hard-sphere crystals, respectively. The result is in excellent agreement with literature data available for the FCC. Further, stabilities of equilateral triangular trivacancies are elucidated and their relative stabilities in FCC and HCP hard-sphere crystals are reported. Overall trivacancy concentration is higher in the HCP HS crystal.