Complex physical or numerical systems may exhibit distinct behaviors in various zones of their design spaces. We present an algorithm that uses multiple cluster-based surrogates for optimizing such box-constrained systems. It partitions the design space into multiple clusters using K-means clustering and develops a separate surrogate for each cluster. It then uses these surrogates to sample additional points in the design space whose function evaluations guide the search for a global optimum. Clustering, surrogate construction, and smart sampling are employed iteratively to add sample points until a pre-defined threshold. The best solution from these points estimates a global optimum. An extensive test bed of 52 box-constrained functions was used to evaluate and compare the algorithm's performance and computational requirements with sixteen derivative-free optimization solvers. The best version of our algorithm surpassed all sixteen solvers in optimization accuracy for a fixed number of evaluations and demanded lower computational effort than fifteen.