Determining the energetically most favorable structure of nanoparticles is a fundamentally important task towards understanding their stability. In the case of bimetallic nanoclusters, their vast configurational space makes it especially challenging to find the global energy optimum via experimental or computational screening. To that end, this work proposes a two-step optimization-based design framework to address this hard combinatorial problem. Given a nanocluster of fixed shape, a rigorous mixed-integer linear programming model is formulated based on a bond-centric cohesive energy function to identify the most cohesive bimetallic configuration for a given composition. This capability is coupled with a metaheuristic strategy that searches over the space of nanocluster shapes to obtain optimal structures. We apply our proposed methodology on AgCu, AuAg and CuAu systems, quantifying how the size and composition of a nanocluster influences its overall cohesion. Furthermore, we observe various synergistic effects between Cu and Au in promoting cohesive energy, while multiple segregation patterns are identified in all three studied binary systems. Our methodology serves as an efficient computational tool for investigating bimetallic nanoclusters stability properties as well as provides model nanoclusters for further investigations.