Architected elastomeric beam networks have great potential for energy absorption, multi-resonant vibration isolation, and multi-bandgap elastic wave control, due to the reconfigurability and programmability of their mechanical buckling instabilities. However, navigating this design space is challenging due to bifurcations between mono- and bistable beam designs, inherent geometric nonlinearities, and the strong dependence of buckling properties on beam geometry. To investigate these challenges, we developed a Bayesian optimization framework to control the equilibrium states of an inclined elastomeric beam, while also tuning the energy to transition between these configurations. Leveraging symmetry to reduce the design space, the beam shape is parameterized using a Fourier series representation. A penalty method is developed to include monostable designs in objective functions with dependencies on bistable features, enabling monostable results to still be incorporated in the Gaussian process surrogate and contribute to the optimization process. Two objectives are optimized in this study, including the position of the second stable equilibrium configuration and the ratio of output to input energy between the two stable states. A scalarized multi-objective optimization is also carried out to study the trade-off between equilibrium position and the energetics of transition between the stable states. The predicted designs are qualitatively verified through experimental testing. Collectively, the study explores a new parameter space for beam buckling, introduces a penalty method to regularize between mono- and bistable domains, and provides a library of beams as building blocks to assemble and analyze in future studies.