Flank milling is one of the most important technologies for machining of complex surfaces. A small change of the tool orientation in the part coordinate system (PCS) may produce a great rotation of the rotary axes of the machine tool. Therefore, this paper proposes a tool path optimization model for flank milling in the machine coordinate system (MCS). The tool path is computed to smooth the variation of the rotary axes while controlling the geometric deviation. The geometric deviation is measured by the signed distance between the design surface and the tool envelope surface in the PCS. The geometric accuracy is not an objective but a constraint in the proposed optimization model. Given a prescribed geometric tolerance, the tool path smoothness optimization model is reformulated as a constrained nonlinear programming problem. The $\epsilon $ constrained differential evolution with gradient-based mutation ($\epsilon $DEg) is adopted to solve this constrained problem. The validity of the proposed approach is confirmed by numerical examples.