Chemotherapy is one of the primary modalities of cancer treatment. Chemotherapy drug administration is a complex problem that often requires expensive clinical trials to evaluate potential regimens. One way to alleviate this burden and better inform future trials is to build reliable models for drug administration. Previous chemotherapy optimization models have mainly relied on optimal control, which does not lend itself to capturing complex and vital operational constraints in chemotherapy planning involving discrete decisions, such as doses via pills and rest periods. In addition, most of the existing models for chemotherapy optimization lack an explicit toxicity measure and impose toxicity constraints primarily through (fixed) limits on drug concentration. The existing stochastic optimization models also focus on maximizing the probability of cure when tumor heterogeneity is uncertain. In this paper, we develop a mixed-integer program for combination chemotherapy (utilization of multiple drugs) optimization that incorporates various important operational constraints and, besides dose and concentration limits, controls treatment toxicity based on its effect on the count of white blood cells. To address the uncertainty of tumor heterogeneity, we propose chance constraints that guarantee reaching an operable tumor size with a high probability in a neoadjuvant setting. We present analytical results pertinent to the accuracy of the model in representing biological processes of chemotherapy and establish its merit for clinical applications through a numerical study of breast cancer.