Radiation with high dose rate (FLASH) has shown to reduce toxicities to normal tissues around the target and maintain tumor control with the same amount of dose compared to conventional radiation. This phenomenon has been widely studied in electron therapy, which is often used for shallow tumor treatment. Proton therapy is considered a more suitable treatment modality for deep-seated tumors. The feasibility of FLASH proton therapy has recently been demonstrated by a series of pre- and clinical trials. One of the challenges is to efficiently generate wide enough dose distributions in both lateral and longitudinal directions to cover the entire tumor volume. The goal of this paper is to introduce a set of automatic FLASH proton beam optimization algorithms developed recently. To develop a fast and efficient optimizer for the design of a passive scattering proton FLASH radiotherapy delivery at The University of Texas MD Anderson Proton Therapy Center, based on the fast dose calculator (FDC). A track-repeating algorithm, FDC, was validated versus Geant4 simulations and applied to calculate dose distributions in various beamline setups. The design of the components was optimized to deliver homogeneous fields with well-defined diameters between 11.0 and 20.5mm, as well as a spread-out Bragg peak (SOBP) with modulations between 8.5 and 39.0mm. A ridge filter, a high-Z material scatterer, and a collimator with range compensator were inserted in the beam path, and their shapes and sizes were optimized to spread out the Bragg peak, widen the beam, and reduce the penumbra. The optimizer was developed and tested using two proton energies (87.0 and 159.5MeV) in a variety of beamline arrangements. Dose rates of the optimized beams were estimated by scaling their doses to those of unmodified beams. The optimized 87.0-MeV beams, with a distance from the beam pipe window to the phantom surface (window-to-surface distance [WSD]) of 550mm, produced an 8.5-mm-wide SOBP (proximal 90% to distal 90% of the maximum dose); 14.5, 12.0, and 11.0-mm lateral widths at the 50%, 80%, and 90% dose location, respectively; and a 2.5-mm penumbra from 80% to 20% in the lateral profile. The 159.5-MeV beam had an SOBP of 39.0mm and lateral widths of 20.5, 15.0, and 12.5mm at 50%, 80%, and 90% dose location, respectively, when the WSD was 550mm. Wider lateral widths were obtained with increased WSD. The SOBP modulations changed when the ridge filters with different characteristics were inserted. Dose rates on the beam central axis for all optimized beams (other than the 87.0-MeV beam with 2000-mm WSD) were above that needed for the FLASH effect threshold (40Gy/s) except at the very end of the depth dose profile scaling with a dose rate of 1400Gy/s at the Bragg peak in the unmodified beams. The optimizer was able to instantly design the individual beamline components for each of the beamline setups, without the need of time intensive iterative simulations. An efficient system, consisting of an optimizer and an FDC have been developed and validated in a variety of beamline setups, comprising two proton energies, several WSDs, and SOBPs. The set of automatic optimization algorithms produces beam shaping element designs efficiently and with excellent quality.