Architected elastomers have demonstrated great potential for energy absorption, multi-resonant vibration isolation, and multi-bandgap acoustic control, due to the reversibility and programmability of their mechanical instabilities. However, computational design tools are needed to explore the large parameter space that regulates buckling-based mechanical instability behavior. In this study, we develop a machine-learning-based design optimization framework to control equilibrium states of a bistable elastomeric beam, while also regulating the required energy to transition between these configurations. Leveraging symmetry to reduce the design space, the research is performed on a single element, an inclined, slender beam, through a Fourier series parameterization. To evaluate the force-displacement response of the bistable beam, a nonlinear finite element analysis (FEA) with an arc-length continuation method is employed in a commercial FEA tool. Due to the highly non-convex bistability response of the beam in this proposed design space and the computational cost of the FEA analysis, a Bayesian optimization is implemented to promote a better trade-off between the number of function evaluations and rate of convergence. Bayesian optimization depends on several optimization parameters, which are systematically tuned in this study to characterize their role on the optimization process. With the proposed method, the equilibrium displacement and the ratio of output to input energy between stable states can be optimized within a few tens of iterations. A multi-objective optimization is also carried out to study the trade-off between equilibrium position and the energetics to transition between the bistabilities.