Task-invariant control paradigms can enable lowerlimb exoskeletons to provide assistance for their users across various locomotor tasks without prescribing to specific joint kinematics. As an energetic control method, energy shaping can alter a human's body energetics in the closed-loop to provide gait benefits. To obtain the energy shaping law for underactuated systems, a set of nonlinear partial differential equations, called the matching condition, needs to be solved to determine the achievable closed-loop dynamics. However, solving matching conditions for high-dimensional nonlinear systems is generally difficult. In addition, how to define parameters for the closed-loop dynamics that render the optimal exoskeleton assistance remains unclear. In this paper, we proposed a twolayer, human-in-the-loop optimization framework for lower-limb exoskeletons to customize their assistance to human users. The inner-layer optimization finds solutions to the matching condition, meanwhile following the energy trajectories of a virtual reference model defined based on the self-selected gaits of humans and a scaled version of their anatomical parameters. The outer-layer incorporates human-in-the-loop Bayesian Optimization to update reference energy's parameters for reducing metabolic costs. Simulation results on two biped models demonstrate that the proposed framework can solve matching conditions numerically at the selected timestamps and the associated energy shaping strategies can reduce human metabolic cost. Moreover, exoskeletons torques calculated using an able-bodied subject's kinematic data well match human biological torques.