The efficiency and scope of aeroelastic wing optimization strategies can be increased using analysis-specific structural idealizations, such as high-fidelity models for detailed stress analyses and low-fidelity models for aeroelastic analyses. In this work, a numerical method is presented that enables efficient and accurate reduction of a high-fidelity finite element model to a Timoshenko beam-based model with lumped masses. The result is a representation based on 13 independent physical beam stiffness parameters per element. The method also yields analytical beam sensitivities with respect to changes in the high-fidelity model. Using these, an approach is suggested for integrating the beam reduction method into a gradient-based multidisciplinary design optimization architecture. The reduction technique is demonstrated on a simplified wing-box and on the University of Bristol Ultra-Green aircraft configuration wing. The effects of unbalanced skin composite laminates, rotated internal ribs, varying wing taper and sweep, and wing boundary constraints on the wing stiffness are shown to be captured with sufficient accuracy for static and dynamic aeroelastic analysis purposes. The accuracy of the analytically calculated gradients is demonstrated by comparison with complex-step and finite difference derivative approximations.