Layout optimization of wind farms constitutes an important and challenging task in complex terrain. This is especially due to the complex interactions of the boundary layer flows in complex terrain and wind turbine wakes, which renders wake modelling in complex terrain difficult. This study tackles this challenge with a new engineering wake model, which is developed by superposing a Gaussian shape wake model on top of the background flow field, assuming that the centerlines of wind turbine wakes follow the streamlines of the background flow field. The model is found to predict wind turbine wakes in complex terrain with good accuracy and at the same time it is computationally cheap to run for optimization applications. Comparisons with high fidelity simulations and field measurements for a real wind farm with 25 turbines in complex terrain demonstrate its effectiveness. A systematic comparison of eight optimization algorithms, which includes two gradient-based and six gradient-free algorithms, is also carried out for the layout optimization problem in complex terrain. To accelerate the optimization process, a double-stage approach is proposed, which optimizes the objective function first neglecting wake effects and then, in the second stage, including them. While all the tested optimization algorithms can improve the original wind farm layout, random search, local search, and pattern search are found to be the top three algorithms in terms of optimization results and computational cost.