Multi-leaf collimators (MLCs) are emerging as the prevalent modality to apply intensity modulated radiotherapy (IMRT). Both the principle and the particular design of MLCs stipulate complex constraints on the practically applicable intensity modulated radiation fields. Most consequentially, the distribution of exposure times across the maximum field outline is either a piecewise constant function in the static mode or a piecewise linear function in the dynamic mode of driving an MLC. In view of clinical utility, the total leaf movement should be minimized, which requires that MLC-related constraints be considered in the dose optimization process. A method is proposed to achieve this for both static MLC fields and dynamic leaf close-in application. The method is an amendment to a generic gradient-based IMRT dose optimization algorithm and solves numerical problems related to the non-convexity of the MLC constraints, which can cause erratic behaviour of a gradient-based algorithm. It employs bistable penalty functions to select preferrable leaf configurations from the configuration space of the MLC, which is limited by specific design features. Together with an 'annealing' escape mechanism from local minima, the algorithm is capable of finding the optimum of an IMRT problem as leaf sequences with minimized leaf travel. In particular, the efficiency of static IMRT can be raised to the levels of unmodulated fields with very few field segments, thereby increasing the utility of IMRT in clinical practice.