Shape optimization is a rich field at the intersection of analysis, optimization, and engineering. It seeks to determine the optimal geometry of structures to minimize performance objectives, subject to physical constraints—often modeled by Partial Differential Equations (PDEs). Traditional approaches commonly assume that these constraints admit a unique solution for each candidate shape, implying a single-valued shape-to-solution map. However, many real-world structures exhibit multistability, where multiple stable configurations exist under identical physical conditions.
This research departs from the single-solution paradigm by investigating shape optimization in the presence of multiple solutions to the same PDE constraints. Focusing on a neo-Hookean hyperelastic model, we formulate an optimization problem aimed at controlling the energy jump between distinct solutions. Drawing on bifurcation theory, we develop a theoretical framework that interprets these solutions as continuous branches parameterized by shape variations. Building on this foundation, we implement a numerical optimization strategy and present numerical results that demonstrate the effectiveness of our approach.