Directed Acyclic Graphs (DAGs) are informative graphical outputs of causal learning algorithms to visualize the causal structure among variables. In practice, different causal learning algorithms are often used to establish a comprehensive analysis pool, which leads to the challenging problem of ensembling the heterogeneous DAGs with diverse and conflicting embedded information. While naive approaches such as MEDIAN and greedy heuristics are presented in the literature, a mathematical framework to build the ensemble DAG is not investigated. Therefore, we propose a two-stage ensemble framework for causality learning with heterogeneous DAGs. In the first stage, we implement a data partitioning procedure to categorize the input data. Then, we apply multiple causal learning algorithms to each class and ensemble the results across the partitions for each method. The results across different approaches in the second stage are then ensembled to produce the final DAG. We define a novel mathematical model based on the marginal contribution concept, reflecting the collective information for each edge from input graphs. We investigate the relationship between the proposed modeling approach and well-known problems such as rank aggregation and the Traveling Salesman Problem. We also present a parametric analysis for the distance function, enriched by insights from their asymptotic behavior. We design a simulation experiment to show the computational advantage of adopting a lazy constraints solution approach. In addition to the T-cell synthetic dataset, we use the dataset of water main breaks from 2017 to 2021 in the City of Tampa, Florida, to demonstrate the framework’s applicability. The results show that the optimal solution of the proposed framework outperforms current approaches in the literature based on causality-based Structural Intervention Distance performance measures.