Consider a service system where incoming tasks are instantaneously dispatched to one out of many heterogeneous server pools. Associated with each server pool is a concave utility function that depends on the class of the server pool and its current occupancy. We derive an upper bound for the mean normalized aggregate utility in stationarity and introduce two load balancing policies that achieve this upper bound in a large-scale regime. Furthermore, the transient and stationary behavior of these asymptotically optimal load balancing policies is characterized on the scale of the number of server pools in the same large-scale regime. Funding: This work was supported by the Netherlands Organization for Scientific Research (NWO) through [Gravitation Grant NETWORKS-024.002.003] and [Gravitation Grant Vici 202.068]. Supplemental Material: The online appendix is available at https://doi.org/10.1287/stsy.2022.0103 .