This paper presents a method to significantly accelerate optimization of film cooling systems. The method combines high-fidelity computational fluid dynamics with scalar tracking implemented, a proxy model (linear superposition model) initialized with the computational fluid dynamics solution, and a multi-objective evolutionary algorithm approach. The proposed method is structured as follows: the computational fluid dynamics solution is used to predict the (generally complex) flow domain for the film cooling system; the scalar tracking method identifies the contributions of individual holes to an overall cooling effectiveness distribution by associating a unique passive scalar variable to the flow associated with each hole, and solving an additional advection–diffusion (scalar transport) equation; the proxy model is a (generally linear) superposition model implemented – for example – in Matlab, which inherits the scalar values from the computational fluid dynamics solution, and allows extrapolation of solutions to new design points as part of an optimization process; the optimization process is handled with a multi-objective evolutionary algorithm approach which iterates the proxy model to optimize for a defined objective function. The process works with inner and outer convergence loops. The inner convergence loop is the multi-objective evolutionary algorithm interfacing with the proxy model, which achieves convergence against a design target. At the end of each inner loop cycle, a high-fidelity computational fluid dynamics simulation is run, and this is used to recalibrate the proxy model. Convergence for a given objective function is typically achieved with six outer-loop iterations (high-fidelity computational fluid dynamics runs) and 10,000 inner-loop iterations per outer-loop iteration. The significant advantage of the proposed method is that for certain optimization problems, the computational cost can be reduced by several orders of magnitude, replacing thousands of high-fidelity computational fluid dynamics runs with approximately six computational fluid dynamics runs. The process is demonstrated by applying the optimization method to the film cooling of a flat plate. In our example we have an objective function which maximizes the component life (related to the difference from an arbitrary target temperature distribution) and minimizes the mixing loss introduced by the films. The flow environment was moderately compressible. The optimization converged after six computational fluid dynamics runs. A 30% reduction in mixing loss, a 11% increase in component life, and a 30% reduction in cooling mass flow rate were achieved. The advantages and limitations of the proposed method are also discussed in detail.