Radiation therapy (RT) aims to deliver tumoricidal doses with minimal radiation-induced normal-tissue toxicity. Compared to conventional RT (of conventional dose rate), FLASH-RT (of ultra-high dose rate) can provide additional normal tissue sparing, which however has created a new nonconvex and nonsmooth optimization problem that is highly challenging to solve. In this paper, we propose a stochastic three-operator splitting (STOS) algorithm to address the FLASH optimization problem. We establish the convergence and convergence rates of the STOS algorithm under the nonconvex framework for both unbiased gradient estimators and variance-reduced gradient estimators. These stochastic gradient estimators include the most popular ones, such as SGD, SAGA, SARAH, and SVRG, among others. The effectiveness of the STOS algorithm is validated using FLASH radiotherapy planning for patients.