The traditional distribution network reliability assessment uses simulation of component outages to quantify the impact of a set of events on power supply reliability. This simulation-based reliability assessment method requires heuristic or meta-heuristic methods to be used for its inclusion in distribution network expansion planning. To overcome the shortcomings of convergence and local optimization in heuristic or meta-heuristic methods, this paper proposes a mixed integer linear optimization model for multi-stage expansion planning of distribution networks considering reliability costs. This model establishes an explicit assessment model for distribution network reliability index (i.e., the system average interruption duration index, the system average interruption frequency index, and expected energy not supplied), and incorporates its related reliability costs into the multi-stage expansion planning objective function of the distribution network, with system operation constraints as constraints. The effectiveness and applicability of the proposed method was verified on four example systems with 24 nodes, 54 nodes, 86 nodes, and 138 nodes.