The proposed sequence, although in principle not unique, does make sense.
Clearly the goal is to evaluate averages properties of a cluster of atoms in the vacuum at finite temperature.
Let start clarifying an important point about $NVT$ canonical ensemble. It describes a system of $N$ particles in a volume $V$ at thermal equilibrium with a thermostat at fixed temperature $T$. Notice that what is fixed is the temperature of the thermostat. The system of $N$ particles is not a system at constant kinetic energy.
Therefore, the plot of the temperature of the cluster, i.e. of a quantity proportional to its kinetic energy, must show fluctuations even in the canonical ensemble. And it couldn't be different, since in a classical canonical ensemble one can evaluate analytically the variance of the kinetic energy, finding a non-zero value. The same for the total energy (sum of kinetic and potential energy of the cluster).
Let's come to the proposed protocol. If we want to study properties of a cluster at finite temperature (corresponding for example to a real cluster thermalized with an environment at fixed temperature), the first step is to start with a given atomic configuration adding (and subtracting) energy to the system until the average kinetic energy of the cluster corresponds to the target temperature within a given error bar. This heating step in principle could be accomplished in many ways. For example, an alternative could be a progressive continuous rescaling of velocities.
However, implementing an algorithm to simulate a canonical ensemble at the target temperature $T$ could be a more efficient way to reach quickly typical phase space configurations for a systema at temperature $T$.
After the equilibration phase, the next important step is to simulate the system at thermodynamic equilibrium at fixed temperature. It would be attempting to continue a canonical simulation and in principle nothing is wrong with such an approach. However, efficiency considerations suggest to use a different protocol.
Indeed, just referring to the starting considerations about the canonical ensemble, it turns out that there could be a problem in the case of long production runs to evaluate averages in the canonical ensemble: since the energy of the system must fluctuate, soon or later a rare large fluctuation would happen, making one or more atoms evaporating from the cluster. Nice effect to reproduce the real behavior of the cluster, but a disaster for the computational physicist who would like to collect data for an $N$ particle cluster.
This is the main rationale for using a microcanonical ensemble during the production run.