The energy-based variational theory has been used to model multi-fracture growth in XFEM, which introduces the energy Hessian matrix to characterize the competition propagation among branch fractures. This framework, relying on Newton method under smoothness assumption, fails in layered heterogeneous structure because energy-gradient discontinuity at interfaces produces unreliable search steps that deviate from the energy descent direction. In this paper, a novel variational XFEM with aggregating historical energy gradient information is developed for multi-fracture propagation in heterogeneous formation. During the iteration process, this method incorporates the gradient information from historical iteration path to construct a piecewise quadratic model, which determines the search step with energy-optimal gradient aggregating weigh in non-smooth regions. Once the iteration leaving the non-smooth regions, the resulting search step can degenerate to the standard Newton form, achieving a Newton-like rapid convergence rate. In addition, the unique null-step mechanism of this method ensures a stable and sufficient total energy decline at each iteration, thereby improving the overall solving efficiency. This method is validated through two benchmark fracture propagation cases for its accuracy, and the results of energy variation and iteration convergence from some general multi-fracture growth cases in layered heterogeneous structure show its superior performance in the total energy minimization with gradient discontinuity. The proposed method provides a generalized XFEM framework for modeling complex geometric and physical interactions among each branch fracture in a multi-fracture system of a realistic engineering structure.