Summary: Shape optimization of mechanical devices is investigated in the context of large, geometrically strongly nonlinear deformations and nonlinear hyperelastic constitutive laws. A weighted sum of the structure compliance, its weight, and its surface area are minimized. The resulting nonlinear elastic optimization problem differs significantly from classical shape optimization in linearized elasticity. Indeed, there exist different definitions for the compliance: the change in potential energy of the surface load, the stored elastic deformation energy, and the dissipation associated with the deformation. Furthermore, elastically optimal deformations are no longer unique so that one has to choose the minimizing elastic deformation for which the cost functional should be minimized, and this complicates the mathematical analysis. Additionally, along with the non-uniqueness, buckling instabilities can appear, and the compliance functional may jump as the global equilibrium deformation switches between different bluckling modes. This is associated with a possible non-existence of optimal shapes in a worst-case scenario. In this paper, the sharp-interface description of shapes is relaxed via an Allen-Cahn or Modica-Mortola type phase-field model, and soft material instead of void is considered outside the actual elastic object. An existence result for optimal shapes in the phase field as well as in the sharp-interface model is established, and the model behavior for decreasing phase-field interface width is investigated in terms of \(\Gamma \)-convergence. Computational results are based on a nested optimization with a trust-region method as the inner minimization for the equilibrium deformation and a quasi-Newton method as the outer minimization of the actual objective functional. Furthermore, a multi-scale relaxation approach with respect to the spatial resolution and the phase-field parameter is applied. Various computational studies underline the theoretical observations.