Unresolved binaries in star clusters can bias stellar and cluster mass estimates, making their proper treatment essential for studying cluster dynamics and evolution. We aim to develop a fast and robust framework for jointly deriving stellar masses and multiplicity statistics of member stars, together with optimal cluster parameters. We use Gaia DR3 parallaxes together with multi-band photometry of open cluster (OC) members to infer stellar masses and binary mass-ratios through simulation-based inference (SBI), while iteratively fitting the cluster parameters. The validation of our SBI framework on simulated clusters demonstrates that the inclusion of infrared photometry significantly improves the detection of low mass-ratio binaries. The minimum mass-ratio threshold for reliably identifying unresolved binaries depends on cluster properties and the available photometry, but typically lies below q=0.5. Applying our method to 42 well-populated OCs, we derive a catalogue of stellar masses and mass-ratios for 27201 stars, achieving typical uncertainties of 0.08 in q and 0.01M_{Sun}_ in the primary stellar mass. We analyse the archetype OCs M67 and NGC 2360 in detail, including mass segregation and mass-ratio distribution among other characteristics, while deriving multiplicity fractions for the rest of the sample. We find evidence that the high mass-ratio (q=>0.6) binary fraction shows a strong correlation with the age and a weak anti-correlation with the cluster metallicity. Furthermore, the variation of the binary fraction with stellar mass in OCs shows strong accordance with the observed dependence for field stars heavier than ~>0.6M_{Sun}_. Our work paves a path for future population-level investigations of multiplicity statistics and precision stellar masses in extended samples of OCs.