  We devise a numerical scheme for the time evolution of matrix product
operators by adapting the time-dependent variational principle for matrix
product states [J. Haegeman et al, Phys. Rev. B 94, 165116 (2016)]. A simple
augmentation of the initial operator $\mathcal{O}$ by the Hamiltonian $H$ helps
to conserve the average energy $\mathrm{tr}[H \mathcal{O}(t)]$ in the numerical
scheme and increases the overall precision. As demonstration, we apply the
improved method to a random operator on a small one-dimensional lattice, using
the spin-1 Heisenberg XXZ model Hamiltonian; we observe that the augmentation
reduces the trace-distance to the numerically exact time-evolved operator by a
factor of 10, at the same computational cost.
