As the size of engineered systems grows, problems in reliability theory can become computationally challenging, often due to the combinatorial growth in the number of cut sets. In this paper we demonstrate how Multilevel Monte Carlo (MLMC) — a simulation approach which is typically used for stochastic differential equation models — can be applied in reliability problems by carefully controlling the bias-variance tradeoff in approximating large system behaviour. In this first exposition of MLMC methods in reliability problems we address the canonical problem of estimating the expectation of a functional of system lifetime for non-repairable and repairable components, demonstrating the computational advantages compared to classical Monte Carlo methods. The difference in computational complexity can be orders of magnitude for very large or complicated system structures, or where the desired precision is lower.