This work aims to address the design optimization of bio-inspired locomotive devices in collective swimming by developing a computational methodology which combines surrogate-based optimization with high fidelity fluid-structure interactions (FSI) simulations of thunniform swimmers. Three main phases highlight the contribution and novelty of the current work. The first phase includes the development and bench-marking of a constrained surrogate-based optimization algorithm which is appropriate to the current design problem. Additionally, new FSI techniques, such as a volume-conservation scheme, has been developed to enhance the accuracy and speed of the simulations. The second phase involves an investigation of the optimized hydrodynamics of a solitary accelerating self-propelled thunniform swimmer during start-up. The third phase extends the analysis to include the optimized hydrodynamics of accelerating swimmers in phalanx schools. Future work includes extending the analysis to the optimized hydrodynamics of steady-state and accelerating swimmers in a diamond-shaped school. The results of the first phase indicate that the proposed optimization algorithm maintains a competitive performance when compared to other gradient-based and gradient-free methods, in dealing with expensive simulations-based black-box optimization problems with constraints. In addition, the proposed optimization algorithm is capable of insuring strictly feasible candidates during the optimization procedure, which is a desirable property in applied engineering problems where design variables must remain feasible for simulations or experiments not to fail. The results of the second phase indicate that the optimized kinematic gait of a solitary accelerating swimmer generates the reverse Karman vortex street associated with high propulsive efficiency. Moreover, the efficiency of sub-optimum modes, in solitary swimming, is found to increase with both the tail amplitude and the effective flapping length of the swimmer, and a new scaling law is proposed to capture these trends. Results of the third phase indicate that the optimal midline kinematics in accelerating phalanx schools resemble those of accelerating solitary swimmers. The optimal separation distance in a phalanx school is shown to be around 2L (where L is the swimmer's total length). Furthermore, separation distance is shown to have a stronger effect, ceteris paribus, on the propulsion efficiency of a school when compared to phase synchronization.