In the present paper, a generalization of the method of partial summation of the expansion of the thermodynamical potential is proposed. This generalization allows one to obtain the corresponding equations for higher-order correlation matrices, as well as to formulate the variational method for their solution. We show that correlation matrices of equilibrium quantum system satisfy a variational principle for thermodynamic potential which is functional of these matrices that provides a thermodynamic consistency of the theory. This result is similar to a variational principle for correlation functions of classical systems.