Ergonomics is a key factor to consider when designing control architectures for effective physical collaborations between humans and humanoid robots. In contrast, ergonomic indexes are often overlooked in the robot design phase, which leads to suboptimal performance in physical human-robot interaction tasks. This paper proposes a novel methodology for optimizing the design of humanoid robots with respect to ergonomic indicators associated with the interaction of multiple agents. Our approach leverages a dynamic and kinematic parameterization of the robot link and motor specifications to seek for optimal robot designs using a bilevel optimization approach. Specifically, a genetic algorithm first generates robot designs by selecting the link and motor characteristics. Then, we use nonlinear optimization to evaluate interaction ergonomy indexes during collaborative payload lifting with different humans and weights. To assess the effectiveness of our approach, we compare the optimal design obtained using bilevel optimization against the design obtained using nonlinear optimization. Our results show that the proposed approach significantly improves ergonomics in terms of energy expenditure calculated in two reference scenarios involving static and dynamic robot motions. We plan to apply our methodology to drive the design of the ergoCub2 robot, a humanoid intended for optimal physical collaboration with humans in diverse environments