Recent advances in computational mechanics and machine learning have enabled increasingly realistic simulations of tumor growth and therapeutic response, opening new possibilities for data-informed treatment optimization. In this work, we present a unified computational framework for the stochastic optimization of immunotherapy design and scheduling, built upon a high-fidelity mechanistic model of tumor–host interactions. Tumor progression and immune response are described using a biphasic porous-media formulation that accounts for tissue mechanics, transport phenomena, and the effects of anti-PD-L1 immunotherapy, while explicitly incorporating parametric uncertainty in tumor and vascular properties. Experimental observations are assimilated through Bayesian inference to probabilistically update key model parameters, providing a quantified description of uncertainty in tumor dynamics. Building on these probabilistic predictions, a stochastic optimization problem is formulated to identify optimal immunotherapy characteristics and treatment schedules. The optimization is solved using a genetic algorithm, enabling efficient exploration of a high-dimensional design space. To overcome the prohibitive computational cost associated with repeated evaluations of the nonlinear multiphysics model, a Deep Operator Network (Deep-O-Net) surrogate is employed to accurately approximate the model response. Numerical results demonstrate that the proposed framework achieves substantial improvements in tumor suppression and treatment efficacy compared to non-optimized protocols. The presented methodology provides a robust and extensible platform for uncertainty-aware, personalized immunotherapy planning and highlights the potential of combining mechanistic modeling, Bayesian learning, and operator-based machine learning in cancer treatment optimization.