This paper formulates a model order reduction method for electromagnetic boundary element analysis and extends it to computer-aided design integrated shape optimization of multi-frequency electromagnetic scattering problems. Firstly, a series expansion technique is adopted to decouple frequency-dependent terms from the integrands in boundary element formulation, and the second-order Arnoldi procedure is utilized to reduce the order of original systems. Secondly, to ensure geometric exactness and avoid re-meshing during shape optimization, the isogeometric boundary element method is employed, which uses the Non-Uniform Rational B-splines to represent the geometry and to discretize boundary integral equations. Lastly, we employ the Grey Wolf Optimization-Artificial Neural Network as a surrogate model for shape optimization in multi-frequency electromagnetic scattering problems, with the radar cross section as the objective function. Several examples are provided to demonstrate the accuracy and efficiency of the proposed algorithm.