We investigate whether viable Hu-Sawicki-like $f(R)$ models can produce deviations from $\Lambda\mathrm{CDM}$ that can be tested against current background cosmological data. We adopt a machine-learning approach based on Genetic Algorithms (GA) to reconstruct analytical perturbations around the Hu-Sawicki class of models. We develop a pipeline that interfaces the \texttt{GATO} GA code with the \texttt{CANDI} cosmology code. Each $f(R)$ function generated by the GA is first tested against theoretical viability conditions, including stability, the recovery of a standard matter-dominated epoch, the General Relativity limit, and chameleon screening mechanism. Viable candidates are then passed to \texttt{CANDI} to reconstruct the corresponding background cosmology and are tested against DESI DR2 BAO measurements and the Pantheon+ Type Ia supernova catalogue. %\newline The deviations we find are largest at late times, where the lower curvature makes modified-gravity effects more relevant, and are rapidly suppressed at higher redshift, in agreement with the imposed matching to the matter-dominated era. To further quantify deviations from the standard cosmological model, we compute the $Om(z)$ diagnostic. It shows only a very small departure from the constant $\Lambda\mathrm{CDM}$ behaviour. The effective dark energy equation of state associated with the reconstructed $f(R)$ function also evolves only weakly, showing a mild transition from an effective quintessence-like nature to an effective phantom-like regime. Overall, our results indicate that, within perturbations around the Hu-Sawicki class of models, current background data allow only limited deviations from $\Lambda\mathrm{CDM}$.