Even if electro-thermal ice protection systems consume less energy operating in de-icing mode than in anti-icing mode, they still need to be optimized for energy usage. This paper defines suitable design variables, objective functions and constraints for such optimization. A derivative-free technique is used to carry out the optimization process, which normally needs a large number of unsteady conjugate heat transfer simulations. To avoid such prohibitive computations, reduced-order modeling (ROM) is used to construct a simplified low-dimensional model at each time step of the simulation. The approach is illustrated through test cases, showing its benefits in bringing energy and safety considerations together in designing de-icing systems. Nomenclature ice C = specific heat capacity of ice w C = specific heat capacity of water E = power density cycle E = total heating energy c h = convective heat transfer coefficient f h = water film thickness k = dimension of design space evap L = latent heat of evaporation fusion L = latent heat of fusion subl L = latent heat of sublimation LWC = liquid water content evap m = mass transfer rate by evaporation ice m = mass transfer rate by ice accretion MVD = mean volumetric diameter = POD basis function = POD basis coefficient M = number of modes used by POD N = number of snapshots IPS Q = Conductive heat transfer rate by IPS T = surface temperature target T = target surface temperature , d T = far-field droplet temperature max,ice T = maximum ice thickness max, , ice cycle T = maximum ice thickness in one cycle d u = droplets velocity f u = mean water film velocity V = snapshot solution U = free-stream velocity = solid emissivity w = density of water