To increase airspace capacity, alleviate flight delay, and improve network robustness, an optimization method of multi-layer air transportation networks is put forward based on Laplacian energy maximization. The effectiveness of taking Laplacian energy as a measure of network robustness is validated through numerical experiments. The flight routes addition optimization model is proposed with the principle of maximizing Laplacian energy. Three methods including the depth-first search (DFS) algorithm, greedy algorithm and Monte-Carlo tree search (MCTS) algorithm are applied to solve the proposed problem. The trade-off between system performance and computational efficiency is compared through simulation experiments. Finally, a case study on Chinese airport network (CAN) is conducted using the proposed model. Through encapsulating it into multi-layer infrastructure via k-core decomposition algorithm, Laplacian energy maximization for the sub-networks is discussed which can provide a useful tool for the decision-makers to optimize the robustness of the air transportation network on different scales.