This paper presents an improved Monte Carlo algorithm for calculating thermoluminescence (TL) curves for complex cluster systems, which contain a large number of localized levels. The new algorithm is based on generation of two random numbers, instead of the many random number generations used in previous approaches. The first random number generation is associated with the choice of which type of energy transition takes place in the model, while the second random number generation is connected with the choice of a cluster group for which this transition occurred. The proposed new algorithm was tested on a new TL model consisting of clusters containing both electron traps and deep competitor traps. The luminescence centers in the model are treated as uniformly distributed defects. Applicability of the Monte Carlo method is demonstrated by simulating for the first time in this type of model the complete TL process, including the excitation, relaxation and heating stages. The new model explains the anomalous heating rate effect of TL, and describes the effect of localized trapping processes within the cluster on the intensity and temperature position of the TL peaks.