Designing icephobic surfaces to delay ice formation is crucial for applications like aviation safety and cryopreservation. While Classical Nucleation Theory (CNT) provides a thermodynamic foundation, real-world stochastic effects and complex wetting states make designing these surfaces difficult. The performance of icephobic surfaces is influenced by surface structure. However, finding the best designs for specific environmental conditions is still a challenge. In this work, we introduce a hybrid approach that combines experiments, CNT, and blackbox optimization to predict and optimize ice nucleation time on micropatterned surfaces. Cylindrical SU-8 micropillar arrays with different heights and spacings were created, and the apparent contact angles and freezing delay times were measured at \(-10~^\circ\mathrm{C}\) and \(-20~^\circ\mathrm{C}\). An analytical model was developed to describe the wetting states between the Wenzel and Cassie Baxter regimes, optimizing its parameters using the Mesh Adaptive Direct Search (MADS) algorithm. This approach allows us to estimate contact angles and ice nucleation times for any surface geometry within the studied design space. The predicted contact angles matched experimental results with a mean absolute percentage error (MAPE) of 2.09% (\(R^2\) = 0.92) and the approximate nucleation times showed a MAPE of 27.3% (\(R^2\) = 0.75). Our method also identified the best micropillar geometries that maximize freezing delay and was validated with an independent dataset, showing strong predictive ability. This work emphasizes the benefit of combining physics-based models with data-driven optimization for rapid design of icephobic surfaces.