Methods of quantum nuclear wave-function dynamics have become very efficient in simulating large isolated systems using the time-dependent variational principle (TDVP). However, a straightforward extension of the TDVP to the density matrix framework gives rise to methods that do not conserve the energy in the isolated system limit and the total system population for open systems where only energy exchange with environment is allowed. These problems arise when the system density is in a mixed state and is simulated using an incomplete basis. Thus, the basis set incompleteness, which is inevitable in practical calculations, creates artificial channels for energy and population dissipation. To overcome this unphysical behavior, we have introduced a constrained Lagrangian formulation of TDVP applied to a non-stochastic open system Schrödinger equation [L. Joubert-Doriol, I. G. Ryabinkin, and A. F. Izmaylov, J. Chem. Phys. 141, 234112 (2014)]. While our formulation can be applied to any variational ansatz for the system density matrix, derivation of working equations and numerical assessment is done within the variational multiconfiguration Gaussian approach for a two-dimensional linear vibronic coupling model system interacting with a harmonic bath.