Even if electro-thermal ice protection systems (IPS) consume less energy when operating in de-icing mode than in anti-icing mode, they still need to be optimized for energy usage. The optimization, however, should also take into account the effect of the de-icing system on the aerodynamic performance. The present work offers an optimization framework in which both thermal and aerodynamic viewpoints are taken into account in formulating various objective and constraint functions by considering the energy consumption, the thickness, the volume, the shape and the location of the accreted ice on the surface as the key parameters affecting the energy usage and the aerodynamic performance. The design variables include the power density and the activation time of the electric heating blankets. A derivative-free technique, called the mesh adaptive direct search (MADS) method, is used to carry out the optimization process, which would normally need a large number of unsteady conjugate heat transfer (CHT) calculations for the IPS simulation. To avoid such prohibitive computations, reduced-order modeling (ROM) is used to construct simplified low-dimensional CHT models. The approach is illustrated through several test cases, in which different combinations of objective and constraint functions, design variables and cycling sequence patterns are examined. In these test cases, the energy consumption is significantly reduced compared to the experiments by improving the spatial and temporal distribution of the thermal energy usage. The results show the benefits of the approach in bringing energy, safety and aerodynamic considerations together in designing de-icing systems.