Abstract: Tuned mass dampers (TMDs) effectively reduce structural vibrations. Traditional research on vortex-induced vibration (VIV) control in flexible structures often considers only a single mode and uses basic aerodynamic force models, which may not accurately capture amplitude-dependent damping. However, structures adopted by long-span bridges display multiple modes with frequencies closed each other, and ignoring secondary modes can lead to errors due to their impact on the stiffness and performance of a designing TMD. This study introduces an optimization framework of multiple TMDs adopted in multi-mode VIV control of flexible structures, using a multi-span bridge as a case study. It incorporates secondary modes and nonlinear aerodynamic effects, deriving force-structure-MTMDs system equations in modal coordinates with a polynomial aerodynamic force model. A novel approach is presented to quantify each mode's contribution and to identify essential secondary modes. The results show that the polynomial model accurately predicts the amplitude after TMD installation. Optimal MTMD parameters are determined through global optimization and compared with mode-by-mode design, showing better control performance.